Leaf Area Index identified as a major source of variability in modelled CO 2 fertilization

Leaf Area Index identified as a major source of variability in modelled CO2 fertilization Qianyu Li1,4, Xingjie Lu2, Yingping Wang3, Xin Huang2, Peter M. Cox4, Yiqi Luo1,2 1Ministry of Education Key Laboratory for Earth System Modeling, Department of Earth System Science, Tsinghua University, Beijing 100084, China 5 2 Center for Ecosystem Science and Society (Ecoss), Northern Arizona University, Flagstaff, AZ 86011, USA 3 CSIRO Oceans and Atmosphere, PMB #1, Aspendale, Victoria 3195, Australia 4College of Engineering, Mathematics and Physical Sciences, University of Exeter, Exeter, EX4 4QF, UK


Introduction
Terrestrial ecosystems absorb roughly 30 % of anthropogenic CO 2 emissions, and is of great uncertainty and vulnerable to global climate change (Cox et al., 2000;Le Quéré et al., 2018).Persistent increase in atmospheric CO 2 concentration will stimulate plant growth and ecosystem carbon storage, forming a negative feedback to CO 2 concentration (Long et al., 2004;Friedlingstein et al., 2006;Canadell et al., 2007).This concentration-carbon feedback (β), also called the CO 2 fertilization effect, has been identified as a major uncertainty in modeling terrestrial carbon-cycle response to historical climate change (Huntzinger et al., 2017).In the Coupled Model Intercomparison Project (C4MIP) and the Coupled Model Intercomparison Project Phase 5 (CMIP5), all models agree that terrestrial carbon sink will gradually saturate in the future but disagree on the magnitude of β (Friedlingstein et al., 2006(Friedlingstein et al., , 2015;;Arora et al., 2013).Some studies pointed out that the contribution of β is 4 to 4.5 times larger, and more uncertain, than climate-carbon feedback (γ ) (Gregory et al., 2009;Bonan and Levis, 2010;Arora et al., 2013).Apart from the substantial uncertainty across different mod-Published by Copernicus Publications on behalf of the European Geosciences Union.els, Smith et al. (2016) suggested that Earth system models (ESMs) in CMIP5 overestimate global terrestrial β values compared with remote sensing data and Free-Air CO 2 Enrichment (FACE) experimental results.Though satellite products they used may underestimate the effect of CO 2 fertilization on net primary productivity (De Kauwe et al., 2016), the large disparity between models and FACE experiments gives us little confidence in making policies to combat global warming.
The response of the ecosystem carbon cycle to elevated CO 2 (eCO 2 ) is primarily driven by stimulation of leaf-level carboxylation rate in plants (Polglase and Wang, 1992;Long et al., 2004;Heimann et al., 2008).The CO 2 stimulation of carboxylation then translates into increasing gross primary production (GPP) and net primary production (NPP), possibly leading to increased biomass and soil carbon storage and slowing down anthropogenically driven increase in atmospheric CO 2 (Canadell et al., 2007;Iversen et al., 2012).The leaf-level CO 2 fertilization for C 3 plants is generally well characterized with models from Farquhar et al. (1980), and the basic biochemical mechanisms have been adopted by most land-surface models, although some models implement variants of Farquhar et al. (1980) (Rogers et al., 2016).Previous research with both theoretical analysis and data synthesis from a large number of experiments has revealed that normalized CO 2 sensitivity of leaf-level photosynthesis, which represents kinetics sensitivity of photosynthetic enzymes, varies little among different C 3 species at a given CO 2 concentration (Luo and Mooney, 1996;Luo et al., 1996).However, the CO 2 fertilization effects are considerably more variable at canopy and ecosystem level than at the leaf level, because a cascade of uncertain processes, such as soil moisture feedback (Fatichi et al., 2016), canopy scaling (Rogers et al., 2016), nutrient limitation (Zaehle et al., 2014), allocation (De Kauwe et al., 2014), and carbon turnover process (Friend et al., 2014), influences the responses of GPP, NPP and carbon storage.Therefore, understanding which processes in ecosystem models amplify the variability in β from biochemical and leaf levels to canopy and ecosystem levels is quite important.
Leaf area index (LAI) largely affects canopy assimilation and plant growth under eCO 2 .Many satellite products exhibit increasing trends of LAI over the past 30 years, although marked disparity still exists among these products (Jiang et al., 2017).Zhu et al. (2016) have attributed global increases in satellite LAI primarily to increased CO 2 concentration.LAI plays a key role in scaling leaf-level biogeophysical and biogeochemical processes to global-scale responses in ecosystem models, and the representation of LAI in models causes large uncertainty (Ewert, 2004;Hasegawa et al., 2017).Models generally predict that LAI dynamics will respond to eCO 2 positively due to enhanced NPP and leaf biomass (De Kauwe et al., 2014).But how the increasing LAI in turn feeds back to ecosystem carbon uptake as a result of more light interception has not been discussed in previous research.The relative contributions of the leaf-level photosynthesis and LAI to modeled β have rarely been quantified and compared.
The CO 2 fertilization effects depend on locations, vegetation types and soil nutrient conditions.The strongest absolute CO 2 fertilization effect has been found in tropical and temperate forests where the larger biomass presents than other regions.In comparison, the weakest response to eCO 2 occurs in boreal forests (Joos et al., 2001;Peng et al., 2014).But with gradual eCO 2 , relative response in tropical forests might not be very high owing to light limitation caused by canopy closure (Norby et al., 2005).In addition, β might be overestimated by the neglect of nitrogen (N) limitations on plant growth (Luo et al., 2004;Thornton et al., 2009;Coskun et al., 2016).Several lines of evidence suggest that N availability also influences decomposition of soil organic matter (Hunt et al., 1988;Neff et al., 2002;Averill et al., 2016).β will be reduced by 50 %-78 % in C-N coupled simulations compared with C-only simulations in land-surface models (Thornton et al., 2007;Sokolov et al., 2008;Zaehle et al., 2010).Inadequate phosphorus (P) will also constrain terrestrial carbon uptake, especially in tropical areas (Aerts and Chapin, 2000;Vitousek et al., 2010).It is reported that N limitation on carbon uptake is significant in boreal ecosystems, while P limitation has a profound influence in tropical ecosystems in the CASA-CNP model (Wang et al., 2010).However, whether N and P limitations affect the variability of β across different vegetation types at different hierarchical levels from biochemistry to ecosystem carbon storage has not been carefully examined.
In this study, we tried to answer the following questions: how does variability, as measured by the coefficient of variation (CV) within and across different plant functional types (PFTs), in the CO 2 fertilization effect change at different hierarchical levels from leaf to canopy GPP, ecosystem NPP and total carbon storage levels?What is the most important process causing the variability of β for different geographical locations and PFTs?How do nutrient limitations influence the variability of β at different hierarchical levels?We used the Community Atmosphere Biosphere Land Exchange model (CABLE) to identify key mechanisms driving diverse β values under the RCP 8.5 scenario.

CABLE model description
CABLE (version 2.0) is the Australian community landsurface model (Kowalczyk et al., 2006) and incorporates CASA-CNP to simulate global carbon (C), nitrogen (N) and phosphorus (P) cycles (Wang et al., 2010(Wang et al., , 2011)).Leaf photosynthesis, stomatal conductance, and heat and water transfer in CABLE are calculated using the two-leaf approach (Wang and Leuning, 1998) for both sunlit leaves and shaded leaves.
The descriptions of the photosynthesis module are in Supplement S1.
Leaf area index (LAI) is calculated as where C leaf is leaf carbon pool and SLA is specific leaf area.
In the CABLE model, leaf growth is divided into four phases.Phase 1 is from leaf budburst to the beginning of steady leaf growth, phase 2 is from the start of steady leaf growth to the start of leaf senescence, phase 3 is the period of leaf senescence, and phase 4 is from the end of leaf senescence to the start of leaf bud burst.During phase 1, allocation of available carbon to leaf is fixed to 0.8, and allocations to wood and root are set to 0.1 for woody biomes and 0 and 0.2, respectively, for non-woody biomes.During steady leaf growth (phase 2), the allocation coefficients are constants but vary from biome to biome, taking their values from Fung et al. (2005).During phases 3 and 4, the leaf allocation is zero and available carbon is divided between wood and root in proportion to their allocation coefficients.For evergreen biomes, leaf phenology remains at phase 2 throughout the year (Wang et al., 2010).SLA is PFT-specific and does not change through time in this study.
GPP is the sum of canopy net photosynthesis rate (A) and day respiration (R d ).NPP is calculated as the difference between GPP and autotrophic respiration (R a ) (including maintenance and growth respiration), and acts as an input to the compartmental nine-pool carbon-cycle model.The network for carbon transfer in the compartmental model is based on the CASA' model (Fung et al., 2005), including three vegetation pools (leaf, wood and root), three litter pools (metabolic litter, structure litter and coarse wood debris), and three soil pools (fast soil pool, slow soil pool and passive soil pools).Heterotrophic soil respiration (R h ) is calculated as the sum of the respired CO 2 from the decomposition of all litter and soil organic carbon pools (Wang et al., 2010).Wang et al. (2012) and Zhang et al. (2013) provided details explaining how nutrient limitations are incorporated into the carbon cycle in the CASA-CNP module in the CABLE model.In brief, NPP is calculated as where LAI represents leaf area index, and V cmax and J max are maximum carboxylation rate and maximum rate of electron transport of the top leaves, respectively; both are linearly dependent on leaf N (g N m −2 ) according to the relationships developed by Kattge et al. (2009) for different plant functional types.R mi is maintenance respiration rate of plant tissue (i: leaf, wood and root), contingent on nitrogen amount in each part of the plant.R g is growth respiration, which is described as a function of leaf nitrogen to phosphorus ratio.Heterotrophic respiration (R h ) is limited by the mineral N pool required for microbial soil carbon decomposition (Wang et al., 2010).Net ecosystem productivity (NEP = GPP − R a − R h ) is the amount of carbon that is either sequestered or lost from ecosystems, and is controlled by N and P availability via abovementioned C-N-P interactions.

Experimental design
CABLE was run from 1901 to 2100 for C-only, C-N and C-N-P modes.C-only simulation was designed to identify the key carbon-cycle processes that influence the variability of the CO 2 fertilization effects.C-N and C-N-P simulations were run to explore how nutrients affect the patterns of and mechanisms underlying the variability of the CO 2 fertilization effects.The respective effects of N and P can be calculated through the difference in the carbon uptake between C-N and C-only or C-N-P and C-N simulations.CABLE was first spun up by using meteorological forcing from Community Climate System Model (CCSM) simulations (Hurrell et al., 2013) during 1901 to 1910 until steady states were achieved for the C-only, C-N and C-N-P cases separately.
Hourly meteorological driving data include temperature, specific humidity, air pressure, downward solar radiation, downward long-wave radiation, rainfall, snowfall, and wind.In order to separate the CO 2 fertilization effect from the effect of climate change, climate forcing was held as the average annual cycle of CCSM meteorological data from 1901 to 2100.Atmospheric CO 2 concentrations from 1901 to 2100 were taken from the CMIP5 dataset, representing global annual averages and the RCP 8.5 scenario after 2010 (Etheridge et al., 1996;MacFarling Meure et al., 2006).The spatial resolution of CABLE used here is 1.9 • × 2.5 • (latitude vs. longitude).N deposition is prescribed from atmospheric transport models (Lamarque et al., 2010(Lamarque et al., , 2011)), spatially explicit but fixed as the average from 1901 to 2100 in time.N fixation is prescribed from a process-based model, spatially explicit but constant in time (Wang and Houlton, 2009).P enters ecosystems through constant rates of weathering and atmospheric deposition (from Mahowald et al., 2008).

Calculation of β values at five hierarchical levels
We aimed to analyze the CO 2 fertilization effects for biochemical reaction (L), leaf photosynthesis rate (p), leaf-tocanopy scaling factor (S), leaf area index (LAI), sunlit leaf GPP (GPP sun ), shaded leaf GPP (GPP sha ), canopy GPP, NPP, and ecosystem carbon storage (cpool) from C-only, C-N and C-N-P simulations of CABLE.Canopy GPP is the sum of sunlit leaf GPP and shaded leaf GPP.Ecosystem carbon storage is the sum of plant, litter and soil carbon stock.Since CO 2 concentration increases on a yearly basis, annual carbon fluxes and storages such as GPP sun , GPP sha , canopy GPP, NPP and ecosystem carbon storage were calculated.Leaf-tocanopy scaling factor and LAI were averaged within a year.β values of these variables were calculated as the normalized sensitivities of those variables to atmospheric CO 2 concentration (C a ) as β V : where V in the denominator represents the average annual value of S sun , S sha , LAI, GPP, GPP sun , GPP sha , NPP and ecosystem carbon storage between 2 consecutive years.Subscripts "sun" and "sha" denote the sunlit and shaded components.dV is the difference of these variables between 2 consecutive years.dC a is the difference of the corresponding C a .
The unit of β V is ppm −1 .It should be noted that β V is the relative response, which is similar to the traditional definition of the β factor by Bacastow and Keeling (1973), but different from the carbon-concentration feedback parameter in Friedlingstein et al. (2006).The relative response facilitates the comparison among PFTs with different initial biomass and the comparison across carbon fluxes and storages with different units.
Leaf biochemical response (L) was first proposed by Luo et al. (1996).The L function is the normalized response of leaf photosynthesis rate to a small change in intercellular CO 2 concentration (C i ) and has been suggested to be an invariant function for C 3 plants grown in diverse environments.The rate of photosynthesis is typically RuBP-regenerationlimited under a high CO 2 concentration.We found photosynthesis rates are increasingly limited by RuBP regeneration under the RCP 8.5 scenario.Besides, theoretical analysis by Luo and Mooney (1996) showed that biochemical responses are similar for either Rubisco-or RuBP-limited photosynthesis.In this study, L can be used to indicate leaf biochemical response to eCO 2 .For sunlit leaf and shaded leaf, formulations of L under RuBP-regeneration limitation are defined as In this study, * sun and * sha are yearly average CO 2 compensation points in the absence of day respiration for sunlit leaf and shaded leaf, respectively.C i varies significantly on sub-daily, intra-annual and inter-annual bases.We are interested in how C i responds to eCO 2 on an inter-annual basis.So, we first outputted hourly C i and then calculated yearly GPP-weighted average C i for sunlit leaf (C isun ) and shaded leaf (C isha ).
Then leaf-level β p is defined as the product of L and dC i dC a .For sunlit leaf and shaded leaf, the formulations are Leaf-to-canopy scaling factor (S) scales fluxes at the single top leaf of the canopy to whole canopy fluxes.The formulations of S for sunlit leaves and shaded leaves are where k b is the extinction coefficient of a canopy of black leaves for direct beam radiation.k n is an empirical parameter used to describe the vertical distribution of leaf nitrogen in the canopy (Kowalczyk et al., 2006).In our simulation, k n is uniformly assigned as 0.001 for different PFTs.
The leaf-to-canopy scaling factor varies with time because k b is the function of Sun angle, and LAI varies seasonally and inter-annually.The annual value of the leaf-to-canopy scaling factor was just calculated as the average of hourly leaf-to-canopy scaling factors in a year.Big-leaf β GPP sun (or β GPP sha ) can be decomposed as the sum of normalized sensitivity of photosynthesis rate, β p sun (or β p sha ), and leaf-to-canopy scaling factor, β S sun (or β S sha ), as shown in Eqs. ( 10) and ( 11).Detailed mathematical derivations are in Supplement S2.
There are 10 patches in each model grid in CABLE.Each patch consists of a certain land use type with a specific fraction.To study the variation of β across different C 3 PFTs, biome-level parameters such as * sun , C isun , S sun and LAI were calculated as mean values based on PFTs, whereas biome-level GPP, GPP sun , GPP sha , NPP and ecosystem carbon storage were integrated sums based on PFTs.Then L sun , L sha , β p sun , β p sha , β GPP , β NPP and β cpool at the year 2023 (relative to 2022) for different C 3 PFTs were calculated and compared.Coefficients of variation (CVs) of β values were calculated across various C 3 PFTs for these hierarchical levels.
The year 2023 was chosen because large oscillations of LAI occurred for shrub after 2025 in the C-N-P simulation (Supplement Fig. S1c).For C-N and C-N-P simulations, the time series of LAI, GPP, and NPP for shrub, C 3 grass and tundra underwent small short-term variability and therefore were smoothed using the "smooth" function in MATLAB software before the calculation of β.We also calculated β values for each patch and CV of β values across different geographical locations within a specific PFT at different hierarchical levels at the year 2023 to explore the variability of β within the same PFTs.All abovementioned calculations were processed in MATLAB R2014b.

Temporal trends of β at ecosystem level for different PFTs
In C-only simulation, β cpool values for different C 3 PFTs all decline with time from 2011 to 2100 under the RCP 8.5 scenario (Fig. 1a).However, the magnitudes of  1 and Fig. S2).
C i /C a values for both sunlit and shaded leaves in C-N-P simulation are very similar to those in C-N simulation (Table 1 and Fig. S3).In all of the simulations, values of CO 2 compensation point in the absence of day respiration ( * ) for a specific PFT do not change over time since air temperature as an input to the model is not affected by the biophysical feedback in the offline model simulations (Figs. 2c,d,S2c,d,S3c,d).But there is a huge variance of * across different C 3 PFTs because of different leaf temperature which * values depend on.

