Abstract

Carbon uptake by land and ocean as a biogeochemical response to increasing atmospheric CO2 concentration is called concentration–carbon feedback and is one of the carbon cycle feedbacks of the global climate. This feedback can have a major impact on climate projections with an uncertain magnitude. This paper focuses on the concentration–carbon feedback in terrestrial ecosystems, analyzing the mechanisms and strength of the feedback reproduced by Earth system models (ESMs) participating in phase 5 of the Coupled Model Intercomparison Project. It is confirmed that multiple ESMs driven by a common scenario show a large spread of concentration–carbon feedback strength among models. Examining the behavior of the carbon fluxes and pools of the models showed that the sensitivity of plant productivity to elevated CO2 is likely the key to reduce the spread, although increasing CO2 stimulates other carbon cycle processes. Simulations with a single ESM driven by different CO2 pathways demonstrated that carbon accumulation increases in scenarios with slower CO2 increase rates. Using both numerical and analytical approaches, the study showed that the difference among CO2 scenarios is a time lag of terrestrial carbon pools in response to atmospheric CO2 increase—a high rate of CO2 increase results in smaller carbon accumulations than that in an equilibrium state of a given CO2 concentration. These results demonstrate that the current quantities for concentration–carbon feedback are incapable of capturing the feedback dependency on the carbon storage state and suggest that the concentration feedback can be larger for future scenarios where the CO2 growth rate is reduced.

1. Background

Global carbon cycling in atmospheric, ocean, and terrestrial ecosystems plays a key role in long-term climate projections because carbon cycle feedbacks can have a great impact on the global climate through changing atmosphere–land and atmosphere–ocean carbon budgets, especially at centennial or longer time scales. Earth system models (ESMs) have been developed to explicitly incorporate climate–carbon cycle interactions by implementing biogeochemical cycle processes in atmosphere–ocean general circulation models (GCMs). These fully coupled, climate–carbon cycle models make it possible to produce climate change projections using internally estimated atmospheric CO2 concentrations (e.g., Friedlingstein et al. 2006) or a diagnosis of anthropogenic CO2 emissions compatible with a prescribed CO2 concentration (e.g., Jones et al. 2013) with the climate–carbon cycle interactions.

There are two primary types of carbon cycle feedbacks: concentration–carbon and climate–carbon (Friedlingstein et al. 2003). The former generates a feedback loop between atmospheric CO2 and the carbon cycle of the land and ocean. This promotes CO2 uptake by terrestrial ecosystems by stimulating photosynthesis and plant growth and producing a large atmosphere–ocean CO2 difference, which intensifies CO2 dissolution and diffusion in the ocean. The climate–carbon feedback also generates another feedback loop between the climate and carbon cycle: global warming can induce greater ecosystem respiration and a resultant increase of CO2 in the atmosphere, further accelerating global warming. The two feedbacks on the global carbon cycle can be briefly summarized by the following equation, using the definition of Friedlingstein et al. (2003):

 
formula

where ΔC represents changes in the carbon amount of the land L and ocean O from the preindustrial state and is expressed as the sum of the two types of carbon cycle feedback: carbon amount change owing to concentration–carbon feedback (first term on the right) and climate–carbon feedback (second term). In this formula, the carbon storage change caused by concentration–carbon feedback is treated as proportional to the CO2 increase ΔCO2 (ppmv) with feedback parameter β (PgC ppmv−1), and the change caused by climate–carbon feedback is represented as a function of global temperature change ΔT (kelvin) with feedback parameter γ (PgC K−1).

To evaluate these carbon cycle feedbacks and their parameters, Friedlingstein et al. (2006) and Arora et al. (2013) conducted multimodel intercomparisons and analyzed carbon cycle feedbacks by compiling results of simulations run systematically under the same experimental protocols. To evaluate the magnitude of the concentration–carbon feedback [parameter β in Eq. (1)], they used results from an experiment called an “uncoupled simulation” or biogeochemically coupled (BGC) simulation. The atmospheric CO2 is allowed to vary but does not produce any radiative imbalance in the atmosphere or resultant global warming (although there can be a slight indirect warming or cooling effect, owing to changes in land surface conditions). From their multimodel intercomparisons, because elevated CO2 levels stimulate carbon accumulation in the land and ocean, this feedback is likely to reduce atmospheric CO2, forming a negative feedback loop (the sign of β should be positive in this case).

In contrast, the strength of the climate–carbon feedback [parameter γ in Eq. (1)] is evaluated via a radiatively coupled simulation (Arora et al. 2013), in which atmospheric CO2 is allowed to increase and induce global warming, but the elevated CO2 does not directly stimulate terrestrial or ocean carbon cycles. Thus, in a radiatively coupled simulation experiment, global warming is the primary driver that alters carbon cycles. Instead of conducting such an experiment, Friedlingstein et al. (2006) evaluated climate–carbon feedback by examining the difference between the results of a fully coupled experiment and BGC simulation. Although there are methodological differences between Friedlingstein et al. (2006) and Arora et al. (2013), both studies show a common trend that global warming is likely to reduce the carbon amount in the land and ocean, which results in more CO2 in the atmosphere. This means that the climate–carbon cycle forms a positive feedback loop (i.e., a negative sign of γ).

These sensitivity simulations have revealed the relative importance of carbon cycle feedbacks within the earth system. In an intermodel comparison, Arora et al. (2013) showed that the contribution of concentration–carbon feedback to the total global carbon budget is about 4.5 times greater than the carbon–climate feedback in an idealized CO2 increase scenario. Gregory et al. (2009) equivalently treated the strength of carbon cycle feedbacks with other physical climate feedbacks and compared them, based on the results of Friedlingstein et al. (2006). Gregory et al. (2009) found that on average, the climate–carbon feedback is similar in magnitude to the cloud feedback. However, the magnitude of the concentration–carbon feedback was around 4 times that of the climate–carbon feedback and comparable in magnitude to the total noncarbon feedback.

In addition to having large magnitude, carbon cycle feedback can have large uncertainties, especially concentration–carbon feedback in terrestrial ecosystems. For example, Gregory et al. (2009) showed that concentration–carbon feedback can be a greater source of projection uncertainty than climate–carbon feedback. Arora et al. (2013) reported that the ocean carbon cycle showed less variation in both concentration– and climate–carbon feedbacks, whereas in terrestrial ecosystems carbon cycle feedbacks (especially concentration–carbon feedback), it differed greatly between models. In addition to this “intermodel” spread of the concentration–carbon feedback in multimodel simulation, the strength of the feedback probably depends on the rate of CO2 increase, that is, a model may produce different magnitudes of concentration–carbon feedback when various CO2 scenarios are applied (Boer and Arora 2009; Gregory et al. 2009). We call this an “interscenario” spread. For example, Gregory et al. (2009) conducted BGC experiments for scenarios with different rates of CO2 increase, showing that concentration–carbon feedback in ESMs was more scenario dependent than climate–carbon feedback. Similar results were also shown by Boer and Arora (2009), who used a new approach to evaluate the feedback strength. Moreover, Hajima et al. (2012) suggested that this dependence of concentration–carbon feedback on CO2 growth rate could greatly vary the airborne fraction among simulations with representative concentration pathways (RCPs; Moss et al. 2010). Thus, the combination of intermodel spread derived from different modeling/parameters and interscenario spread arising from the rate of CO2 increase assumed in the scenarios can produce large uncertainties in the anthropogenic CO2 emissions that are compatible with a target pathway, being affected by other factors such as land use change scenarios (Jones et al. 2013).

In this study, we investigated mechanisms underlying terrestrial concentration–carbon feedback to reveal uncertain processes within that feedback. To this end, we posed two questions. First, what is the cause of variation in terrestrial concentration–carbon feedback among simulations by multiple ESMs? Arora et al. (2013) showed that this feedback in the models from phase 5 of the Coupled Model Intercomparison Project (CMIP5; Taylor et al. 2012) had a large spread, but the cause of that spread remains unknown. Although Anav et al. (2013) and Todd-Brown et al. (2013) showed in their model intercomparison studies that plant production and the amount of soil carbon might be key factors accounting for variation among CMIP5 ESMs. However, because they analyzed historical simulations, it is difficult to separate the effects of concentration–carbon feedback from those of climate–carbon feedback. To reveal the cause of intermodel variation in concentration–carbon feedback on the terrestrial carbon cycle, we analyzed the feedback strength and related variables of CMIP5 ESMs in simulations with a common CO2 scenario.

The second question is what is the cause of variations in terrestrial concentration–carbon feedback among scenarios with different rates of CO2 increase? It has been clearly demonstrated by ESMs that terrestrial carbon change caused by concentration–carbon feedback differs greatly depending on the CO2 increase rate of the scenario. Consequently, the feedback strength should vary among the scenarios. This suggests that current feedback quantities are incapable of capturing the dependency of the feedback on the system state (Boer and Arora 2009; Gregory et al. 2009). Therefore, to more adequately evaluate the feedback strength and make more precise climate projections, it is necessary to make clear the detailed mechanisms underlying this interscenario difference. To advance the findings by Boer and Arora (2009) and Gregory et al. (2009), we used a common ESM and applied the model to several types of scenarios with various CO2 pathways.

