Hydrological and biogeochemical constraints on terrestrial carbon cycle feedbacks

2.1.1. CMIP5 models

2.1.2. Reference products for ET, NBP, precipitation and temperature

2.2.1. Carbon cycle feedbacks

We calculated the carbon cycle feedbacks in the CMIP5 models following the methodology introduced by Wenzel et al (2014), which differs little from that used in the Coupled Climate Carbon Cycle Model Intercomparison project (C⁴MIP; Friedlingstein et al (2006)). Based on this approach, the 'biogeochemically' coupled ('uncoupled') simulations were used to compute the concentration-carbon feedback, which can be translated into a model's CO₂ fertilization effect. In these simulations, the change in the land carbon storage (${\rm{\Delta }}\,{{\rm{C}}}^{{\rm{bgc}}-{\rm{coupled}}}$) is proportional to the change in atmospheric CO₂ concentration in these experiments (${\rm{\Delta }}\,{{\rm{CO}}}_{2}^{{\rm{bgc}}-{\rm{coupled}}}$) multiplied by the concentration-carbon feedback (β):

$$\begin{eqnarray}&&{\rm{\Delta }}\,{{\rm{C}}}^{{\rm{bgc}}-{\rm{coupled}}}=\beta \,\cdot \,{\rm{\Delta }}{{\rm{CO}}}_{2}^{{\rm{bgc}}-{\rm{coupled}}}.\end{eqnarray} \tag{ 1 }$$

The climate-carbon feedback was estimated from the fully coupled (1pctCO₂) experiments. Since the change in atmospheric CO₂ concentration is indistinguishable in the biogeochemically coupled and fully coupled experiments, the change in the land carbon storage in the fully coupled simulations (${\rm{\Delta }}{{\rm{C}}}^{{\rm{fullyc}}}$) equals the change in atmospheric CO₂ concentration in these simulations (${\rm{\Delta }}\,{{\rm{CO}}}_{2}^{{\rm{fullyc}}}$) multiplied by the concentration-carbon feedback (β) plus the temperature change (${\rm{\Delta }}\,{{T}}^{{\rm{fullyc}}}$) in the coupled experiment multiplied by the climate-carbon feedback (γ):

$$\begin{eqnarray}&&{\rm{\Delta }}{{\rm{C}}}^{{\rm{fullyc}}}=\beta \,\cdot \,{\rm{\Delta }}{{\rm{CO}}}_{2}^{{\rm{fullyc}}}+\gamma \,\cdot \,{\rm{\Delta }}\,{T}^{{\rm{fullyc}}}.\end{eqnarray} \tag{ 2 }$$

In all models except one, the two feedbacks were calculated during the period between the year 1 and year 140 of the idealized simulations. In the case of GFDL-ESM2M the feedbacks were calculated for a different time window, since the model hold the atmospheric CO₂ concentration after the year of doubling (year 70) constant. Therefore, we calculated the feedbacks for this model for the period between the year 28 (atmospheric CO₂ concentration corresponds to 2000 levels in the historical experiments) and the year 70 (atmospheric CO₂ concentration corresponds to 2045 levels based on the representative concentration pathway (RCP) 8.5 experiment; van Vuuren et al (2011)).

2.2.2. Response of NBP to climate variations

The response of NBP to interannual variations in temperature and precipitation was derived empirically based on a multiple regression approach as in Piao et al (2013). Based on this approach, NBP can be derived as a linear combination of anomalies in temperature and precipitation:

$$\begin{eqnarray}&&{\rm{NBP}}={\gamma }_{{\rm{reg}}}\,\cdot \,{x}_{T}+\delta \,\cdot \,{x}_{{\rm{PR}}}+\epsilon ,\end{eqnarray} \tag{ 3 }$$