Comparison of β at different hierarchical levels
To further trace the cause of the divergence of β across PFTs as shown in Fig. 1 at a specific time, L sun , L sha , β p sun , β p sha , β GPP , β NPP and β cpool at the year 2023 for different C 3 PFTs in all simulations were plotted in Fig. 3. CV is marked above data points for each variable to indicate degree of variation across different C 3 PFTs.In C-only simulation (Fig. 3a), results show that at leaf biochemical level, L values for sunlit leaf and shaded leaf range from 0.00055 to 0.00097 ppm −1 .Variations of L sun and L sha among PFTs are small (CV = 0.15 and 0.13).At leaf photosynthesis level, β p sun and β p sha for the seven PFTs vary from 0.00041 to 0.00072 ppm −1 , and the variations among different PFTs are not significant (CV = 0.18 and 0.12).But β values are diverging when scaled up to GPP level, with CV jumping to 0.49 among PFTs.β values of deciduous broadleaf forest and shrub greatly increase from leaf level to GPP level.However, canopy scaling effects do not significantly amplify β values at canopy levels (β GPP ) for deciduous needleleaf forest, tundra and evergreen broadleaf forest.Magnitudes and variance of β NPP are similar to those of β GPP because NPP linearly correlates with GPP for all C 3 PFTs (Fig. S4).
Magnitudes of β cpool for all PFTs are decreased compared with those of β NPP and β GPP .Deciduous broadleaf forest and shrub have the highest β GPP and β NPP values (around 0.0026 ppm −1 ).Deciduous broadleaf forest has the greatest β cpool value (around 0.0018 ppm −1 ) among all.Deciduous needleleaf forest has the lowest β GPP , β NPP and β cpool values.CV of β cpool among different PFTs reaches the highest (0.58) compared with CV of β values at other levels.
In C-N and C-N-P simulations, magnitudes and variations of β at leaf biochemical and photosynthetic levels are comparable to those in C-only simulation because C i and * values only slightly change under nutrient limitations (Figs.3b, c, S2, S3).Nutrient-limited β GPP values are smaller than those in C-only simulation, except for evergreen broadleaf forest.There is a large divergence of nutrientlimited β GPP across different PFTs, which is similar to Conly simulation.However, unlike in C-only simulation, β NPP values in nutrient-coupled simulations are reduced for most C 3 PFTs and diverge more compared with β GPP values.Coefficients of variation (CVs) of β cpool in nutrient-coupled simulations exceed 0.8, larger than that in C-only simulation.
Within-PFT variations of β in C-only simulation were listed in Table 2, including CVs for biochemical response L, leaf-level β p , β GPP , β NPP and β cpool across different geographical locations within each PFT.Variations of biochemical and leaf-level responses are relatively smaller than those at canopy and ecosystem levels within all C 3 PFTs.β GPP values are greatly differentiated across different geographical locations.Variations of β NPP are very similar to those of β GPP within all PFTs except the evergreen needleleaf forest.CVs of β cpool are lower than those of β NPP within most PFTs except evergreen broadleaf forest and tundra.Within-PFT variations of β in C-N and C-N-P simulations are similar to those in C-only simulation (data not shown).
To further explore why β values at canopy and ecosystem levels are diverging across different C 3 PFTs, the correlations between β GPP and β LAI , β NPP and β LAI , β cpool and β LAI for C-only, C-N and C-N-P simulations were plotted at the year 2023.Results show that β GPP , β NPP and β cpool all have significant linear correlations with β LAI across different C 3 PFTs (Fig. 4).Results also show that β LAI linearly correlates with β GPP , β NPP and β cpool across patches within the same PFT, although there are some discontinuous points within evergreen broadleaf forest where the canopy of many patches closes (Figs.S5-S7).Therefore variations of β values from leaf to ecosystem scale can be well explained by β LAI or the LAI response to increasing CO 2 .