We present details of the multiple ESM simulations with a common scenario in section 2a. Multiscenario simulations using a single ESM are summarized in section 2b. Other methods used to analyze ESMs are summarized in section 3. Simulation results and their analyses are described in section 4a for multimodel simulations, section 4b for multiscenario simulations, and their simultaneous comparisons are in section 4c. In section 5, we discuss nonintuitive feedback behaviors by examining analytical solutions obtained by a simple carbon cycle model. Subsequently, we offer suggestions for reducing uncertainty. Finally, we summarize our findings in section 6.

2. Methods: Simulations by ESMs

a. Biogeochemically coupled simulations by CMIP5 ESMs with a common scenario

To investigate the intermodel spread of the terrestrial concentration–carbon feedback, we used outputs from eight CMIP5 ESMs in an experiment specifically designed for analyzing the sensitivity of the concentration–carbon feedback (Table 1, see table for full model expansions). This experiment was called “esmFixClim1” in the CMIP5 experimental protocols, and the model settings were the same as those in the BGC experiment by Arora et al. (2013) described in section 1. In this experiment, atmospheric CO2 increases at 1% yr−1 (ppy) for 140 yr, from preindustrial levels to a final concentration of about 1140 ppmv, but the increase does not influence atmospheric radiative transfer and thus does not cause global warming. By downloading data from the CMIP5 Earth System Grid (http://cmip-pcmdi.llnl.gov/cmip5/), we obtained outputs from experiments performed with the following eight ESMs with BGC experimental settings: BCC-CSM1.1, CanESM2, IPSL-CM5A-LR, MIROC-ESM, HadGEM2-ES, MPI-ESM-LR, NorESM1-ME, and CESM1(BGC). The basic simulation results with biogeochemically coupled setting using these models, especially for terrestrial carbon change (Table 2; Fig. 1a), have been published by Arora et al. (2013). However, our analysis does not include results from the University of Victoria Earth System Climate Model, version 2.9 (UVic ESCM2.9).

Table 1.

Biogeochemically coupled simulation settings and participating ESMs and all model expansions. For final CO2 level, the CO2 concentration in preindustrial state (284.6 ppmv) is represented by CO2PI.

Biogeochemically coupled simulation settings and participating ESMs and all model expansions. For final CO2 level, the CO2 concentration in preindustrial state (284.6 ppmv) is represented by CO2PI.
Biogeochemically coupled simulation settings and participating ESMs and all model expansions. For final CO2 level, the CO2 concentration in preindustrial state (284.6 ppmv) is represented by CO2PI.
Table 2.

Increments of carbon pools [CL; CV; and CS (litter + soil + other carbon pools)] and fluxes (GPP, NPP, and HR) in ESMs during a 140-yr simulation with 1-ppy CO2 increment scenario. Each value is calculated as the difference between means of the last and first decades.

Increments of carbon pools [CL; CV; and CS (litter + soil + other carbon pools)] and fluxes (GPP, NPP, and HR) in ESMs during a 140-yr simulation with 1-ppy CO2 increment scenario. Each value is calculated as the difference between means of the last and first decades.
Increments of carbon pools [CL; CV; and CS (litter + soil + other carbon pools)] and fluxes (GPP, NPP, and HR) in ESMs during a 140-yr simulation with 1-ppy CO2 increment scenario. Each value is calculated as the difference between means of the last and first decades.
Fig. 1.

Changes in (a) CL, (b) NPP, and (c) HR in simulations by eight ESMs: BCC-CSM1.1 (black), CanESM2 (green), HadGEM2-ES (purple), IPSL-CM5A-LR (red), MIROC-ESM (blue), MPI-ESM-LR (gray), NorESM1-ME (brown), and CESM1(BGC) (pink). Biogeochemically coupled simulations were conducted using a scenario with a 1-ppy CO2 increment.

Fig. 1.

Changes in (a) CL, (b) NPP, and (c) HR in simulations by eight ESMs: BCC-CSM1.1 (black), CanESM2 (green), HadGEM2-ES (purple), IPSL-CM5A-LR (red), MIROC-ESM (blue), MPI-ESM-LR (gray), NorESM1-ME (brown), and CESM1(BGC) (pink). Biogeochemically coupled simulations were conducted using a scenario with a 1-ppy CO2 increment.

b. Biogeochemically coupled simulations by MIROC-ESM with different rate of CO2 increase scenarios

To investigate the scenario dependence of the concentration–carbon feedback in terrestrial ecosystems, we ran BGC simulations with MIROC-ESM (Watanabe et al. 2011) for two groups of CO2 scenarios (Table 1). The first scenario group is named TRN and was designed to investigate the transient state of the terrestrial carbon cycle. It consists of three CO2 concentration scenarios: CO2 concentration monotonically increased from preindustrial levels (CO2PI), with rates of 1) 0.5 ppy for 280 yr, 2) 1.0 ppy for 140 yr, and 3) 2.0 ppy for 70 yr. At the end of these three simulations, the CO2 concentration reaches 1140 ppmv (4 × CO2PI).

The second scenario group, called STB, was designed to examine the carbon storage behavior after stabilization of CO2 concentration and is composed of four scenarios: 1) abrupt CO2 doubling from the preindustrial state and stabilized, 2) CO2 increase at 2.0 ppy from the preindustrial state and stabilized at 2 × CO2PI level, 3) CO2 increase at 1.0 ppy and stabilized at 2 × CO2PI level, and 4) CO2 increase at 0.5 ppy until it reaches 2 × CO2PI level. All four scenarios in the STB group have identical simulation years (140 yr) and the CO2 concentration ultimately reached the 2 × CO2PI level. For both scenario groups, other external anthropogenic forcings, such as non-CO2 greenhouse gas (GHG) concentrations, aerosol emissions, and land-cover changes, were fixed at preindustrial levels.

3. Analysis methods

a. Basic analysis: Relationships among changes in land carbon, net primary production, and heterotrophic respiration

In the BGC simulation, the strength of the concentration–carbon feedback in terrestrial ecosystems was quantified in terms of the change in total terrestrial carbon ΔCL. Because in the BGC experiment there was neither land use change nor resulting disturbance of land carbon, ΔCL must equal the time integration of net ecosystem production NEP(t):

 
formula

where NPP(t) and HR(t), respectively, represent net primary production and heterotrophic respiration as a function of time t, and T represents an integration period. The term ΔCL would also include carbon changes derived from biomass burning, so this equation is not strictly correct. However, because our objective here is to examine carbon dynamics on global and centennial time scales, we can regard ΔCL as identical to the time integral of NEP. The values NPP(t) and HR(t) can be rewritten as

 
formula

and

 
formula

where the subscript 0 denotes the variable in the steady state, which usually appears in the final phase of model spinup, and ΔX represents the changes in variable X.

If the spinup is long enough for the equilibrium state to be reached and the interannual variability of NPP and HR is sufficiently small, then the relationship should hold and thus ΔCL in Eq. (2) should be

 
formula

This equation with some simple assumptions indicates that strength of the concentration–carbon feedback represented in terms of ΔCL can be evaluated from the balance between ΔNPP and ΔHR and does not depend on initial values of NPP and HR.

b. An identity to decompose land carbon change into factors

To determine the origin of intermodel and interscenario spreads of concentration–carbon feedback, we adopted a method similar to the “growth analysis” technique commonly used in ecological studies to analyze individual plant growth, in which the growth rate is expressed as the product of several factors (e.g., Garnier 1992). A similar method has been used in the development of future emissions scenarios (Nakicenovic et al. 2000) by representing GHG emissions as the product of various factors such as economic development and energy intensity.

In our analyses, we first applied the growth analysis technique to ΔCL/ΔCO2, which is identical to β in Eq. (1), to decompose it into component factors and obtain

 
formula

where GPP and HR are carbon fluxes of gross primary production and heterotrophic respiration, respectively. The variables CV and CS represent carbon stored in vegetation and soil, respectively. The first factor on the right-hand side of this equation indicates the GPP response to increased atmospheric CO2. The second factor represents the part of ΔGPP used for plant dry-matter production, not for plant respiration. The third factor represents the increase in vegetation carbon per unit assimilation, likely related to litter fall or the mortality rate of plants. The fourth factor shows the ratio of ΔCS and ΔCV, the carbon allocation in terrestrial ecosystems. The ΔHR/ΔCS as the fifth factor is analogous to the decomposition rate in the form of a first-order decay function. Finally, ΔCL/ΔHR can be regarded as the bulk turnover rate of terrestrial carbon in a transitional state, because the numerator is the total carbon change and the denominator is the carbon efflux from the terrestrial pool.

In Eq. (6), the response of plant production to elevated CO2 concentration is summarized by the first factor (ΔGPP/ΔCO2). However, in this form, it is unclear whether the GPP increment is brought about by sensitive modeling of leaf-scale photosynthesis or by greater allocation of carbon into the leaf component, which usually expands the total leaf area. Indeed, global gross primary productivity is calculated by

 
formula

where P is the gross photosynthetic rate at leaf scale depending on CO2, L is the leaf area density, and A is the unit land area. In this form, it is clear that GPP is dependent on the magnitude of P and on the leaf area index (LAI) that appears at the upper limit of the finite integration. Although it is desirable to separately compare ESMs of these processes, limited data availability prevents detailed analysis of GPP. Alternatively, by introducing a new factor ΔLAI (the change in a global annual mean of leaf area index), we can relate the leaf area index to the factor ΔGPP/ΔCO2 using

 
formula

