Leaf area index identified as a major source of variability in modeled CO₂ fertilization

The concentration–carbon feedback (β), also called the CO₂ fertilization effect, is a key unknown in climate–carbon-cycle projections. A better understanding of model mechanisms that govern terrestrial ecosystem responses to elevated CO₂ is urgently needed to enable a more accurate prediction of future terrestrial carbon sink. We conducted C-only, carbon–nitrogen (C–N) and carbon–nitrogen–phosphorus (C–N–P) simulations of the Community Atmosphere Biosphere Land Exchange model (CABLE) from 1901 to 2100 with fixed climate to identify the most critical model process that causes divergence in β. We calculated CO₂ fertilization effects at various hierarchical levels from leaf biochemical reaction and leaf photosynthesis to canopy gross primary production (GPP), net primary production (NPP), and ecosystem carbon storage (cpool) for seven C₃ plant functional types (PFTs) in response to increasing CO₂ under the RCP 8.5 scenario. Our results show that β values at biochemical and leaf photosynthesis levels vary little across the seven PFTs, but greatly diverge at canopy and ecosystem levels in all simulations. The low variation of the leaf-level β is consistent with a theoretical analysis that leaf photosynthetic sensitivity to increasing CO₂ concentration is almost an invariant function. In the CABLE model, the major jump in variation of β values from leaf levels to canopy and ecosystem levels results from divergence in modeled leaf area index (LAI) within and among PFTs. The correlation of β_GPP, β_NPP, or β_cpool each with β_LAI is very high in all simulations. Overall, our results indicate that modeled LAI is a key factor causing the divergence in β in the CABLE model. It is therefore urgent to constrain processes that regulate LAI dynamics in order to better represent the response of ecosystem productivity to increasing CO₂ in Earth system models.

1 Introduction

Terrestrial ecosystems absorb roughly 30 % of anthropogenic CO₂ 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₂ concentration will stimulate plant growth and ecosystem carbon storage, forming a negative feedback to CO₂ concentration (Long et al., 2004; Friedlingstein et al., 2006; Canadell et al., 2007). This concentration–carbon feedback (β), also called the CO₂ 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, 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 models, 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₂ Enrichment (FACE) experimental results. Though satellite products they used may underestimate the effect of CO₂ 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₂ (eCO₂) 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₂ 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₂ (Canadell et al., 2007; Iversen et al., 2012). The leaf-level CO₂ fertilization for C₃ 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₂ sensitivity of leaf-level photosynthesis, which represents kinetics sensitivity of photosynthetic enzymes, varies little among different C₃ species at a given CO₂ concentration (Luo and Mooney, 1996; Luo et al., 1996). However, the CO₂ 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₂. 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₂ 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₂ 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₂ fertilization effects depend on locations, vegetation types and soil nutrient conditions. The strongest absolute CO₂ 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₂ occurs in boreal forests (Joos et al., 2001; Peng et al., 2014). But with gradual eCO₂, 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₂ 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.

2 Materials and methods

2.1 CABLE model description

CABLE (version 2.0) is the Australian community land-surface model (Kowalczyk et al., 2006) and incorporates CASA-CNP to simulate global carbon (C), nitrogen (N) and phosphorus (P) cycles (Wang et al., 2010, 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

$$\begin{array}{}\text{(1)}& \mathrm{LAI}={\mathrm{C}}_{\mathrm{leaf}}\cdot \mathrm{SLA},\end{array}$$

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₂ 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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{NPP}=\mathrm{GPP}\left(\mathrm{LAI},V_{\mathrm{cmax}}\left({\mathrm{N}}_{\mathrm{l}}\right),J_{\max }\left({\mathrm{N}}_{\mathrm{l}}\right)\right)\\ \text{(2)}& {\displaystyle }& {\displaystyle }-{\sum }_i{R}_{\mathrm{mi}}\left({\mathrm{N}}_i\right)-R_{\mathrm{g}}\left({\displaystyle \frac{{\mathrm{N}}_{\mathrm{l}}}{{\mathrm{P}}_{\mathrm{l}}}}\right),\end{array}$$

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⁻²) 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.