β of sunlit and shaded leaves
To understand the in-depth mechanism for the influence of LAI on canopy GPP, we investigated the response of sunlit and shaded leaf GPP separately from C-only simulation.Temporal trends of sunlit leaf GPP (GPP sun ) and shaded leaf GPP (GPP sha ) were plotted for each type of C 3 PFTs from  1901 to 2100 in Fig. 5. From the beginning of the simulation, GPP sha is higher than GPP sun for all C 3 PFTs.With significant increases in CO 2 concentration from 2011, GPP sha responds more drastically than GPP sun .Shaded leaf GPP of deciduous broadleaf forest and shrub responds to eCO 2 more significantly than other PFTs.However, a single sunlit leaf has a higher photosynthesis rate (p sun ) than a shaded leaf (p sha ) because of more radiation absorbed.Thus, the LAIdependent canopy scaling factor of shaded leaves (S sha ) contributes more to the magnitude and sensitivity of canopy GPP than photosynthesis rate.Then temporal trends were plotted for β GPP sun (β GPP sha ) and decomposing factors β p sun (β p sha ) and β S sun (β S sha ) for each PFT as Eqs.( 10) and ( 11) to further evaluate the above inference.Results show that both of the sensitivities of GPP sun and GPP sha tend to approach zero through time because the decomposing factors β p sun , β p sha , β S sun and β S sha all decline with time (Fig. 6).β p sun and β p sha overlap through time for each PFT.Magnitudes of β GPP sha are higher than those of β GPP sun for all C 3 PFTs because the values of β S sha are higher than those of β S sun .For deciduous needleleaf forest and tundra, both β p sun (β p sha ) and β S sun (β S sha ) contribute to the magnitudes and trends of β GPP sun (β GPP sha ).For evergreen needleleaf forest, deciduous broadleaf forest, shrub and C 3 grass, β S sun (β S sha ) dominates the magnitude and change in β GPP sun (β GPP sha ).For evergreen broadleaf forest, β S sha predominates the magnitude and change in β GPP sha before 2035.