By combining Eqs. (6) and (7), we obtain

 
formula

By summarizing ΔGPP/ΔLAI as Δgpp (and thus ΔLAI = ΔGPP/Δgpp), Eq. (8a) can be expressed as

 
formula

which we hereafter use in our analysis.

We investigated the origin of intermodel and interscenario spreads of the concentration–carbon feedback by applying this equation to the ESM outputs. Note that the factors used in the growth analysis are not necessarily independent. For example, in Eq. (8b), the factors ΔHR/ΔCS and ΔCSCV are not independent because a higher decomposition rate gives larger values of HR, which reduces the proportion of soil carbon in land ecosystems. By contrast, factors ΔHR/ΔCS and ΔNPP/ΔGPP in Eq. (8b) may be independent because the former relates to soil decomposition and the latter to autotrophic respiration, which in most cases are modeled independently of each other.

c. Simplified terrestrial carbon model

We used a global terrestrial carbon cycle model with reduced complexity to mimic and analyze behaviors of the ESM. This simple type of carbon cycle model has long been useful for interpreting carbon dynamics reproduced in more complex models (e.g., Friedlingstein et al. 1995; Thompson et al. 1996; Todd-Brown et al. 2013) or analyzing the global carbon budget (Taylor and Lloyd 1992). The simplified terrestrial carbon model used here (hereafter STC) consists of only two global carbon pools (CV and CS), the sizes of which are determined by the balance between influx and efflux. In the case of vegetation carbon, because the influx is GPP and efflux is the sum of autotrophic respiration (AR) and litter fall rate (LF), the carbon balance in this pool can be expressed as

 
formula

Because the focus here is the response of the terrestrial carbon cycle to elevated CO2 concentration, GPP is modeled simply with one explanatory variable ΔCO2, and the form of the function is derived from the Michaelis–Menten equation:

 
formula

where GPP0 denotes the initial value of GPP, for which we used the GPP value simulated by the ESM whose behavior is mimicked by the STC model. The terms kGPP1 and kGPP2 are analogous to the kinetic parameters of the Michaelis–Menten equation. The AR and LF are modeled by a first-order decay function that describes the respiration response to changes in the carbon pool:

 
formula

and

 
formula

where kAR and kLF are constants. An equation of this form is commonly used in soil decomposition submodels and for calculation of maintenance respiration in plants.

Similarly for soil carbon, the carbon balance equation should be

 
formula

and

 
formula

where kHR is a constant for the soil decomposition process.

The initial values of vegetation and soil carbon are needed to calculate the temporal evolution of CV and CS. Thus, initial conditions are calculated by the following equations, under the assumption that the left-hand sides of Eqs. (9) and (13) are zero in a steady-state balance:

 
formula

and

 
formula

We applied this simplified global carbon model [Eqs. (9)(16)] to the simulation results of the ESMs.

This kind of simplified model can be used to produce a number of simulation ensembles with perturbing parameters (e.g., Meinshausen et al. 2009) or to conduct simulations over periods longer than a millennium (e.g., Raupach et al. 2011). In this study, as in Todd-Brown et al. (2013), we used the simplified carbon cycle model to mimic ESM behavior, and we analyzed the terrestrial concentration–carbon feedback by solving Eqs. (9)(16) numerically or analytically. BGC experiments by ESMs may project slight warming on regional or global scales. Thus, some climate–carbon feedback might even occur in the BGC simulation and modify the strength of the concentration–carbon feedback by perturbing the terrestrial water cycle (e.g., by modifying water use efficiency). However, we expected the effect of this process to be limited, so we omitted it from the model described by Eqs. (9)(16).

4. Results

a. Concentration–carbon feedbacks among CMIP5 ESMs

1) Land carbon, net primary production, and heterotrophic respiration among ESMs

We first examined the temporal evolution of the concentration–carbon feedback measured by changes in the total amount of terrestrial carbon (Fig. 1a; Fig. 2 of Arora et al. 2013) and then decomposed it into NPP and HR, following Eq. (5) (Figs. 1b,c). In the eight ESMs, elevated CO2 caused continuous carbon increase in terrestrial ecosystems (Fig. 1a), because in the BGC experiment, the elevated CO2 concentration had a “fertilization effect” on the terrestrial carbon cycle but did not produce significant global warming, which can reduce terrestrial carbon by intensifying ecosystem respiration. As expected, the elevated CO2 concentration also increased NPP (Fig. 1b). Because the stimulated productivity increased biomass (Table 2; ΔCV) and additional plant growth means more litter input to the soil, carbon accumulates in the soil (Table 2; ΔCS). Thus, in the BGC experiments, heterotrophic respiration also increased (Fig. 1c), because it is usually modeled as proportional to the size of the soil carbon pool. The ΔCS response was slow, producing a time lag after the stimulation of plant productivity by the elevated CO2 concentration. Thus, ΔHR was smaller than ΔNPP, ΔNEP was positive, and carbon accumulated in terrestrial ecosystems.

These trends in CL, NPP, and HR were shared by all the models, but model sensitivity to elevated CO2 concentration varied greatly among them (Table 2; Fig. 1). At the end of the simulation, ΔCL ranged from 186 PgC (NorESM1-ME) to 1206 PgC (MPI-ESM-LR), with a standard deviation (SD) of 380 PgC. The minimum–maximum ranges of ΔNPP and ΔHR among the ESMs were 87 and 82 PgC yr−1, respectively. These intermodel spreads are larger than the currently observed or estimated absolute value of NPP of 56.2 ± 14.3 PgC yr−1 (Ito 2011). Among the ESMs, the order in which ΔCL increased was very similar to that of increase in ΔNPP (Fig. 2a), and the rank correlation coefficient was 0.91. In other words, models with larger ΔNPP had greater ΔCL. Although HR sensitivity somewhat differed among the models because of the various HR schemes, ΔHR was strongly correlated (R2 = 0.99) with ΔNPP (Fig. 2b). Most carbon-based land ecosystem models do not account for nutrient uptake, and so the soil decomposition process cannot affect plant productivity. Therefore, the ΔHR in each ESM was strongly controlled by ΔNPP in this BGC experiment, not vice versa. Two models [NorESM1-ME and CESM1(BGC)] incorporate the nitrogen (N) cycle, and thus the size of their ΔNPP might be partly controlled by the size of ΔHR, because increased soil decomposition can accompany the nutrient release available for plants (e.g., Thornton et al. 2009). However, as shown in Fig. 2b, we did not find that these models exhibited clearly different behaviors relative to other nonnitrogen models. Thus, although these two models might have a feedback from ΔHR on ΔNPP to some extent, this does not seem to change the conclusion that the magnitude of ΔNPP strongly controls the size of ΔHR in this BGC simulation.

Fig. 2.

Relationship between simulated NPP increments and (a) ΔCL and (b) HR. Increments were calculated as the difference between averages of the last decade (years 131–140) and first decade (years 1–10). Lines in (a) for ΔCL were also calculated by an analytical method in which kL was assumed to depend (solid line) or not depend (dotted line) on NPP increments. In (b), regressed line passing through the origin is shown (y = 0.91x, with R2 = 0.99).

Fig. 2.

Relationship between simulated NPP increments and (a) ΔCL and (b) HR. Increments were calculated as the difference between averages of the last decade (years 131–140) and first decade (years 1–10). Lines in (a) for ΔCL were also calculated by an analytical method in which kL was assumed to depend (solid line) or not depend (dotted line) on NPP increments. In (b), regressed line passing through the origin is shown (y = 0.91x, with R2 = 0.99).

2) Intermodel spread found in decomposed factors

We evaluated the intermodel spread of ΔCL/ΔCO2 by examining the spread of its seven factors represented in Eq. (8b); the results are shown in Table 3 and Fig. 3. To emphasize the relative strength of intermodel spread across the factors, each variable in Fig. 3 is normalized by its multimodel mean; the normalized variable x′ is calculated as , where x is the original value, and is the multimodel mean (see also Fig. S1 in the supplemental material for the comparison of the decomposed factors with nonnormalized values). In Fig. 3, ΔCSCV, which represents the land carbon allocation between vegetation and soil, showed the largest spread, with normalized SD of 0.87. This was followed by factor Δgpp/ΔCO2 with slightly less SD (0.79) and ΔCV/ΔNPP (SD = 0.75). These were followed by factors with intermediate spread: ΔCL/ΔHR (SD = 0.53); ΔHR/ΔCS (SD = 0.44); and ΔGPP/Δgpp (SD = 0.40). The smallest spread was for ΔNPP/ΔGPP, which relates to plant respiration, with SD = 0.38. The decomposed factors represented by the absolute magnitude (Table 3) showed the most factors at the end of simulation change from their initial state. This means that, even in biogeochemically coupled simulations where only CO2 increases cause atmospheric environmental changes of land ecosystems, CO2 stimulates various processes in the terrestrial carbon cycle. The extents of the responses are varied among models.

Table 3.

Comparison of ΔCL/ΔCO2 (=β) and decomposed factors with absolute values. Decomposed factors for start are calculated by , where and are 10-yr averages at the beginning of 1-ppy BGC simulation for variable and . For end, the variables are calculated by , and and represent 10-yr average at the end of simulation. For all, the variables are calculated by , which follows the definition shown in Eq. (8b).