2.2 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₂ 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₂ 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₂ 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₂ 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 $\mathrm{1.9}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$ (latitude vs. longitude). N deposition is prescribed from atmospheric transport models (Lamarque et al., 2010, 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).

2.3 Calculation of β values at five hierarchical levels

We aimed to analyze the CO₂ fertilization effects for biochemical reaction (L), leaf photosynthesis rate (p), leaf-to-canopy 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₂ 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-to-canopy 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₂ concentration (C_a) as β_V:

$$\begin{array}{}\text{(3)}& {\mathit{\beta }}_V={\displaystyle \frac{\mathrm{1}}{V}}\cdot {\displaystyle \frac{\mathrm{d}V}{\mathrm{d}{C}_{\mathrm{a}}}},\end{array}$$

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⁻¹. 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₂ concentration (C_i) and has been suggested to be an invariant function for C₃ plants grown in diverse environments. The rate of photosynthesis is typically RuBP-regeneration-limited under a high CO₂ 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₂. For sunlit leaf and shaded leaf, formulations of L under RuBP-regeneration limitation are defined as

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }{L}_{\mathrm{sun}}={\displaystyle \frac{\mathrm{3}\cdot {\mathrm{\Gamma }}_{\ast \mathrm{sun}}}{\left(C_{i\mathrm{sun}}+\mathrm{2}\cdot {\mathrm{\Gamma }}_{\ast \mathrm{sun}}\right)\left(C_{i\mathrm{sun}}-{\mathrm{\Gamma }}_{\ast \mathrm{sun}}\right)}},\text{(5)}& {\displaystyle }& {\displaystyle }{L}_{\mathrm{sha}}={\displaystyle \frac{\mathrm{3}\cdot {\mathrm{\Gamma }}_{\ast \mathrm{sha}}}{\left(C_{i\mathrm{sha}}+\mathrm{2}\cdot {\mathrm{\Gamma }}_{\ast \mathrm{sha}}\right)\left(C_{i\mathrm{sha}}-{\mathrm{\Gamma }}_{\ast \mathrm{sha}}\right)}}.\end{array}$$

In this study, Γ_∗sun and Γ_∗sha are yearly average CO₂ 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₂ 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 $\frac{\mathrm{d}{C}_i}{\mathrm{d}{C}_{\mathrm{a}}}$. For sunlit leaf and shaded leaf, the formulations are

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sun}}}=L_{\mathrm{sun}}\cdot {\displaystyle \frac{\mathrm{d}{C}_{i\mathrm{sun}}}{\mathrm{d}{C}_{\mathrm{a}}}},\text{(7)}& {\displaystyle }& {\displaystyle }{\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sha}}}=L_{\mathrm{sha}}\cdot {\displaystyle \frac{\mathrm{d}{C}_{i\mathrm{sha}}}{\mathrm{d}{C}_{\mathrm{a}}}}.\end{array}$$

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

$$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle }{S}_{\mathrm{sun}}={\displaystyle \frac{\mathrm{1}-\exp [-\mathrm{LAI}\left(k_{\mathrm{n}}+k_{\mathrm{b}}\right)]}{k_{\mathrm{n}}+k_{\mathrm{b}}}},\text{(9)}& {\displaystyle }& {\displaystyle }{S}_{\mathrm{sha}}={\displaystyle \frac{\mathrm{1}-\exp \left(-k_{\mathrm{n}}\mathrm{LAI}\right)}{k_{\mathrm{n}}}}-{\displaystyle \frac{\mathrm{1}-\exp \left[-\mathrm{LAI}\left(k_{\mathrm{n}}+k_{\mathrm{b}}\right)\right]}{k_{\mathrm{n}}+k_{\mathrm{b}}}},\end{array}$$

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.

Figure 1 Temporal trends of β_cpool from 2011 to 2100 for C₃ PFTs from CABLE-C-only (a), CABLE-CN (b), and CABLE-CNP (c) simulations. β_cpool values for different C₃ 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₃ 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₃ grass (dark blue stars); TUN, tundra (orange diamonds).