Variations of biochemical and photosynthetic responses to eCO 2
The direct CO 2 fertilization effect occurs at leaf level and is determined by kinetic sensitivity of Rubisco enzymes to internal leaf CO 2 concentration.In fact, the normalized short-term sensitivity of leaf-level photosynthesis to CO 2 is mainly regulated by C i and slightly influenced by leaf temperature, regardless of light, nutrient availability, and species characteristics (Luo et al., 1996;Luo and Mooney, 1996).In our study, the modeled C i /C a ratio is approximately constant with eCO 2 for a specific PFT, and varies little within and across PFTs in all simulations.This is in line with FACE experimental results which show almost constant C i /C a values for different PFTs under CO 2 fertilization (Drake et al., 1997;Long et al., 2004).* varies little for different species and only depends on leaf temperature (Luo and Mooney, 1996).Sensitivity analysis in a previous study has shown that a ±5 • of leaf temperature changes caused approximately ±7 ppm changes in * , leading to variation of 0.12 to leaflevel β (Luo and Mooney, 1996).The overall variation of leaf-level β caused by variation in leaf temperature is still quite small compared with that of β GPP .Therefore, biochemical and leaf-level β values vary little within and among PFTs in this study.Our results also illustrate that nutrient effects do not significantly change C i and * , leading to similar biochemical and leaf-level β values in all simulations, which is in accordance with Luo et al. (1996).To identify the source of uncertainty of β in CMIP5 models, Hajima et al. ( 2014) decomposed β into several carbon-cycle components.They used GPP divided by LAI (GPP / LAI) as a proxy to represent leaf-level photosynthesis for CMIP5 models, since there are no leaf-level process outputs of these models.They found the sensitivities of GPP / LAI to eCO 2 diverged a lot among models.One possible issue of this calculation is that it ignores different canopy structures used by each CMIP5 model, such as bigleaf, two-leaf or multiple-layer.Our results just show that the sensitivities of GPP / LAI are different from our mechanistic calculation of leaf-level β for different PFTs in a two-leaf model.β values estimated from GPP / LAI formulation are greatly underestimated for woody trees and slightly overestimated for C 3 grass and tundra, but best match for shrub if compared with our calculation (Fig. S8).Therefore diagnostics such as C i and * for leaf-level β are more desirable for woody trees.Another advantage of our calculation of leaflevel β is that the reason for the divergence of leaf-level β across PFTs can be traced back to the difference from C i and leaf temperature as shown in Fig. 2.