Comparison of ΔCL/ΔCO2 (=β) and decomposed factors with absolute values. Decomposed factors for start are calculated by , where  and  are 10-yr averages at the beginning of 1-ppy BGC simulation for variable  and . For end, the variables are calculated by , and  and  represent 10-yr average at the end of simulation. For all, the variables are calculated by , which follows the definition shown in Eq. (8b).
Comparison of ΔCL/ΔCO2 (=β) and decomposed factors with absolute values. Decomposed factors for start are calculated by , where  and  are 10-yr averages at the beginning of 1-ppy BGC simulation for variable  and . For end, the variables are calculated by , and  and  represent 10-yr average at the end of simulation. For all, the variables are calculated by , which follows the definition shown in Eq. (8b).
Fig. 3.

Intermodel spread of ΔCL/ΔCO2 (=β) and its decomposed factors, consisting of Δgpp/ΔCO2, ΔGPP/Δgpp, ΔNPP/ΔGPP, ΔCV/ΔNPP, ΔCSCV, ΔHR/ΔCS, and ΔCL/ΔHR. Variables shown in panels were normalized by multimodel means [, where x, , and x′ represent original value, multimodel mean, and normalized value, respectively]. Crosses represent nonnitrogen cycle models and circles represent carbon–nitrogen cycle models. Gray shades show SD normalized by multimodel means. Increments of each factor were calculated as the difference between averages of the last decade (years 131–140) and first decade (years 1–10).

Fig. 3.

Intermodel spread of ΔCL/ΔCO2 (=β) and its decomposed factors, consisting of Δgpp/ΔCO2, ΔGPP/Δgpp, ΔNPP/ΔGPP, ΔCV/ΔNPP, ΔCSCV, ΔHR/ΔCS, and ΔCL/ΔHR. Variables shown in panels were normalized by multimodel means [, where x, , and x′ represent original value, multimodel mean, and normalized value, respectively]. Crosses represent nonnitrogen cycle models and circles represent carbon–nitrogen cycle models. Gray shades show SD normalized by multimodel means. Increments of each factor were calculated as the difference between averages of the last decade (years 131–140) and first decade (years 1–10).

In addition to the distribution of intermodel spreads among factors, this analysis provided information on the relative location of each ESM. For example, NorESM1-ME and CESM1(BGC), which include an explicit nitrogen cycle that may greatly alter land carbon cycle dynamics, produced distinctly lower values of Δgpp/ΔCO2. Since the N cycle in these models acts to constrain the potential GPP, depending on available nitrogen for plants (Tjiputra et al. 2013), the resultant concentration–carbon feedback can be significantly reduced by N limitation (Bonan and Levis 2010; Thornton et al. 2009). Thus, the smaller Δgpp/ΔCO2 values from NorESM-M and CESM1(BGC) may be attributable to N limitation of the carbon cycle. Although these two models also gave smaller values of ΔNPP/ΔGPP, it is not likely that this is related to the inclusion of the nitrogen cycle, because smaller values of ΔNPP/ΔGPP were also evident in the “carbon-only” simulations (Gotangco Castillo et al. 2012).

The implementation of dynamic global vegetation models (DGVMs) in land ecosystem models may also produce the large difference in Fig. 3. In our research, three ESMs (HadGEM2-ES, MPI-ESM-LR, and MIROC-ESM) used DGVMs, but we can neither find a systematic trend that is likely derived from DGVM implementation nor any reasonable DGVM-specific mechanism that can cause variation of the decomposed factors.

3) Effective parameters in STC model

In the previous section, we demonstrated the intermodel spread among decomposed factors, and the analysis using Eq. (8b) successfully depicted what occurs in carbon cycle processes in ESMs. However, it should be noted that a factor with a greater intermodel spread does not necessarily indicate the primary factor that determines the strength of concentration–carbon feedback. As described in section 2, this is because the decomposed factors in Eq. (8b) sometimes depend on each other; changes in one factor may alter other factors located “downstream.” To solve this problem, we used the STC model composed by Eqs. (9)(16) as follows. First, this model was fit to each ESM’s behavior in the BGC simulation, by automatically selecting parameters that minimize errors in GPP for kGPP1 and kGPP2, AR for kAR, LF for kLF, and HR for kHR. The fitted parameters are shown in Table 4. Next, to find the parameters that control the concentration–carbon feedback strength, we conducted BGC simulations using the STC model, increasing the number of ESM-derived parameters to be applied to the STC model in a stepwise fashion.

Table 4.

Parameters obtained by fitting STC model against ESMs and estimation biases by STC models. Parameters are obtained by fitting formulas of STC model represented by Eqs. (9)(14) to automatically minimize the squared biases. Correlation coefficient R is evaluated by comparing variables (ΔGPP, ΔAR, ΔLF, and ΔHR) from ESMs and STC model in 140-yr, 1-ppy BGC simulation.

Parameters obtained by fitting STC model against ESMs and estimation biases by STC models. Parameters are obtained by fitting formulas of STC model represented by Eqs. (9)–(14) to automatically minimize the squared biases. Correlation coefficient R is evaluated by comparing variables (ΔGPP, ΔAR, ΔLF, and ΔHR) from ESMs and STC model in 140-yr, 1-ppy BGC simulation.
Parameters obtained by fitting STC model against ESMs and estimation biases by STC models. Parameters are obtained by fitting formulas of STC model represented by Eqs. (9)–(14) to automatically minimize the squared biases. Correlation coefficient R is evaluated by comparing variables (ΔGPP, ΔAR, ΔLF, and ΔHR) from ESMs and STC model in 140-yr, 1-ppy BGC simulation.

Figure 4a shows that, although only the GPP-related parameters (kGPP1 and kGPP2) of the STC model are tuned to each ESM and the other parameters (kAR and kHR) are obtained using multimodel means, the reproducibility of the concentration–carbon feedback strength was reasonably good. It had a correlation coefficient R of 0.92 and root-mean-square error (RMSE) of 0.22. By further applying ESM-specific parameters to AR (kAR) in the STC model, the reproducibility was slightly improved, with 0.95 for R and 0.17 for RMSE (Fig. 4b). Further application of ESM-specific parameters to LF (Fig. 4c) also demonstrates a slight improvement in the STC model capability, with 0.95 for R and 0.17 for RMSE. Finally, all ESM-specific parameters were used in the model (Fig. 4d), resulting in very good reproducibility with R = 0.98 and RMSE = 0.086. It is possible, however, that the correlation coefficient can depend on the sequence that ESM-specific parameters are added to the STC model. To confirm this, the reproducibility of the STC model was further examined by only applying the (i) GPP-related parameters (kGPP1 and kGPP2), (ii) AR-related parameter (kAR), (iii) LF-related parameter (kLF), or (iv) HR-related parameter (kHR). From this analysis, we confirmed that the GPP-related parameters have the biggest influence on the STC models ability to reproduce the concentration–carbon feedback of ESMs, with the highest correlation coefficient (0.92). Separately including the ESM-derived parameters of kAR, kLF, or kHR resulted in far smaller correlation coefficients of 0.47, −0.77, and 0.15, respectively (Fig. S2 in the supplemental material).

Fig. 4.

Comparison of ESMs and STC models for β (PgC ppmv−1) in 140-yr simulations of BGC experiment (symbols: colors the same as in Fig. 3). STC model was fit against each ESM (a) only for GPP-related parameters (kGPP1 and kGPP2); (b) GPP- and AR-related parameters (kGPP1, kGPP2, and kAR); (c) GPP-, AR-, and LF-related parameters (kGPP1, kGPP2, kAR, and kLF); and (d) GPP-, AR-, LF-, and HR-related parameters (kGPP1, kGPP2, kAR, kLF, and kHR). Black line represents regression line of y = ax + b and gray line is y = x. RMSE was calculated by , where N (=8) is total number of models.

Fig. 4.

Comparison of ESMs and STC models for β (PgC ppmv−1) in 140-yr simulations of BGC experiment (symbols: colors the same as in Fig. 3). STC model was fit against each ESM (a) only for GPP-related parameters (kGPP1 and kGPP2); (b) GPP- and AR-related parameters (kGPP1, kGPP2, and kAR); (c) GPP-, AR-, and LF-related parameters (kGPP1, kGPP2, kAR, and kLF); and (d) GPP-, AR-, LF-, and HR-related parameters (kGPP1, kGPP2, kAR, kLF, and kHR). Black line represents regression line of y = ax + b and gray line is y = x. RMSE was calculated by , where N (=8) is total number of models.

This analysis clearly shows that the GPP response to elevated CO2 concentration is predominant in controlling the strength of concentration–carbon feedback in ESMs, and the contributions of other processes to the feedback strength are likely smaller. These results are consistent with that shown in Fig. 2; the response of plant productivity to elevated CO2 concentration may be the major control on concentration–carbon feedback strength in the ESMs. This is supported by Tachiiri et al. (2012), who stated that parameters related to photosynthetic capacity may be one of the most influential in determining the strength of concentration–carbon feedback.

b. Concentration–carbon feedbacks among CO2 scenarios

1) Land carbon and β among scenarios: MIROC-ESM and STC model