Big-leaf ${\mathit{\beta }}_{{\mathrm{GPP}}_{\mathrm{sun}}}$ (or ${\mathit{\beta }}_{{\mathrm{GPP}}_{\mathrm{sha}}}$) can be decomposed as the sum of normalized sensitivity of photosynthesis rate, ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sun}}}$ (or ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sha}}}$), and leaf-to-canopy scaling factor, ${\mathit{\beta }}_{{\mathrm{S}}_{\mathrm{sun}}}$ (or ${\mathit{\beta }}_{{\mathrm{S}}_{\mathrm{sha}}}$), as shown in Eqs. (10) and (11). Detailed mathematical derivations are in Supplement S2.

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{\mathit{\beta }}_{{\mathrm{GPP}}_{\mathrm{sun}}}={\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sun}}}+{\mathit{\beta }}_{{\mathrm{S}}_{\mathrm{sun}}}\text{(11)}& {\displaystyle }& {\displaystyle }{\mathit{\beta }}_{{\mathrm{GPP}}_{\mathrm{sha}}}={\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sha}}}+{\mathit{\beta }}_{{\mathrm{S}}_{\mathrm{sha}}}\end{array}$$

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₃ 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, ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sun}}}$, ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sha}}}$, β_GPP, β_NPP and β_cpool at the year 2023 (relative to 2022) for different C₃ PFTs were calculated and compared. Coefficients of variation (CVs) of β values were calculated across various C₃ 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₃ 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.

3 Results

3.1 Temporal trends of β at ecosystem level for different PFTs

In C-only simulation, β_cpool values for different C₃ PFTs all decline with time from 2011 to 2100 under the RCP 8.5 scenario (Fig. 1a). However, the magnitudes of β_cpool differ among different PFTs, with the highest values occurring in deciduous broadleaf forest from 2011 to 2075 and in shrub after 2075, and the lowest values occurring in deciduous needleleaf forest and tundra. β_cpool values for deciduous needleleaf forest and tundra nearly overlap over time. As compared with C-only simulation, values of β_cpool are reduced when N limitation is included as in C–N simulation for all C₃ PFTs except evergreen broadleaf forest (Fig. 1b). Deciduous 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 similar to those in C–N simulation (Fig. 1c) as P limitation is quite weak under present conditions in the current version of CABLE (Zhang et al., 2011).

Table 1 The ratio of intercellular CO₂ concentration (C_i) to atmospheric CO₂ concentration (C_a) for different C₃ PFTs, and mean and coefficient of variation (CV) across these PFTs of C_i∕C_a in C-only, 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.

3.2 Variations of intercellular CO₂ concentration and CO₂ compensation point

To reveal which processes cause the large disparity of β across PFTs as shown in Fig. 1, we first compared intercellular CO₂ concentration (C_i) and CO₂ 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₂ 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₃ PFTs (Table 1). 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₃ PFTs (Table 1). 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 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).

Figure 2 Responses of yearly intercellular CO₂ concentration (C_i) to eCO₂ of a single sunlit leaf (a) and shaded leaf (b) for C₃ PFTs from CABLE-C-only simulation. Temporal trends of CO₂ 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₂ 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.

In all of the simulations, values of CO₂ 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₃ PFTs because of different leaf temperature which Γ_∗ values depend on.

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₃ PFTs. β values at biochemical (L_sun and L_sha for sunlit and shaded leaves, respectively) and leaf levels (${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sun}}}$ and ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{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.

3.3 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, ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sun}}}$, ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sha}}}$, β_GPP, β_NPP and β_cpool at the year 2023 for different C₃ 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₃ 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⁻¹. Variations of L_sun and L_sha among PFTs are small (CV = 0.15 and 0.13). At leaf photosynthesis level, ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sun}}}$ and ${\mathit{\beta }}_{p_{\mathrm{sha}}}$ for the seven PFTs vary from 0.00041 to 0.00072 ppm⁻¹, 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₃ 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⁻¹). Deciduous broadleaf forest has the greatest β_cpool value (around 0.0018 ppm⁻¹) 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.