Variations of β at canopy and ecosystem levels
The two-leaf scaling scheme in CABLE is widely employed by many land-surface models, such as Community Land Model version 4.5 (CLM4.5,Oleson et al., 2013) and the Joint UK Land Environment Simulator version 4.5 (JULES4.5,Best et al., 2011;Clark et al., 2011;Harper et al., 2016).We found the responses of ecosystem carbon cycle to eCO 2 diverge primarily because the responses of LAI diverge within and among PFTs in all simulations.Besides, GPP of shaded leaves responds to eCO 2 more strongly than GPP of sunlit leaves for all C 3 PFTs.This is because the portion of shaded leaves increases exponentially with increasing LAI (Fig. S9), leading to a rapid change in shaded leaf GPP, while for sunlit leaves, GPP shows a saturating response because of the decreasing portion of sunlit leaves with increasing LAI (Dai et al., 2004).Our results also indicate that saturation of GPP is not only regulated by the leaf-level photosynthetic response, but also by the response of the LAI-dependent scaling factor to eCO 2 .For shaded leaves, the sensitivity of the LAI-dependent scaling factor contributes more to the magnitude and trend of β GPP sha than that of photosynthesis rate.GPP sha is higher than GPP sun for all PFTs.With significant increase in CO 2 concentration from 2011, GPP sha responds more drastically than GPP sun .Abbreviations are the same as Fig. 1.
The evidence all suggests LAI is a key process in modeling the response of ecosystem carbon cycle to climate change.
It has been reported that different CMIP5 models have simulated diverse LAI during 1985-2006.And modeled LAI values in most CMIP5 models have been overestimated according to satellite products (Anav et al., 2013).Many global vegetation models simulated increasing LAI trends globally in response to eCO 2 during the historical period (Zhu et al., 2016).Our modeling study also shows that LAI responds positively to eCO 2 for all C 3 PFTs in all simulations.But experimental results are not consistent.In one review paper with 12 FACE experimental results, trees had a 21 % increase in LAI and herbaceous C 3 grasses did not show a significant change in LAI (Ainsworth and Long, 2005).Some studies reported that LAI dynamics did not significantly change in specific FACE experiments, such as in a closed-canopy deciduous broadleaf forest (ORNL FACE, Norby et al., 2003) and in a mature evergreen broadleaf forest (EucFACE, Duursma et al., 2016).The negligible change in LAI at the Euc-FACE probably leads to an insignificant response of productivity at this site, even though leaf photosynthesis rate significantly increases under eCO 2 (Ellsworth et al., 2017).Besides the impact of LAI on the global carbon cycle, the increasing trend of LAI exerts profound biophysical impacts on climate by altering the energy and water cycles on the Earth's surface (Forzieri et al., 2017;Zeng et al., 2017).But there is a great uncertainty in the relationships between LAI and biophysical processes among land-surface models (Forzieri et al., 2018).
In this study, modeled nutrient-unlimited β GPP and β NPP values are higher than leaf photosynthetic responses for all C 3 PFTs in C-only simulation (Fig. 3a).Nutrient-limited β NPP are still higher than photosynthetic responses for many PFTs in C-N and C-N-P simulations (Fig. 3b, c).However, it is generally observed in experiments that the leaf-level response is consistently larger than the whole plant response (Long et al., 2006;Leuzinger et al., 2011).One possible reason is that models overestimate the response of LAI to eCO 2 , as this study has shown that LAI is an important factor in driving ecosystem response to CO 2 fertilization.And it is also likely the overestimation of the response of LAI to eCO 2 is responsible for the overestimation of CO 2 fertilization in ESMs reported by previous studies (Smith et al., 2015;Mystakidis et al., 2017).
The overall response of LAI to eCO 2 depends on several processes in this study: ( 1  carbon is allocated to woods or roots at higher CO 2 concentration (De Kauwe et al., 2014).Unfortunately, CABLE has fixed allocation coefficients and likely overestimates LAI response, leading to overestimated responses of GPP, NPP and total carbon storage.Third, we fixed SLA to calculate LAI in CABLE.But a reduction in SLA is a commonly observed response in eCO 2 experiments (Luo et al., 1994;Ainsworth and Long, 2005;De Kauwe et al., 2014).Tachiiri et al. (2012) also found SLA and β values are most effectively constrained by observed LAI to smaller values in a model.Therefore, the fixed SLA may also lead to over-prediction of the response of canopy cover to eCO 2 .Fourth, in our results, LAI values for most C 3 PFTs are below the maximum LAI limits with eCO 2 in C-only simulation.With only one exception, LAI values of many evergreen broadleaf forest patches saturate at the prescribed maximum value under high CO 2 concentration (Fig. S1a and Table S1).That is why the sensitivity of LAI for evergreen broadleaf forest is low and thus leads to small relative GPP enhancements.If the preset LAI upper limits are narrowed, β values are expected to be significantly reduced.Hence model parameters related to LAI need to be better calibrated according to experiments and observations in order to better represent the response of ecosystem productivity to eCO 2 (De Kauwe et al., 2014;Qu and Zhuang, 2018).
In this study, the almost identical values and variance of β NPP to those of β GPP within and across C 3 PFTs in Conly simulation suggest carbon use efficiency (CUE) does not change with eCO 2 , as autotrophic respiration is calculated from GPP and plant carbon.In C-N and C-N-P simulations, magnitudes of β NPP for all C 3 PFTs except evergreen broadleaf forest all decline compared with those of β GPP , indicating CUE also declines with eCO 2 under nutrient limitations.However, FACE experimental results indicate that CUE values under eCO 2 are not changed in the N-limited Duke site (Hamilton et al., 2002;Schäfer et al., 2003), increase in the fertile POPFACE site (Gielen et al., 2005) or decrease in the fertile ORNL site (DeLucia et al., 2005).Thus, representations of nutrient effects on GPP and autotrophic respiration in land-surface models should be carefully calibrated with experimental data (DeLucia et al., 2007).Our results also show that β NPP values diverge more than β GPP values across different PFTs in nutrient-coupled simulations, because the different nutrient-limiting effects on autotrophic respiration introduce additional variation across different PFTs.Although β values at ecosystem levels are more variable with nutrient effects, LAI responses are still linearly correlated well with β GPP , β NPP and β cpool across C 3 PFTs in nutrient-coupled simulations as in C-only simulation, confirming the dominant role of LAI in regulating carbon-cycle response under CO 2 fertilization.

Q. Li et al.: Leaf area index identified as a major source
The reduced magnitudes of β cpool compared with those of β GPP and β NPP in all simulations indicate carbon turnover processes make ecosystems respond to eCO 2 less sensitively due to the slow allocation and carbon turnover processes.A previous study using seven global vegetation models identified carbon residence time as the dominant cause of uncertainty in terrestrial vegetation responses to future climate and atmospheric CO 2 change (Friend et al., 2014).The response of soil carbon storage to eCO 2 also depends on soil carbon residence time (Harrison et al., 1993).In this study and many other models, allocation coefficients were fixed over time (Walker et al., 2014).But allocation patterns to plant organs with different lifespans have been reported to change in response to eCO 2 in experiments, thereby altering carbon residence time in plants and soil (De Kauwe et al., 2014).Therefore, the fixed allocation scheme we adopted in this study might lead to some biases in simulating the response of carbon residence time to eCO 2 .In our study, soil decomposition rate is assumed not to be affected by CO 2 level, as in most other conventional soil carbon models (Friedlingstein et al., 2006;Luo et al., 2016).However, recent synthesis of experimental data suggested that replenishment of new carbon into soil due to eCO 2 increases turnover rate of soil carbon (Van Groenigen et al., 2014;Van Groenigen et al., 2017).Within a certain PFT, the variation of β cpool across different geographical locations is usually smaller than that of β NPP , while the greater variation of β cpool than that of β NPP across different C 3 PFTs in C-only simulation suggests other processes such as different carbon allocation patterns, plant carbon turnover, and the soil carbon dynamics of various PFTs are responsible for the additional divergence.In nutrient-coupled simulations, the variations of β cpool across different C 3 PFTs are only slightly larger than those of β NPP , indicating that nutrients do not bring many differential effects on carbon turnover processes for different PFTs.

Implication for understanding β in other models
Although we analyzed a single land-surface model in detail, the patterns of and mechanisms underlying the variability of β we found may be generally applicable to other models.The basic Farquhar photosynthesis model and two-leaf scaling scheme in the CABLE model are shared by many land-surface models.Some models use variants of the Farquhar photosynthesis model such as the co-limitation approach described by Collatz et al. (1991).Inflection point from Rubisco-to RuBP-limited processes is an important control of the absolute photosynthetic response to eCO 2 (Rogers et al., 2016).However, the relative photosynthetic responses for different ecosystems will converge to a small range because the normalized photosynthetic response to eCO 2 only depends on estimates of intercellular CO 2 concentration (C i ), Michaelis-Menten constants (K c , K o ) and CO 2 compensation point ( * ), and the relative photosynthetic responses are similar for either Rubisco-or RuBP-limited photosynthesis (Luo et al., 1996;Luo and Mooney, 1996).Soil moisture availability is another key constraint on photosynthetic response.Water stress on plants is generally alleviated under eCO 2 due to reduced stomatal conductance (Leuzinger and Körner, 2007;Fatichi et al., 2016).Different models simulate diverse levels of water stress on productivity (De Kauwe et al., 2017).Water stress is simulated in many models to regulate stomatal conductance (Rogers et al., 2016;Wu et al., 2018).For example, the CABLE model represents water stress by an empirical relationship based on soil texture and limits the slope of the coupled relationship between photosynthesis rate and stomatal conductance as Eq.(S11).The influence of water stress is reflected by C i .Synthesis of many empirical study results and our results in this study all show that the ratio of C i to C a is relatively constant, probably due to homeostatic regulations through photosynthetic rate and stomatal conductance (Pearcy and Ehleringer, 1984;Evans and Farquhar, 1991).Wong et al. (1979) showed plant stomata could maintain a constant C i /C a ratio across a wide range of environmental conditions, including the water stress condition.Land-surface models might simulate relatively constant C i /C a ratios under water stress as well since photosynthesis and stomatal conductance are theoretically depicted based on experimental results.Moreover, Luo and Mooney (1996) found that changing the C i /C a ratio from 0.6 to 0.8 caused a variation of less than 0.08 in the sensitivity of leaf photosynthesis to a unit of increase in C a .K c and K o are variable among species, but only slightly affect leaf-level response (Luo and Mooney, 1996).Different leaf temperature will exert a limited influence on the variability of leaf-level β, as we discussed above.Therefore, leaf-level β values for different C 3 PFTs are more likely to converge in other landsurface models.
A recent study used 16 crop models to simulate rice yield at two FACE sites (Hasegawa et al., 2017).These models have diverse representations of primary productivity.Their results showed that the variation of yield response across models was not much associated with model structure or magnitude of primary photosynthetic response to eCO 2 , but was significantly related to the estimations of leaf area.This is consistent with our conclusion and highlights the great need to improve prognostic LAI modeling.Other landsurface modeling groups may benefit from a similar analysis to identify major causes of variability of β across the hierarchical levels from biochemistry to land carbon storage.Candidate causes that can make substantial contributions to the variability include changes in leaf area index, changes in carbon use efficiency and changes in land carbon residence times.If modeling groups can add leaf-level diagnostics in the next inter-model comparison project, it will greatly help disentangle the uncertainty of concentration-carbon feedback.

Conclusions
Exploring the variability of β at different hierarchical levels within and across different C 3 PFTs helps unravel model mechanisms that govern terrestrial ecosystem responses to elevated CO 2 .Our study shows that the sensitivities of biochemistry and leaf-level photosynthesis to eCO 2 are very similar within and across C 3 PFTs in C-only, C-N and C-N-P simulations of CABLE, in accordance with previous theoretical analysis, while β values of GPP, NPP and ecosystem carbon storage diverge primarily because the sensitivities of LAI significantly differ within and across different PFTs in all simulations.After decomposing β into photosynthetic and LAI components, we find LAI contributes more than photosynthesis to the magnitudes and trends of model responses.
Our results indicate that processes related to LAI need to be better constrained with results from experiments and observations in order to better represent the responses of ecosystem carbon cycle processes to changes in CO 2 and climate.
Author contributions.QL and YL designed the study.QL and XL conducted the model simulations.XL, YL, YW and XH helped to solve the problems of coding.QL, XL, YL and PMC analyzed the results together.All the co-authors were involved in writing the paper and contributed to the review process.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.Temporal trends of β cpool from 2011 to 2100 for C 3 PFTs from CABLE-C-only (a), CABLE-CN (b), and CABLE-CNP (c) simulations.β cpool values for different C 3 PFTs all decline with time from 2011 to 2100 under the RCP 8.5 scenario, but the magnitudes of β cpool differ across them in all simulations.In C-N and C-N-P simulations, magnitudes of β cpool are reduced compared with those in C-only simulation for all C 3 PFTs except evergreen broadleaf forest.ENF, evergreen needleleaf forest (light green squares); EBF, evergreen broadleaf forest (red circles); DNF, deciduous needleleaf forest (dark blue triangles); DBF, deciduous broadleaf forest (pink triangles); SHB, shrub (dark green diamonds); C3GRAS, C 3 grass (dark blue stars); TUN, tundra (orange diamonds).

Figure 2 .
Figure 2. Responses of yearly intercellular CO 2 concentration (C i ) to eCO 2 of a single sunlit leaf (a) and shaded leaf (b) for C 3 PFTs from CABLE-C-only simulation.Temporal trends of CO 2 compensation point in the absence of day respiration ( * ) for sunlit leaf (c) and shaded leaf (d) from 2011 to 2100 from CABLE-C-only simulation.The ratio of C i to C a (C i /C a ) is approximately constant with eCO 2 for each PFT and varies little across PFTs.* values vary across different PFTs, but do not change over time for each PFT.Abbreviations and symbols are the same as Fig. 1.

Figure 3 .
Figure 3. Biome-level β values at different levels at the year 2023 from CABLE-C-only (a), CABLE-CN (b), and CABLE-CNP (c) simulations.CV means coefficient of variation of biome-level β across C 3 PFTs.β values at biochemical (L sun and L sha for sunlit and shaded leaves, respectively) and leaf levels (β p sun and β p sha for sunlit and shaded leaves, respectively) are very similar across PFTs, but greatly diverge at canopy level (β GPP ) and ecosystem levels (β NPP and β cpool ) in all simulations.Unlike in C-only simulation, β NPP diverges more than β GPP across different PFTs in nutrient-coupled simulations.Abbreviations and symbols are the same as Fig. 1.

Figure 4 .
Figure 4. Correlations between β GPP and β LAI , β NPP and β LAI , and β cpool and β LAI at the year 2023 across C 3 PFTs from CABLE-Conly (a-c), CABLE-CN (d-f) and CABLE-CNP (g-i) simulations.β GPP , β NPP and β cpool all have significant linear correlations with β LAI in all simulations.Abbreviations and symbols are the same as Fig. 1.

Figure 5 .
Figure 5. Temporal trends of GPP sun (red points) and GPP sha (black points) for C 3 PFTs from 1901 to 2100 from CABLE-C-only simulation.GPP sha is higher than GPP sun for all PFTs.With significant increase in CO 2 concentration from 2011, GPP sha responds more drastically than GPP sun .Abbreviations are the same as Fig.1.
) NPP increase, (2) change in allocation of NPP to leaf, (3) change in specific leaf area (SLA) in response to eCO 2 , and (4) PFT-specific minimum and maximum LAI values prescribed in the model.First, the low responses of LAI to eCO 2 for deciduous needleleaf forest and tundra can be attributed to smaller NPP enhancements in cold areas.The large divergence of the response of LAI within PFTs is mainly due to the large range of NPP increment across different geographical locations.The reduced magnitude of β LAI under nutrient limitations is the direct outcome of reduced β NPP .Accurate estimate of response of GPP and NPP is therefore fundamental to realistic LAI modeling.Second, diverse allocation schemes influence the responses of LAI for different PFTs.And, results from two FACE (Duke Forest and Oak Ridge) experiments indicate that the carbon allocated to leaves is decreased and more