To examine concentration–carbon feedback in various scenarios, we carried out simulations with two groups of scenarios (TRN and STB) and analyzed CL and β in those scenarios (Figs. 5a–d for TRN and Figs. 5e–h for STB). We used the MIROC-ESM and STC model with parameters tuned to MIROC-ESM for these simulations. The STC model was set with three global carbon pools ( appendix A). The reason for using these pools is described later. In TRN simulations, as CO2 increases more rapidly at 2.0 ppy and the CO2 fertilization effect strengthens, ΔCL shows a faster increase than in the other two scenarios (Figs. 5a,b). In contrast, as CO2 slowly increases in a 0.5-ppy simulation, ΔCL shows the slowest increase. However, when ΔCL was plotted against CO2 concentration, the 0.5-ppy simulation showed higher carbon accumulation relative to the other scenarios (Fig. 5c). The reason for this in the slower rate of the CO2 increase scenario can be explained by the degree to which terrestrial carbon in the transient state is close to the equilibrium carbon amount at a given CO2 concentration. To confirm this, the equilibrium carbon amount at a given CO2 level was estimated by extrapolating the STC model simulation to an equilibrium state (gray line in Fig. 5 and see  appendix A for details). The equilibrium carbon amount always showed a larger ΔCL than those of other TRN simulations across all CO2 concentration ranges, and the 0.5-ppy simulation was relatively close to that of the equilibrium state. These results clearly indicated that land carbon can be effectively accumulated when CO2 concentration slowly increases, because the state of carbon in a slower scenario is closer to the equilibrium carbon amount at a given CO2 level. This mechanism produces the interscenario spread of β, producing a larger β in slower rate of CO2 increase scenarios, as shown by Gregory et al. (2009).

Fig. 5.

Terrestrial carbon change and concentration–carbon feedback parameter in multiscenario simulations conducted by MIROC-ESM (crosses) and a simple model with three carbon pools (lines). (top) The ESM and the simple model were forced with (a) monotonically CO2 increase scenario groups from TRN. The scenarios consist of 2.0-ppy CO2 increase from preindustrial levels (CO2PI) for 70 yr (red), 1.0 ppy for 140 yr (black), and 0.5 ppy for 280 yr (blue), in which CO2 levels ultimately reach about 1140 ppmv (4 × CO2PI) in all simulations. (bottom) The ESM and simple model were forced with the STB scenario group, consisting of (e) four CO2 scenarios: abrupt CO2 doubling from preindustrial state and stabilized (green), stabilized scenario at 2 × CO2PI level following a 2.0-ppy (red), 1.0-ppy (black), and 0.5-ppy (blue) CO2 scenario, in which total simulation years by the ESM and simple model were set to 140 and 400 yr, respectively. Equilibrium carbon state and concentration–carbon feedback parameter at any CO2 level were estimated by the simple model and shown by gray lines. For simplified mode, the parameters for vegetation are shown in Table 4, and soil parameters were manually fitted against MIROC-ESM and given as kSF = 0.25, kSS = 0.012, and fA = 0.7.

Fig. 5.

Terrestrial carbon change and concentration–carbon feedback parameter in multiscenario simulations conducted by MIROC-ESM (crosses) and a simple model with three carbon pools (lines). (top) The ESM and the simple model were forced with (a) monotonically CO2 increase scenario groups from TRN. The scenarios consist of 2.0-ppy CO2 increase from preindustrial levels (CO2PI) for 70 yr (red), 1.0 ppy for 140 yr (black), and 0.5 ppy for 280 yr (blue), in which CO2 levels ultimately reach about 1140 ppmv (4 × CO2PI) in all simulations. (bottom) The ESM and simple model were forced with the STB scenario group, consisting of (e) four CO2 scenarios: abrupt CO2 doubling from preindustrial state and stabilized (green), stabilized scenario at 2 × CO2PI level following a 2.0-ppy (red), 1.0-ppy (black), and 0.5-ppy (blue) CO2 scenario, in which total simulation years by the ESM and simple model were set to 140 and 400 yr, respectively. Equilibrium carbon state and concentration–carbon feedback parameter at any CO2 level were estimated by the simple model and shown by gray lines. For simplified mode, the parameters for vegetation are shown in Table 4, and soil parameters were manually fitted against MIROC-ESM and given as kSF = 0.25, kSS = 0.012, and fA = 0.7.

In STB simulations, the spread of ΔCL among scenarios simulated by MIROC-ESM expanded in the first half of the 140-yr simulations (Fig. 5f). However, as CO2 concentration stabilized after reaching 2 × CO2PI, the interscenario spread of ΔCL decreased. In fact, extended simulations by the STC model out to 400 yr showed that the spread ultimately diminished. As a result, when we examined the relationships of ΔCL–CO2 or β–CO2 in STB simulations, ΔCL in all scenarios gradually approached its equilibrium state after CO2 concentration stabilized (Figs. 5g,h).

Here, we selected the STC model with the three global carbon pools. This is because the STC model with two carbon pools had difficulty in capturing nonlinear behaviors of a carbon cycle with more complex CO2 pathways (Fig. S3 in the supplemental material), most likely since it was unable to represent the transitional allocation change among soil carbon pools. However, as shown in section 4a(2), two carbon pools are adequate and effective for capturing multi-ESMs behaviors in a single scenario, especially in one with a simple CO2 pathway such as the 1-ppy CO2 increase simulation.

2) Interscenario spread in decomposed factors

In the previous section, we described how concentration–carbon feedback changes under different scenarios, yielding different values of ΔCL or β. To discover what occurs within the terrestrial carbon cycle, we applied the identity in Eq. (8) to the multiscenario simulations. Here, we selected TRN simulations by MIROC-ESM for the analyses. The differential periods for the analysis were 0–280 yr for the 0.5-ppy scenario simulation, 0–140 yr for 1.0 ppy, and 0–70 yr for 2.0 ppy (i.e., final CO2 level was identical between scenarios). The results are depicted in Fig. 6 and show that ΔCL/ΔCO2 (=β) was higher in a scenario where the rate of CO2 increase is slower. However, the first four decomposed factors (Δgpp/ΔCO2, ΔGPP/Δgpp, ΔNPP/ΔGPP, and ΔCV/ΔNPP), which consist of carbon flux or pool with high turnover rate (ΔCV), did not change much between scenarios. This suggests that these processes simply follow the CO2 concentration at a given time, with less time lag. In contrast, the last three factors (ΔCSCV, ΔHR/ΔCS, and ΔCL/ΔHR), which include variables related to carbon pools with a slow turnover rate (ΔCS or ΔCL), showed clear spreads among scenarios; more carbon was allocated to soils in slower scenarios (higher value of ΔCSCV), along with a low decomposition rate (lower value of ΔHR/ΔCS) and longer residence time of carbon in the land ecosystem (higher ΔCL/ΔHR). Because those carbon pools with low turnover rate slowly respond to CO2 increases, these carbon pools lag the forcing and consequently produce the interscenario spread.

Fig. 6.

As in Fig. 3, but for results of BGC simulation by MIROC-ESM under TRN scenarios. TRN simulation consists of 140-yr simulation with 1.0-ppy CO2 increase scenario (crosses), 280 yr for 0.5 ppy (circles), and 70 yr for 2.0 ppy (squares). In the three scenarios, CO2 concentration ultimately reached about 1140 ppmv.

Fig. 6.

As in Fig. 3, but for results of BGC simulation by MIROC-ESM under TRN scenarios. TRN simulation consists of 140-yr simulation with 1.0-ppy CO2 increase scenario (crosses), 280 yr for 0.5 ppy (circles), and 70 yr for 2.0 ppy (squares). In the three scenarios, CO2 concentration ultimately reached about 1140 ppmv.

From this analysis, we found that several carbon cycle processes within the ecosystems can depend on scenarios, especially in the processes with slow responses to CO2 increases. This is analogous to the conclusion in the previous section that total carbon storage change depends on scenarios because it has a time lag to CO2 increases. Here, we note two points. First, considering that the vegetation pool in MIROC-ESM has a relatively quick response to forcing (confirmed by relatively lower value in ΔCV/ΔNPP in Table 3 and higher value in kLF in Table 4), the spread found in ΔCV/ΔNPP in Fig. 6 is likely to be underestimated, and thus the spread in ΔCSCV might be overestimated. Second, Fig. 6 displays the carbon cycle behavior in transient states. The varying responses of factors between scenarios (e.g., ΔCSCV) gradually diminished after the CO2 concentration stabilized (such as in STB simulations; data not shown). This leads to several carbon processes changing from their initial condition.

c. Concentration–carbon feedback in RCP simulations