Table 2 Coefficients of variation of L, β_p, β_GPP, β_NPP and β_cpool across different geographical locations within each C₃ PFT at the year 2023 in the CABLE-C-only simulation. The two numbers in the same unit are for sunlit leaves and shaded leaves, respectively. Values for shaded leaves are in brackets. Abbreviations are the same as Fig. 1.

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 nutrient-limited β_GPP across different PFTs, which is similar to C-only simulation. However, unlike in C-only simulation, β_NPP values in nutrient-coupled simulations are reduced for most C₃ 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.

Figure 4 Correlations between β_GPP and β_LAI, β_NPP and β_LAI, and β_cpool and β_LAI at the year 2023 across C₃ PFTs from CABLE-C-only (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.

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₃ 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₃ 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₃ 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₂.

Figure 5 Temporal trends of GPP_sun (red points) and GPP_sha (black points) for C₃ 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₂ concentration from 2011, GPP_sha responds more drastically than GPP_sun. Abbreviations are the same as Fig. 1.

3.4 β 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₃ PFTs from 1901 to 2100 in Fig. 5. From the beginning of the simulation, GPP_sha is higher than GPP_sun for all C₃ PFTs. With significant increases in CO₂ concentration from 2011, GPP_sha responds more drastically than GPP_sun. Shaded leaf GPP of deciduous broadleaf forest and shrub responds to eCO₂ 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 LAI-dependent 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 ${\mathit{\beta }}_{{\mathrm{GPP}}_{\mathrm{sun}}}\left({\mathit{\beta }}_{{\mathrm{GPP}}_{\mathrm{sha}}}\right)$ and decomposing factors ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sun}}}\left({\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sha}}}\right)$ and ${\mathit{\beta }}_{{\mathrm{S}}_{\mathrm{sun}}}\left({\mathit{\beta }}_{{\mathrm{S}}_{\mathrm{sha}}}\right)$ 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 ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sun}}}$, ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sha}}}$, ${\mathit{\beta }}_{{\mathrm{S}}_{\mathrm{sun}}}$ and ${\mathit{\beta }}_{{\mathrm{S}}_{\mathrm{sha}}}$ all decline with time (Fig. 6). ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sun}}}$ and ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sha}}}$ overlap through time for each PFT. Magnitudes of ${\mathit{\beta }}_{{\mathrm{GPP}}_{\mathrm{sha}}}$ are higher than those of ${\mathit{\beta }}_{{\mathrm{GPP}}_{\mathrm{sun}}}$ for all C₃ PFTs because the values of ${\mathit{\beta }}_{{\mathrm{S}}_{\mathrm{sha}}}$ are higher than those of ${\mathit{\beta }}_{{\mathrm{S}}_{\mathrm{sun}}}$. For deciduous needleleaf forest and tundra, both ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sun}}}$ (${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{sha}}}$) and ${\mathit{\beta }}_{{\mathrm{S}}_{\mathrm{sun}}}\left({\mathit{\beta }}_{{\mathrm{S}}_{\mathrm{sha}}}\right)$ contribute to the magnitudes and trends of ${\mathit{\beta }}_{{\mathrm{GPP}}_{\mathrm{sun}}}\left({\mathit{\beta }}_{{\mathrm{GPP}}_{\mathrm{sha}}}\right)$. For evergreen needleleaf forest, deciduous broadleaf forest, shrub and C₃ grass, ${\mathit{\beta }}_{{\mathrm{S}}_{\mathrm{sun}}}\left({\mathit{\beta }}_{{\mathrm{S}}_{\mathrm{sha}}}\right)$ dominates the magnitude and change in ${\mathit{\beta }}_{{\mathrm{GPP}}_{\mathrm{sun}}}\left({\mathit{\beta }}_{{\mathrm{GPP}}_{\mathrm{sha}}}\right)$. For evergreen broadleaf forest, ${\mathit{\beta }}_{{\mathrm{S}}_{\mathrm{sha}}}$ predominates the magnitude and change in ${\mathit{\beta }}_{{\mathrm{GPP}}_{\mathrm{sha}}}$ before 2035.