where ${\gamma }_{{\rm{reg}}}$ and δ are the partial regression coefficients that represent the responses of NBP to temperature (PgC yr^–1 K^–1) and precipitation (PgC yr^–1 mm^–1), respectively. The symbols x_T and x_PR are the detrended anomalies in annual temperature and precipitation, respectively, and is the error term. We calculated the observation-based responses of NBP to temperature and precipitation (${\gamma }_{{\rm{reg}}}$ and δ) based on a total of 108 different combinations of the reference products (6 for NBP x 6 for precipitation x 3 for temperature) for the 1989–2005 period and using the 'regmultlin' function in NCL (https://ncl.ucar.edu/Document/Functions/Built-in/reg_multlin.shtml).

2.2.3. Observational constraints approach

We constrain the full ensemble of CMIP5 models (prior ensemble) based on the emergent relationships we establish in the result section (figure 1) and using several observation-based products as constraints. The constraints are applied globally based on the annual means for the period 1989–2005, which is determined by the availability of the observation-based products. The 'constrained' ensemble is constructed by removing models not lying within the uncertainty range of the observation-based constraints. Since we applied several constraints, we decided to exclude only models that are outside the range of all constraints (i.e., a model can be retained even if it complies with only one constraint).

Zoom In Zoom Out Reset image size

The climate-carbon feedback is constrained based on global contemporary variations in ET (ETvarcon) and NBP (NBPvarcon), as well as the response of NBP to variations in temperature (${\gamma }_{{\rm{reg}}}$con) and precipitation (δcon). A summary of the various constraints employed for γ is shown in table 2. A common characteristic in all these approaches is that the observed interannual variability of the terrestrial carbon fluxes can be used to constrain the long-term response of these fluxes, implying that the same processes operate at these different time scales. In the case of ET, the uncertainty range is defined based on the IQR of the different ET datasets contained in the LandFlux-Eval ET product. In the case of NBP and its responses to interannual variations in temperature (γ) and precipitation (δ), the uncertainty range is defined based on the full range of the various observation-based estimates.

Table 2.  Summary of the various constraint approaches used for the climate-carbon feedback (γ).

Ensemble	Rationale	Reference products	Definition of the uncertainty range
ETvarcon	ET is a proxy for GPP.	LandFlux-Eval synthesis	IQR of the different ET datasets
	Interannual variations in GPP are	dataset (Mueller et al 2013)	contained in the LandFlux-Eval
	a response to climate variations		synthesis product
	(figures S5, S6) and drive the variations
	in the carbon sink (figure 2(b), S3)
NBPvarcon	Interannual variations in NBP	NBP from the Global Carbon	Full range of the different reference
	(correlated with GPP variations;	Project and 5 atmospheric CO2	products for NBP
	figures 2(b), S3) are a response to	inversion products	(5 inversions and GCP estimate)
	climate variations
δcon and ${\gamma }_{{\rm{reg}}}$ con	The interannual variations in NBP	NBP from the Global Carbon	Full range of the 108 combinations
	are related to variations in climate	Project and 5 atmospheric CO2	between the various
	(e.g. temperature and precipitation;	inversion products, 6 precipitation	observation-based products
	figures S7, S8). The sensitivity of NBP	datasets and 3 temperature datasets
	to temperature and precipitation
	variations is related to γ

In the absence of direct observational constraints, the concentration-carbon feedback is constrained based on the strong inter-model relationship we identify between the two feedbacks in the CMIP5 models (figure 1(f)).

We find in the CMIP5 models a statistically significant inter-model correlation between the global climate-land carbon feedback strength (γ) and the magnitude of the interannual variations in global ET (r = −0.76, P = 0.029; figure 1(a)), NBP (r = −0.79, P = 0.019; figure 1(b)), GPP (r = −0.82, P = 0.013; figure 1(e)) and the response of NBP to variations in temperature (r = 0.94, P = 0.0005; figure 1(c)) and precipitation (r = −0.9, P = 0.0057; figure 1(d)) over the period 1989–2005. These different relationships can be used as a basis for constraining γ. While these correlations do not automatically imply a causal mechanism, we argue, in line with previous studies (Cox et al 2013, Wenzel et al 2014), that they emerge because the short-term (interannual) and the long-term (decadal to centennial) sensitivities of the terrestrial carbon fluxes to climate are highly linked.

More specifically, interannual variations in GPP, ET and NBP are essentially driven by variations in climate (figure S5,S6), as also shown in previous observational and modeling studies (Tian et al 1998, Luyssaert et al 2007, Reichstein et al 2007, Ciais et al 2009, Weber et al 2009, Friedlingstein et al 2010, Le Maire et al 2010, Jung et al 2011, Poulter et al 2014, Ahlström et al 2015, Anav et al 2015). The short-term (e.g. ET, NBP, GPP interannual variations) and the long-term (γ) sensitivities are correlated in the models (figures 1(a)–(e)) presumably because simulated GPP, ET and NBP respond to climate variations in a similar way whether these variations occur from year to year or over longer time scales. In other words, models that are particularly sensitive to interannual climate variations (i.e., larger anomalies in ET, GPP and NBP-also highly correlated each other, figures S1–S4) will also tend to have a higher sensitivity to future warming (higher γ).

GPP variations are strongly correlated with ET variations, in the CMIP5 models, both at the regional (figures 2(a), S2) and the global scale (r = 0.77, P = 0.025; figure S1(a)), reflecting the well-established coupling between plant photosynthesis and transpiration. As a consequence, the year-to-year variations in ET during the historical period are also highly correlated with γ. However, current observation-based products for global GPP do not offer robust estimates for the interannual variability (Anav et al 2015), preventing the emergent relationship between γ and the interannual variations in GPP to be used as a constraint at this time. But, the high relationship between GPP and ET variations offers the opportunity to use observation-based estimates of interannual ET variations as a proxy for interannual GPP variations and thus as a constraint on γ.

Zoom In Zoom Out Reset image size

In order to better understand the strong inter-model relationships between γ and the response of NBP to interannual variations in temperature (${\gamma }_{{\rm{reg}}}$) and precipitation(δ) (figures 1(c), (d)), we further investigate the link between ecosystems and climate variations, which are related to the climate-carbon feedback, by analyzing in both models and observations the relationship between contemporary interannual anomalies in NBP and in the potential drivers temperature (${\gamma }_{{\rm{reg}}}$) and precipitation (δ), respectively. The analysis focuses on temperature and precipitation since these two climate variables are known to be key determinants of the productivity of ecosystems (Nemani et al 2003, Poulter et al 2014, Ahlström et al 2015, Anav et al 2015, Murray-Tortarolo et al 2016, Seddon et al 2016).

During the period 1989–2005, the different combinations between the observation-based products for NBP, temperature and precipitation agree on a negative response of NBP to variations in temperature (${\gamma }_{{\rm{reg}}}$) ranging from about −3.8 to −2.5 PgC yr^–1 K^–1. In the case of precipiation (δ), the results range between −0.8 and 3.4 PgC yr^–1 mm^–1 with most of the combinations showing a positive response of NBP to interannual variations in precipitation. The negative relationship between NBP and temperature translates into a lower terrestrial carbon sink during warm years, caused either by increased respiration in response to warmer temperatures (Kätterer et al 1998, Davidson and Janssens 2006, Bond-Lamberty and Thomson 2010, Conant et al 2011, Lu et al 2013) and/or by lower photosynthesis in regions where plants operate close to their optimal temperature. The positive relationship between NBP and precipitation suggests that wet years are favorable for plant productivity, resulting in enhanced terrestrial carbon uptake. This mechanism is valid in several ecosystems in the mid-latitudes and in the tropics. In particular, in tropical ecosystems the relationship between NBP and precipitation (as well as temperature) anomalies is more pronounced during ENSO events (Jones and Cox 2001, Gurney et al 2003). Warm and dry conditions during these events cause substantial reductions in plant production and enhance respiration resulting in lower storage on land in tropical land ecosystems (Gatti et al 2014) and in a stronger climate-carbon feedback that relates to interannual variations in climate. Eventually, these anomalies in the carbon sink are mirrored in the atmospheric CO₂ growth rate.

Consistent with the observation-based products, the CMIP5 models show negative and positive responses of NBP to interannual variations in temperature and precipitation, respectively, both globally (figures 1(c), (d), S9, S10) and regionally (figure S7,S8). However, they simulate a large spread in the strength of these relationships, with most of the models tending to be oversensitive to interannual variations in precipitation (figure 1(d)) compared to the reference products (as already shown in Piao et al (2013) for the global averages in the offline TRENDY simulations). The relationship between γ and the responses of NBP to variations in temperature and precipitation are strong and statistically significant also for a longer time period in the observation record where fewer reference data products are available (e.g., 1985–2009, see figures S1(c), (d)). We note here that the high soil moisture-dependence of soil respiration in MPI-ESM-LR (Wenzel et al 2014) decouples the relationship between NBP and precipitation (negative) in this model mainly owing to enhanced decomposition in wet years, which overwhelms the possible enhancement of plant productivity in response to higher soil water availability. For this reason, we have excluded MPI-ESM-LR when calculating the best fit linear regression between γ and δ (figure 1(d)).

These 'responses' of NBP to climate variations need to be interpreted carefully since the multi-linear regression approach does not take into account the effect of synergistic/antagonistic effects between the two drivers and also neglects other confounding mechanisms and drivers of plant productivity. Thus, this method for calculating ${\gamma }_{{\rm{reg}}}$ and δ might describe reasonably well the response of NBP to temperature and precipitation in the models, since most of them lack many other processes that influence the carbon cycle, while it might not always work for the observation-based products. Indeed, this seems to be the case in the models since there is a strong and significant correlation between γ and the response of NBP to interannual variations in temperature (${\gamma }_{{\rm{reg}}};$ figure 1(c)). Thus, the assumption that the same processes control γ in the short-term and long-term scales is valid at least in the model world.

The magnitude of the climate-carbon feedback varies markedly across the CMIP5 models (table 3), as already shown earlier (Arora et al 2013, Ciais et al 2013). Several reasons might explain this large range including, for instance, biases in the underlying simulated climate or the diverse and sometimes inadequate representation of land biogeochemical and biophysical processes in the different ESMs (Plattner et al 2008, Hawkins and Sutton 2009, Anav et al 2013, Piao et al 2013, Hoffman et al 2014, Friedlingstein 2015). All models simulate a negative γ ranging from about −3.7 (NorESM1-ME) to −89PgC/K (GFDL-ESM2M) with a multi-model mean of −45.9 and a model range of 85 PgC K^–1, respectively. The low sensitivity in CESM1-BGC (also seen in NorESM1-ME) was attributed earlier to increased nitrogen mineralization, caused by increased respiration that provides additional nitrogen inputs to vegetation and compensates as a result, the carbon losses through the process of respiration (Arora et al 2013)—an effect not considered in the other carbon-only ESMs. The simulated magnitude of the climate-carbon feedback in the CMIP5 models lies within the range of the C⁴MIP models (Friedlingstein et al 2006). The different observation-based constraints suggest a substantially lower response of the terrestrial carbon cycle to climate change (γ) and a significant decrease in the inter-model range. Specifically, the magnitude of the climate-carbon feedback for the constrained multi-model mean is lowered by 40% (figure 3), from −45.9 to −27.4 PgC K^–1, while the observational constraints half the inter-model spread from 85 to 40 PgC K^–1. Another indication that γ might be overestimated in the CMIP5 models relates to the tendency of some models to simulate high total carbon in vegetation and soils (Todd-Brown et al 2012, Anav et al 2013, Carvalhais et al 2014) and thus high amounts of carbon, contained in these pools, are exposed to the warming.

Zoom In Zoom Out Reset image size

We find a statistically highly significant negative correlation between γ and β in the ESMs (r = −0.92, P = 0.0001; figure 1(f)) with models having a more negative γ having a more positive β. This strong relationship is likely the consequence of a pool effect. In other words, models with a higher β, i.e., stronger fertilization effect, tend to accumulate more carbon in vegetation and soils, thereby increasing the risk of carbon losses due to climate change. Keeping in mind that this relationship is a model-derived feature not directly verifiable based on observations, we examined the consequences of the constrained γ on β. All CMIP5 models simulate a positive β, but with substantial inter-model variation, ranging from about 0.24 PgC ppm^–1 for CESM1-BGC and NorESM1-ME to about 1.76 PgC ppm^–1 for GFDL-ESM2M, resulting in a multi-model mean of 0.91 PgC ppm^–1. The much lowest response to CO₂ in CESM1-BGC and NorESM1-ME originates likely from the strong limitation that nitrogen imposes on the strength of CO₂ fertilization compared to carbon-only models (Thornton et al 2009, Zaehle et al 2010). In the case of the MIROC-ESM model, the relatively low response to CO₂ (0.49 PgC ppm^–1) can be attributed to the empirical treatment of photosynthesis in its land surface module, implicitly including the role of nutrient limitations on CO₂ fertilization (Ito and Oikawa 2002, Arora et al 2013). Selecting only the models complying with the observation-based constraints described in the previous section, results in a substantially lower estimate of γ. Specifically, the constraint reduces the concentration-carbon feedback by 30% from 0.9 PgC ppm^–1 to 0.6 PgC ppm^–1 (figure 3) and also reduces the range by 40% (0.2–1.1 PgC ppm^–1).

By the end of the fully coupled experiments corresponding to a quadrupling of the preindustrial atmospheric CO₂ concentration, all models simulate a high cumulative carbon sink in the terrestrial biosphere (table 3, figure S11(b)), ranging from 152PgC (MIROC-ESM) to 818 PgC (MPI-ESM-LR). The change in the annual net land-atmosphere flux, calculated as the difference in NBP between the last and the first 20 years of the simulation, also shows a large range, with all but two models (GFDL-ESM2M and MIROC-ESM) simulating an increased sink by the end of the simulation (figure S11(d)). In those models that simulate an increase in NBP, this change is driven by the increase in GPP, that more than compensates the increased losses stemming from the faster turnover time of total carbon in soils and vegetations (table 3). The latter is mainly driven by faster soil decomposition and to a lesser extent by faster autotrophic respiration rates in response to the warming (table 3). The decrease in NBP in GFDL-ESM2M and MIROC-ESM is due to the faster turnover time of soil carbon, that more than compensates the increase in GPP. For GFDL-ESM2M, this is a consequence of the model's setup since, as mentioned earlier, the forcing in this model stabilized at the year of CO₂ doubling. Hence, GPP in GFDL-ESM2M stays constant in these experiments after this moment, allowing respiration to regulate NBP and reduce the sink strength in this model.

Obviously, the constrained feedbacks have an impact on the future change in the land carbon storage, since their balance determines NBP (see equation (2)). The observation-based constraints suggest a multi-model mean NBP that is about 19% lower (figure 3) and they reduce the inter-model range of the change in the land carbon storage by about 15% to a range of 151–718 PgC. In the constrained models, the lower CO₂ fertilization effect still more than compensates the lower response to climate change, resulting in a lower terrestrial carbon uptake. These results concerning the magnitude of the future changes in NBP are not easy to extrapolate to other scenarios, since the numbers are likely quite sensitive to the specific emission/concentration scenario applied. However, the lower sink in the constrained ensembles can be related to results from other scenarios (e.g., RCP8.5 experiment) since models simulating large changes in NBP in the idealized 1%/yr increase in CO₂ experiments simulate also large changes in future NBP in the RCP8.5 simulations (table 3, r = 0.96, P = 0.0012). By the end of the RCP8.5 experiment, the models simulate, on average, a cumulative change in NBP of 102 PgC (table 3) while for the same experiment the constrained models have a storage on land that is about 50% lower (53 PgC).