We have described a series of idealized simulations designed to conduct sensitivity analysis. However, in multimodel simulations with multiple scenarios, such as those shown by Jones et al. (2013) in which more complex CO2 pathways in RCP scenarios (Fig. 7a) were used for ESMs, the impacts of both interscenario and intermodel differences on the simulation results must be considered simultaneously. In the MIROC-ESM, concentration–carbon feedback at the end of the TRN simulation was 0.55 PgC ppmv−1 in the 2.0-ppy simulation, 0.74 PgC ppmv−1 for 1.0 ppy, and 0.93 PgC ppmv−1 for 0.5 ppy (black symbols in Fig. 7b). Thus, the magnitude of β in the MIROC-ESM increased about 25% upon halving the CO2 increase rate and decreased about 25% upon doubling that rate. Because the CO2 growth rates in these idealized scenarios were higher than those of the RCPs (Fig. 7a), β should be larger in the RCP simulations than in our sensitivity analysis simulations (β for the RCP simulation defined here assumes a situation in which only CO2 forcing acts on the carbon cycle). In a BGC simulation by MIROC-ESM in which the CO2 pathways of RCP4.5 was applied and other forcings (such as the land use fraction) were fixed at the preindustrial state, β (filled red circle in Fig. 7b) was about twice that in the 0.5-ppy scenario, as expected. Furthermore, by assuming that the sensitivity of climate–carbon feedback (parameter γ) is constant throughout the simulations and independent of scenario, we can approximate β in a fully coupled experiment [open symbols in Fig. 7b; see  appendix B for details and Hajima et al. (2012)]. In simulations with MIROC-ESM, this indirect estimation of β appeared to work well, and β showed an exponentially decreasing trend versus the rate of CO2 increase. Scenario differences produced a spread in β comparable in magnitude to that in simulations by multiple ESMs (gray symbols in Fig. 7b). Considering that the simplified STC model also displays a dependence of concentration–carbon feedback on CO2 growth rate (Fig. 5) and that the basic structure of that model is shared with ESMs, other ESMs are likely to give larger values of β in RCP scenarios than those in 1-ppy CO2 increase scenarios.

Fig. 7.

(a) CO2 concentration trajectories for a historical simulation (gray), the four RCPs (colors), and TRN simulations with CO2 increase rates 0.5, 1.0, and 2.0 ppy (dotted black, dashed black, and solid black lines, respectively). (b) Concentration–carbon feedback parameter β plotted against rate of CO2 increase in different scenarios. Rate of each scenario was calculated by assuming that CO2 increased linearly toward the CO2 level at the end of each simulation. Filled symbols represent β obtained from BGC simulations with MIROC-ESM. Open symbols represent β as approximated by fully coupled simulations by MIROC-ESM with assumption that γ (climate–carbon feedback parameter) was constant and independent of scenario [β = (ΔCLγΔT)/ΔCA; see  appendix B for derivations]. Data required to approximate γ, ΔCL, and ΔCA were from Hajima et al. (2012). Filled gray squares represent values of β obtained from CMIP5 ESMs, except for MIROC-ESM, in a 1-ppy BGC simulation. Curve was fit against MIROC-ESM results (open and closed symbols, excluding gray squares) and estimated as y = 1.57x−0.45 (R2 = 0.97). All β shown in this figure were calculated by taking the difference of terrestrial carbon amount between the beginning and end of each simulation.

Fig. 7.

(a) CO2 concentration trajectories for a historical simulation (gray), the four RCPs (colors), and TRN simulations with CO2 increase rates 0.5, 1.0, and 2.0 ppy (dotted black, dashed black, and solid black lines, respectively). (b) Concentration–carbon feedback parameter β plotted against rate of CO2 increase in different scenarios. Rate of each scenario was calculated by assuming that CO2 increased linearly toward the CO2 level at the end of each simulation. Filled symbols represent β obtained from BGC simulations with MIROC-ESM. Open symbols represent β as approximated by fully coupled simulations by MIROC-ESM with assumption that γ (climate–carbon feedback parameter) was constant and independent of scenario [β = (ΔCLγΔT)/ΔCA; see  appendix B for derivations]. Data required to approximate γ, ΔCL, and ΔCA were from Hajima et al. (2012). Filled gray squares represent values of β obtained from CMIP5 ESMs, except for MIROC-ESM, in a 1-ppy BGC simulation. Curve was fit against MIROC-ESM results (open and closed symbols, excluding gray squares) and estimated as y = 1.57x−0.45 (R2 = 0.97). All β shown in this figure were calculated by taking the difference of terrestrial carbon amount between the beginning and end of each simulation.

We note that anthropogenic non-CO2 forcings included in RCP scenarios, such as land cover, greatly impact the land carbon cycle, and so those forcings can substantially change the β value (Jones et al. 2013). In contrast to this “apparent β” that reflects impacts from both CO2 concentration and other non-CO2 forcings, the β shown in Fig. 7 is the concentration–carbon feedback parameter in the BGC experiment, in which only CO2 forcing acts on the terrestrial carbon cycle. Thus, if we consider both impacts from CO2 and non-CO2 forcing in the simulation, the apparent β can be greatly changed from the β shown in this study. It is also noted that the β shown in Fig. 7 represents the concentration–carbon feedback parameter in the transient state. As shown by STB simulations (Fig. 5), if CO2 concentration is stabilized at a certain level, β should settle in a different state.

5. Discussion

a. Analytical determination of concentration–carbon feedback

So far, we have used numerical results to infer the behavior of concentration–carbon feedback. Here, we investigate the nonlinear and nonintuitive behavior of this feedback by analytically solving the global carbon budget, as used in Friedlingstein et al. (1995) or Thompson et al. (1996).

First, to solve analytically, we assume that ΔNPP is proportional to time:

 
formula

here R is the CO2 increase rate of the scenario (ppmv yr−1), S is the NPP sensitivity to elevated CO2 concentration (PgC yr−1 ppmv−1), t is time (yr), and NPP0 is the initial value of NPP. A simplification of the NPP response to time in Eq. (17) was assumed in some studies, such as Friedlingstein et al. (1995) and Thompson et al. (1996). The former work showed that this linear assumption for that response does not qualitatively change the behavior of terrestrial carbon uptake.

Next, for simplicity, we treat the global carbon pool with one component CL. Thus, the global carbon balance becomes

 
formula

where carbon flux out of the biosphere is assumed proportional to CL, with a fixed parameter kL (in yr−1).

By substituting Eq. (17) into Eq. (18), we obtain

 
formula

Because this formula has the form of a first-order linear ordinary differential equation (dy/dx = axby), its analytical solution is as follows (see  appendix C for derivation):

 
formula

Changes in CL after T years are represented as

 
formula

We can use Eq. (20) to examine whether slower CO2 increase scenarios effectively accumulate carbon. First, we assume ΔCL in two simulations: ΔCL,1 and ΔCL,2 in which CO2 growth rates are R1 and R2 and integration times are T1 and T2, respectively. Second, by taking the difference of ΔCL,1 and ΔCL,2, we obtain

 
formula

Finally, assuming that CL,2 is a result of a simulation in which the CO2 growth rate is halved (R2 = 0.5R1) and the integration time is doubled (T2 = 2T1), the difference should be

 
formula

(see  appendix C for detailed derivations). From Eq. (21), it is clear that more carbon accumulates in a slow scenario than in a faster one. The difference between the two is proportional to NPP sensitivity to elevated CO2 S and inversely proportional to the square of the decomposition rate kL2, with the T-dependent factor .

We can use Eq. (20) to verify the reason for the nonlinear form of the ΔNPP–ΔCL relationship between models (Fig. 2a). In the BGC simulations by multiple ESMs with a common scenario, integration time T and CO2 growth rate R are constant and shared among ESMs. Thus, cumulative changes in terrestrial carbon in models with the same integration time and CO2 increase rate scenario can be summarized as follows (see  appendix C):

 
formula

From this equation, it is clear that ΔCL in each ESM is proportional to ΔNPP or S, but this proportionality is violated if kL varies systematically with ΔNPP in the models. In fact, the parameter shows a trend in which kL is greater in ESMs with large ΔNPP, showing a relationship (R2 = 0.88) of

 
formula

(see Fig. S4 in the supplemental material for the scatterplot). By substituting this empirical relationship into Eq. (22), we obtain a continuous function that relates ΔCL with ΔNPP, which as expected successfully captures the nonlinear trend of concentration–carbon feedback (solid line in Fig. 2a). If kL is instead assumed constant (with kL = 0.02 assumed), then the relationship is no longer nonlinear (dotted line in Fig. 2a). Thus, the results of this analytical approach mean that ESMs with large ΔNPP can have a smaller concentration–carbon feedback than expected from their large productivity, because kL is systematically greater in models with larger ΔNPP. It is noted that the analytical solutions shown in Eqs. (20)(22) hold in a transient situation, such as TRN simulations.

b. Suggestions to reduce uncertainty of terrestrial concentration–carbon feedback

Our results clearly show that the magnitude of the concentration–carbon feedbacks reproduced by multiple ESMs are strongly regulated by NPP sensitivity to elevated CO2 concentration and are weakly affected by different modeling schemes and parameters related to soil decomposition. Therefore, to reduce the intermodel spread of that feedback, the important issue is how to improve the predictability of the global ΔNPP. As described in Eq. (8b), ΔNPP in ESMs can be summarized by three factors: 1) photosynthetic sensitivity to elevated CO2 concentration, which is related in our analysis to the factor Δgpp/ΔCO2; 2) the total mass or area of assimilation organ (leaf), which is implied here by the term ΔGPP/Δgpp (this term is identical to ΔLAI because we defined Δgpp = ΔGPP/ΔLAI); and 3) the autotrophic respiration rate represented by ΔNPP/ΔGPP. In the present work, ΔNPP/ΔGPP was the factor with the smallest spread among models, whereas Δgpp/ΔCO2 produced one of the largest spreads. Thus, an accurate estimation of ΔNPP depends on the capability of the photosynthesis response to increasing CO2 concentration, which is summarized by the factor Δgpp/ΔCO2.