Figure 6 .
Figure 6.Temporal trends of β GPP sun (sensitivity of sunlit leaf GPP; red squares), β GPP sha (sensitivity of shaded leaf GPP; green squares), β S sun (sensitivity of scaling fatcor for sunlit leaf; pink triangles), β S sha (sensitivity of scaling fatcor for shaded leaf; dark blue triangles), β p sun (photosynthetic response for sunlit leaf; purple diamonds) and β p sha (photosynthetic response for shaded leaf; sky blue diamonds) for C 3 PFTs from CABLE-C-only simulation.The sensitivities of GPP sun and GPP sha tend to approach zero through time because the decomposing factors β p sun , β p sha , β S sun and β S sha all decline with time.β S sha determines the magnitudes and trends of β GPP sha for almost all PFTs.Abbreviations are the same as Fig. 1.

Table 1 .
The ratio of intercellular CO 2 concentration (C To reveal which processes cause the large disparity of β across PFTs as shown in Fig.1, we first compared intercellular CO 2 concentration (C i ) and CO 2 compensation point in the absence of day respiration ( * ), which are critical parameters for leaf-level biochemical response.In C-only simulation, the ratio of C i to C a (C i /C a ) is approximately constant with eCO 2 for each PFT (Fig.2a, b).For sunlit leaf, C i /C a values range from 0.64 to 0.70 with CV = 0.03 across different C 3 PFTs (Table1).C i /C a values for shaded leaf are higher than those for sunlit leaf, and the range is 0.68 to 0.76 with CV = 0.03 across different C 3 PFTs (Table1).Evergreen broadleaf forest has the greatest C i /C a value, while deciduous needleleaf forest has the lowest C i /C a value.In C-N simulation, C i /C a values for sunlit leaf are lower than those for the same PFT in C-only simulation, while C i /C a values for shaded leaf change little as compared with those for the same PFT in C-only simulation (Table ciduous broadleaf forest and evergreen broadleaf forest have the greatest β cpool values, while deciduous needleleaf forest and tundra still have the lowest β cpool values in C-N simulation.When both N and P limitations are taken into account as in C-N-P simulation, magnitudes and trends of β cpool are i ) to atmospheric CO 2 concentration (C a ) for different C 3 PFTs, and mean and coefficient of variation (CV) across these PFTs of C i /C a in Conly, C-N, and C-N-P simulations of CABLE under the RCP 8.5 scenario.Values for shaded leaves are in brackets.Abbreviations are the same as Fig. 1.