4 Discussion

4.1 Variations of biochemical and photosynthetic responses to eCO₂

The direct CO₂ fertilization effect occurs at leaf level and is determined by kinetic sensitivity of Rubisco enzymes to internal leaf CO₂ concentration. In fact, the normalized short-term sensitivity of leaf-level photosynthesis to CO₂ 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₂ 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₂ 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 $\pm \mathrm{5}{}^{\circ }$ of leaf temperature changes caused approximately ±7 ppm changes in Γ_∗, leading to variation of 0.12 to leaf-level β (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₂ 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 big-leaf, 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₃ 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 leaf-level β 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.

4.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₂ diverge primarily because the responses of LAI diverge within and among PFTs in all simulations. Besides, GPP of shaded leaves responds to eCO₂ more strongly than GPP of sunlit leaves for all C₃ 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₂. For shaded leaves, the sensitivity of the LAI-dependent scaling factor contributes more to the magnitude and trend of ${\mathit{\beta }}_{{\mathrm{GPP}}_{\mathrm{sha}}}$ than that of photosynthesis rate. 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₂ during the historical period (Zhu et al., 2016). Our modeling study also shows that LAI responds positively to eCO₂ for all C₃ 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₃ 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 EucFACE probably leads to an insignificant response of productivity at this site, even though leaf photosynthesis rate significantly increases under eCO₂ (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₃ 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₂, as this study has shown that LAI is an important factor in driving ecosystem response to CO₂ fertilization. And it is also likely the overestimation of the response of LAI to eCO₂ is responsible for the overestimation of CO₂ fertilization in ESMs reported by previous studies (Smith et al., 2015; Mystakidis et al., 2017).

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

The overall response of LAI to eCO₂ depends on several processes in this study: (1) NPP increase, (2) change in allocation of NPP to leaf, (3) change in specific leaf area (SLA) in response to eCO₂, and (4) PFT-specific minimum and maximum LAI values prescribed in the model. First, the low responses of LAI to eCO₂ 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 carbon is allocated to woods or roots at higher CO₂ 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₂ 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₂. Fourth, in our results, LAI values for most C₃ PFTs are below the maximum LAI limits with eCO₂ 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₂ 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₂ (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₃ PFTs in C-only simulation suggest carbon use efficiency (CUE) does not change with eCO₂, as autotrophic respiration is calculated from GPP and plant carbon. In C–N and C–N–P simulations, magnitudes of β_NPP for all C₃ PFTs except evergreen broadleaf forest all decline compared with those of β_GPP, indicating CUE also declines with eCO₂ under nutrient limitations. However, FACE experimental results indicate that CUE values under eCO₂ 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₃ PFTs in nutrient-coupled simulations as in C-only simulation, confirming the dominant role of LAI in regulating carbon-cycle response under CO₂ fertilization.

The reduced magnitudes of β_cpool compared with those of β_GPP and β_NPP in all simulations indicate carbon turnover processes make ecosystems respond to eCO₂ 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₂ change (Friend et al., 2014). The response of soil carbon storage to eCO₂ 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₂ 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₂. In our study, soil decomposition rate is assumed not to be affected by CO₂ 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₂ 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₃ 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₃ 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.

4.3 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₂ (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₂ only depends on estimates of intercellular CO₂ concentration (C_i), Michaelis–Menten constants (K_c, K_o) and CO₂ 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₂ 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₃ PFTs are more likely to converge in other land-surface 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₂, 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 land-surface 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.

5 Conclusions

Exploring the variability of β at different hierarchical levels within and across different C₃ PFTs helps unravel model mechanisms that govern terrestrial ecosystem responses to elevated CO₂. Our study shows that the sensitivities of biochemistry and leaf-level photosynthesis to eCO₂ are very similar within and across C₃ 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₂ and climate.