It is interesting that the factor Δgpp/ΔCO2 had negative correlation with factor ΔGPP/Δgpp (R = −0.79). This means that ΔGPP/Δgpp, which is linked to ΔLAI, had a lower value in models that had a large Δgpp/ΔCO2. These opposing trends of the two factors suggest that large Δgpp/ΔCO2 in some models might be partly balanced by their small ΔLAI and vice versa. This fact suggests that CMIP5 ESMs increase their GPP in a different manner: by increasing LAI or by being highly sensitive to CO2 in photosynthetic processes. Thus, to reduce the spread of GPP responses to elevated CO2 found in multi-ESMs simulations, it is necessary to give simultaneous constraints to the LAI response (ΔGPP/Δgpp), the sensitivity in non-LAI-related processes (Δgpp/ΔCO2), and the overall sensitivity of plant productivity to elevated CO2 concentration (ΔGPP/ΔCO2 = Δgpp/ΔCO2 × ΔGPP/Δgpp). For this purpose, it is important that the modeling of the relationship between photosynthesis, LAI, and plant productivity be supported by more evidence from free-air CO2 enrichment experiments (FACE; e.g., Ainsworth and Long 2005; Norby and Zak 2011) and a growing number of inventory experiments (e.g., Lewis et al. 2009).

We have shown herein that the strength of the concentration–carbon feedback can vary with the sensitivity of plant productivity to elevated CO2 concentration, increased CO2 rate, and the integration time used in the simulation. Because these factors alter the strength of the terrestrial concentration–carbon feedback determined by simulation, care must be taken when analyzing and interpreting projected climate or related variables obtained by multimodel simulations with several scenarios. For example, Jones et al. (2013) compared compatible fossil fuel emissions between ESMs and integrated assessment models (IAMs) using four RCPs. They found systematic discrepancies in the fossil fuel emissions between those models, with ESMs tending to simulate lesser emissions than IAMs, especially with high-end RCP scenarios. Considering that concentration–carbon feedback by land and ocean could have greater impact than climate–carbon feedback (as shown by Arora et al. 2013), and considering that concentration–carbon feedback can be substantially reduced with an higher rate of CO2 increase, the systematic discrepancies between ESMs and IAMs in RCP simulations may be attributable to features of concentration–carbon feedback in land and ocean.

6. Summary and conclusions

This research represents a first attempt at comprehensively investigating the underlying mechanisms and uncertain processes of concentration–carbon feedback in the terrestrial carbon cycle, using outcomes of CMIP5 ESMs simulations. We evaluated the strength of the concentration–carbon feedback in terrestrial ecosystems among the ESMs using a 1-ppy CO2 increment scenario and biogeochemically coupled settings, focusing on changes in the amount of terrestrial carbon. By examining the changes in the carbon pools and fluxes of the models in detail, we demonstrated that CO2 increases stimulate several carbon cycle processes (such as plant production, litter fall, and heterotrophic respiration), and the degree of the responses are different among the models. However, the large spread in concentration–carbon feedback among the ESMs was explained by the sensitivity of plant productivity to the elevated CO2 concentration reproduced by each model. Thus, to reduce this spread in the strength of the concentration–carbon feedback among ESMs, more investigations are required into plant production sensitivity to elevated CO2 concentration. Because the plant productivity responses to CO2 increases depend on the response of both photosynthesis and leaf area, it is also necessary to simultaneously constrain the model’s behavior on these two processes.

By conducting simulations with an ESM and simplified global carbon models using several types of CO2 pathways, we confirmed that higher carbon accumulation is achieved in scenarios with slower CO2 increase rates for a given CO2 concentration. Using numerical and analytical analyses, we clearly demonstrated that a time lag of terrestrial carbon with atmospheric CO2 increase is the crucial process for the large interscenario spread of concentration–carbon feedback; a high rate of CO2 increase makes the terrestrial carbon amount much less than the equilibrium carbon amount for a given CO2 concentration. Similarly, a time lag to CO2 forcing is also found in the terrestrial ecosystem, especially in the processes related to carbon pools with lower turnover rate (soil or total terrestrial carbon pools). These time lags result in different magnitudes of the feedback strength quantified by parameter β: the multiscenario spread of β in a single ESM across a broad range of prescribed CO2 concentration pathways is comparable with the spread between different ESMs for a single-concentration scenario. This result strengthens the conclusions of previous studies that quantities for evaluating the feedback (e.g., parameter β) are incapable of capturing the land carbon storage state and has direct policy implications, in that the concentration feedback parameter β has a larger value for future scenarios where emissions are reduced.

Acknowledgments

We acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for CMIP5, and we thank the climate modeling groups for producing and making available their model output. We greatly appreciate helpful comments from Drs. Vivek Arora and Atsuhiro Iio. We also thank our editor and the three reviewers, who encouraged us and gave constructive comments. This work was supported by the Program for Risk Information on Climate Change, MEXT, Japan, and by a Grant-in-Aid for Scientific Research on Innovative Areas (21114010).

APPENDIX A

Global Carbon Cycle with Three Carbon Pools

If the carbon pool is divided into three components, a vegetation pool CV and soil pools of fast CSF and slow CSS turnover rates, carbon allocation in the terrestrial ecosystem is expressed as

 
formula

Carbon balance in the vegetation pool can be described by

 
formula

Carbon balance in the fast pool can be described by

 
formula

and

 
formula

where DSF represents the total carbon flow out of the fast pool. This flow is composed of two fluxes: efflux by respiration and flow to the slow carbon pool. By defining fA as the ratio of respiration efflux to total efflux DSF, we obtain

 
formula

Thus, (1 − fA)DSF must be the influx into the slow pool. Similarly, the carbon balance equation for the slow pool is

 
formula

and

 
formula

Here, all carbon decomposed in the slow pool is assumed to be consumed by respiration:

 
formula

Thus, total heterotrophic respiration from the terrestrial carbon pools is

 
formula

By solving these equations, we can compute the global carbon cycle with three carbon pools. For the numerical computation, we need initial values of CV, CSF, and CSS. Because under steady-state conditions, the left-hand sides of Eqs. (A1)(A3) must equal zero, we obtain initial values for the three pools:

 
formula
 
formula
 
formula

In section 4b, to mimic the behavior of MIROC-ESM, the parameter for vegetation in Table 4 (MIROC-ESM) is applied to the STC model, but soil parameters for the simulations are manually fit against the STB and TRN simulation results from MIROC-ESM: kSF = 0.25, kSS = 0.012, and fA = 0.7.

Similarly, we can estimate the carbon amounts in a new steady state after CO2 concentration is perturbed. By representing GPP in the new steady state as GPP′, the carbon amount for three carbon pools in the state (, , and ) can be written as

 
formula
 
formula
 
formula

Consequently, the total carbon amount of land ecosystem can be represented as

 
formula

By substituting Eq. (10) into this, becomes a function of :

 
formula

and changes from the initial condition should be

 
formula

In this equation, quickly increases when is relatively small and gradually saturates as becomes larger, following the plant productivity response to CO2 concentration represented by the Michaelis–Menten function. Consequently, the concentration–carbon feedback parameter in the new equilibrium condition () should be larger for small and gradually decrease with larger . That is, in the following equation,

 
formula

where approaches a finite value when (but should not be 0 because of the definition of ), and becomes smaller as becomes larger.

APPENDIX B

Approximating β from a Fully Coupled Experiment and γ

Friedlingstein et al. (2003) identified two major carbon cycle feedbacks, concentration–carbon and climate–carbon. They expressed terrestrial carbon change as a function of these carbon cycle feedbacks:

 
formula

where CL and CA represent change in terrestrial and atmospheric carbon, respectively; β is the concentration–carbon feedback parameter (PgC ppmv−1); and γ is the climate–carbon feedback parameter (PgC K−1). By rearranging this equation, β can be expressed as

 
formula

Because interscenario variation in γ is smaller than that in β as shown by Gregory et al. (2009), we can approximate β in a fully coupled experiment by fixing the value of γ, and using for other variables (ΔCA, ΔCL, and ΔT) values obtained from a fully coupled experiment using RCPs.

APPENDIX C

Solving the Global Carbon Balance Equation Analytically

a. Derivation of the analytical solution for CL

The global carbon balance with one carbon pool, in which NPP and HR are respectively regarded as influx and efflux of the carbon pool, can be described as

 
formula

For simplicity, we assume that NPP responds linearly to time t (yr) for a scenario with a fixed CO2 increase rate R (ppmv yr−1) and NPP sensitivity to elevated CO2 concentration S (PgC yr−1 ppmv−1):

 
formula

By substituting Eq. (C2) for Eq. (C1), we obtain

 
formula

By replacing CL, t, RS, NPP0, and k with y, x, a, a0, and b, respectively, it is clear that this formula has the form of a first-order ordinary differential equation:

 
formula

The analytical solution to this equation is

 
formula

where c0 is an arbitrary constant. Thus, we obtain the following analytical solution for Eq. (C3):

 
formula

Considering that CL should be in a steady state at t = 0, the initial value of CL, derived by using Eq. (C3), is NPP0/k. By comparing this result with Eq. (C4) for t = 0, we find that c0 = RS/kL2 and Eq. (C4) becomes

 
formula
b. CL difference between slow and fast scenarios

We can use Eq. (C5) to confirm that CL in a slower scenario can be larger than CL in a faster scenario. Changes of terrestrial carbon in a scenario (ΔCL,1), in which integration time is T1 and CO2 increase rate is R1, can be represented as

 
formula

If we assume a slower scenario in which the increase rate of CO2 is halved (R2 = R1/2) and integration time is doubled (T2 = 2T1), we replace R with 0.5R and t with 2t in Eq. (C5):

 
formula

By taking the difference of ΔCL,1 and ΔCL,2, we obtain

 
formula
c. Nonlinear behavior of CL

We can use Eq. (C5) to calculate the terrestrial carbon change in a BGC simulation as

 
formula

where t0 is the initial time, assumed here as 0. Thus, this equation becomes

 
formula

and because there is a relationship ,

 
formula

REFERENCES

REFERENCES
Ainsworth
,
E. A.
, and
S. P.
Long
,
2005
:
What have we learned from 15 years of free-air CO2 enrichment (FACE)? A meta-analytic review of the responses of photosynthesis, canopy properties and plant production to rising CO2
.
New Phytol.
,
165
,
351
372
,
doi:10.1111/j.1469-8137.2004.01224.x
.
Anav
,
A.
, and Coauthors
,
2013
:
Evaluating the land and ocean components of the global carbon cycle in the CMIP5 Earth system models
.
J. Climate
, 26, 6801–6843, doi:10.1175/JCLI-D-12-00417.1.
Arora
,
V.
, and Coauthors
,
2013
: Carbon–concentration and carbon–climate feedbacks in CMIP5 Earth system models. J. Climate,26, 5289–5314,
doi:10.1175/JCLI-D-12-00494.1
.
Boer
,
G. J.
, and
V.
Arora
,
2009
:
Temperature and concentration feedbacks in the carbon cycle
.
Geophys. Res. Lett.
,
36
, L02704, doi:10.1029/2008GL036220.
Bonan
,
G. B.
, and
S.
Levis
,
2010
:
Quantifying carbon-nitrogen feedbacks in the Community Land Model (CLM4)
.
Geophys. Res. Lett.
,
37
, L07401, doi:10.1029/2010GL042430.
Friedlingstein
,
P.
,
I.
Fung
,
E.
Holland
,
J.
John
,
G.
Brasseur
,
D.
Erickson
, and
D.
Schimel
,
1995
:
On the contribution of CO2 fertilization to the missing biospheric sink
.
Global Biogeochem. Cycles
,
9
,
541
556
,
doi:10.1029/95GB02381
.
Friedlingstein
,
P.
,
J.
Dufresne
,
P.
Cox
, and
P.
Rayner
,
2003
:
How positive is the feedback between climate change and the carbon cycle?
Tellus
,
55B
,
692
700
,
doi:10.1034/j.1600-0889.2003.01461.x
.
Friedlingstein
,
P.
, and Coauthors
,
2006
:
Climate–carbon cycle feedback analysis: Results from the C4MIP model intercomparison
.
J. Climate
,
19
,
3337
3353
, doi:10.1175/JCLI3800.1.
Garnier
,
E.
,
1992
:
Growth analysis of congeneric annual and perennial grass species
.
J. Ecol.
,
80
,
665
675
,
doi:10.2307/2260858
.
Gotangco Castillo
,
C. K.
,
S.
Levis
, and
P.
Thornton
,
2012
:
Evaluation of the new CNDV option of the Community Land Model: Effects of dynamic vegetation and interactive nitrogen on CLM4 means and variability
.
J. Climate
,
25
,
3702
3714
,
doi:10.1175/JCLI-D-11-00372.1
.
Gregory
,
J. M.
,
C. D.
Jones
,
P.
Cadule
, and
P.
Friedlingstein
,
2009
:
Quantifying carbon cycle feedbacks
.
J. Climate
,
22
,
5232
5250
,
doi:10.1175/2009JCLI2949.1
.
Hajima
,
T.
,
T.
Ise
,
K.
Tachiiri
,
E.
Kato
,
S.
Watanabe
, and
M.
Kawamiya
,
2012
:
Climate change, allowable emission, and Earth system response to representative concentration pathway scenarios
.
J. Meteor. Soc. Japan
,
90
,
417
433
,
doi:10.2151/jmsj.2012-305
.
Ito
,
A.
,
2011
: A historical meta-analysis of global terrestrial net primary productivity: Are estimates converging? Global Change Biol.,17, 3161–3175, doi:10.1111/j.1365-2486.2011.02450.x.
Jones
,
C.
, and Coauthors
,
2013
: Twenty-first-century compatible CO2 emissions and airborne fraction simulated by CMIP5 Earth system models under four representative concentration pathways. J. Climate, 26, 4398–4413,
doi:10.1175/JCLI-D-12-00554.1
.
Lewis
,
S. L.
, and Coauthors
,
2009
:
Increasing carbon storage in intact African tropical forests
.
Nature
,
457
,
1003
1006
,
doi:10.1038/nature07771
.
Meinshausen
,
M.
,
N.
Meinshausen
,
W.
Hare
,
S. C. B.
Raper
,
K.
Frieler
,
R.
Knutti
,
D. J.
Frame
, and
M. R.
Allen
,
2009
:
Greenhouse-gas emission targets for limiting global warming to 2°C
.
Nature
,
458
,
1158
1162
,
doi:10.1038/nature08017
.
Moss
,
R. H.
, and Coauthors
,
2010
:
The next generation of scenarios for climate change research and assessment
.
Nature
,
463
,
747
756
,
doi:10.1038/nature08823
.
Nakicenovic
,
N.
, and Coauthors
,
2000
:
IPCC Special Report on Emissions Scenarios.
Cambridge University Press, 599 pp.
Norby
,
R. J.
, and
D. R.
Zak
,
2011
:
Ecological lessons from free-air CO2 enrichment (FACE) experiments
.
Annu. Rev. Ecol. Evol. Syst.
,
42
,
181
203
,
doi:10.1146/annurev-ecolsys-102209-144647
.
Raupach
,
M. R.
,
J. G.
Canadell
,
P.
Ciais
,
P.
Friedlingstein
,
P. J.
Rayner
, and
C. M.
Trudinger
,
2011
:
The relationship between peak warming and cumulative CO2 emissions, and its use to quantify vulnerabilities in the carbon–climate–human system
.
Tellus
,
63B
,
145
164
,
doi:10.1111/j.1600-0889.2010.00521.x
.
Tachiiri
,
K.
,
A.
Ito
,
T.
Hajima
,
J. C.
Hargreaves
,
J. D.
Annan
, and
M.
Kawamiya
,
2012
:
Nonlinearity of land carbon sensitivities in climate change simulations
.
J. Meteor. Soc. Japan
,
90A
,
259
274
,
doi:10.2151/jmsj.2012-A13
.
Taylor
,
J. A.
, and
J.
Lloyd
,
1992
:
Sources and sinks of atmospheric CO2
.
Aust. J. Bot.
,
40
,
407
418
,
doi:10.1071/BT9920407
.
Taylor
,
K. E.
,
R. J.
Stouffer
, and
G. A.
Meehl
,
2012
:
An overview of CMIP5 and the experiment design
.
Bull. Amer. Meteor. Soc.
,
93
,
485
498
,
doi:10.1175/BAMS-D-11-00094.1
.
Todd-Brown
,
K. E. O.
,
J. T.
Randerson
,
W. M.
Post
,
F. M.
Hoffman
,
C.
Tarnocai
,
E. A. G.
Schuur
, and
S. D.
Allison
,
2013
:
Causes of variation in soil carbon simulations from CMIP5 Earth system models and comparison with observations
.
Biogeosciences
,
10
,
1717
1736
,
doi:10.5194/bg-10-1717-2013
.
Thompson
,
M. V.
,
J. T.
Randerson
,
C. M.
Malmstrom
, and
C. B.
Field
,
1996
:
Change in net primary production and heterotrophic respiration: How much is necessary to sustain the terrestrial carbon sink?
Global Biogeochem. Cycles
,
10
,
711
726
,
doi:10.1029/96GB01667
.
Thornton
,
P. E.
, and Coauthors
,
2009
:
Carbon-nitrogen interactions regulate climate-carbon cycle feedbacks: Results from an atmosphere-ocean general circulation model
.
Biogeosciences
,
6
,
2099
2120
,
doi:10.5194/bg-6-2099-2009
.
Tjiputra
,
J. F.
,
C.
Roelandt
,
M.
Bentsen
,
D. M.
Lawrence
,
T.
Lorentzen
,
J.
Schwinger
,
Ø.
Seland
, and
C.
Heinze
,
2013
:
Evaluation of the carbon cycle components in the Norwegian Earth System Model (NorESM)
.
Geosci. Model Dev.
,
6
,
301
325
,
doi:10.5194/gmd-6-301-2013
.
Watanabe
,
S.
, and Coauthors
,
2011
:
MIROC-ESM 2010: Model description and basic results of CMIP5-20c3m experiments
.
Geosci. Model Dev.
,
4
,
845
872
,
doi:10.5194/gmd-4-845-2011
.

Footnotes

*

Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/JCLI-D-13-00177.1.s1.

Supplemental Material