Greenhouse Gas Concentration and Volcanic Eruptions Controlled the Variability of Terrestrial Carbon Uptake Over the Last Millennium.

The terrestrial net biome production (NBP) is considered as one of the major drivers of interannual variation in atmospheric CO2 levels. However, the determinants of variability in NBP under the background climate (i.e., preindustrial conditions) remain poorly understood, especially on decadal-to-centennial timescales. We analyzed 1,000-year simulations spanning 850-1,849 from the Community Earth System Model (CESM) and found that the variability in NBP and heterotrophic respiration (RH) were largely driven by fluctuations in the net primary production (NPP) and carbon turnover rates in response to climate variability. On interannual to multidecadal timescales, variability in NBP was dominated by variation in NPP, while variability in RH was driven by variation in turnover rates. However, on centennial timescales (100-1,000 years), the RH variability became more tightly coupled to that of NPP. The NBP variability on centennial timescales was low, due to the near cancellation of NPP and NPP-driven RH changes arising from climate internal variability and external forcings: preindustrial greenhouse gases, volcanic eruptions, land use changes, orbital change, and solar activity. Factorial experiments showed that globally on centennial timescales, the forcing of changes in greenhouse gas concentrations were the largest contributor (51%) to variations in both NPP and RH, followed by volcanic eruptions impacting NPP (25%) and RH (31%). Our analysis of the carbon-cycle suggests that geoengineering solutions by injection of stratospheric aerosols might be ineffective on longer timescales. ©2019. The Authors.

Introduction

The climate system has experienced a global warming of about 0.85°C since 1850, with changes in radiative forcings due to increased long‐lived atmospheric greenhouse gases and changes in aerosols and other forcings (Intergovernmental Panel on Climate Change, 2013). The degree of climate warming relates to anthropogenic emissions as well as feedbacks, for example, of the carbon‐cycle‐climate system (Heimann & Reichstein, 2008). For example, the global terrestrial ecosystem acted as a negative carbon‐climate feedback that effectively slowed down global warming by continually sequestering about 30% of anthropogenic CO2 emissions during the past 50 years (e.g., Le Quéré et al., 2009, 2014; Schimel et al., 2015). Currently, the global terrestrial carbon uptake (hereafter referred to as the net biome production or NBP) is a net sink for atmospheric CO2; the NBP is affected by many factors related to human activity or climate change, including CO2 levels, surface weather (also termed “climate drivers”), land use change, and nutrient availability.

The estimated global NBP since the 1960s (from modeling and carbon inventories) shows significant influences of both natural variability and anthropogenic forcings, as can be seen in the long‐term trend and the residual variability (e.g., Piao et al., 2013; Schimel et al., 2015; Shevliakova et al., 2013; Zhang et al., 2018). One complicating factor is that the superposition of natural and anthropogenic influences on the terrestrial carbon cycle makes it difficult to disentangle those two influences, and the more readily available instrumental observations and model simulations over the last 100 years are too short to draw strong conclusions on multidecadal and longer‐term variability (Lehner et al., 2015). Therefore, there is a need to study the causes of variability in global NBP under the climate state without the strong anthropogenic forcings (i.e., the background climate state). During the 1,000 years ending in 1850, the variability of NBP was governed by internal variability of carbon uptake and losses and their responses to natural external forcings (e.g., volcanic eruptions and solar activity) with a limited contribution of early land‐use change to CO2 and CH4 (Ferretti et al., 2005; MacFarling Meure et al., 2006). The preindustrial last millennium (PILM, 850‐1849) thus constitutes a test‐bed to characterize multiannual timescales of variability in global NBP and its driving mechanisms without the complications of anthropogenic influences. This allows us to study the timescale‐dependent response of the land carbon cycle to natural climate variability and external natural forcings, and how they together influenced the long‐term variability of atmospheric CO2 levels in the past 1,000 years.

The NBP depends on the imbalance between the carbon input to ecosystem pools from net primary production (NPP) and output from heterotrophic respiration (RH) and disturbances (D) such as wild fires (ignoring losses from fluvial export etc.); all of these fluxes vary on different timescales and across land regions. They are also coupled via the terrestrial carbon pools. In general, terrestrial carbon cycle processes respond either rapidly (e.g., photosynthesis) or slowly (e.g., soil carbon accumulation) to climate (Gregory et al., 2009; Huntzinger et al., 2017; Jung et al., 2017; Luo & Weng, 2011; Zhang et al., 2016). The ecosystem input and output fluxes are coupled through internal carbon processes such as turnover and allocation, so that NPP, RH, and D will respond differently on interannual to centennial timescales to atmospheric CO2 concentrations and climate (solar radiation, humidity, temperature, and precipitation, etc.). Interannual timescales of about 5‐7 years represent a prominent mode of NBP variability, associated with key climate modes such as the El Niño‐Southern Oscillation (ENSO) and the North Atlantic Oscillation, which influence regional climate, therefore carbon cycle (Bastos et al., 2016; Fang et al., 2017; Wang et al., 2013; Zeng et al., 2005; Zhang et al., 2018; Zhu et al., 2017). Decadal to multidecadal variability in NBP is associated with other climate‐ocean modes, namely, the Pacific Decadal Oscillation and the Atlantic Multi‐decadal Oscillation (Zhang et al., 2018; Zhu et al., 2017). Furthermore, sporadic volcanic eruptions may also significantly influence terrestrial land carbon uptake (e.g., Gu et al., 2003; Mercado et al., 2009). How terrestrial carbon cycle and its internal processes responded to those different natural external forcings and different time scales remains largely unresolved.

Key external forcings on the climate system include volcanic aerosols (Gu et al., 2003; Mercado et al., 2009; Miller et al., 2012), solar variation and change in orbital parameters (Kaufman et al., 2009), and land use and land cover change (LULCC) (Devaraju et al., 2015; Erb et al., 2016; He et al., 2014; Kloster et al., 2012; Wang et al., 2015); the influence of some of these forcings on NBP may be difficult to study as they occur either sporadically (e.g., volcanoes) or on timescales of decades or longer (e.g., land‐use change). An ecosystem's carbon cycle may also be limited by nutrient supply such as soil nitrogen and phosphorus content (Fleischer et al., 2013; Goll et al., 2012; Greaver et al., 2016; Wang et al., 2015; Zaehle & Friend, 2010; Zhang et al., 2011, 2013). In order to isolate the response of the carbon cycle to a particular external forcing, one requires a series of model simulations, controlling one forcing at a time (Arora et al., 2013; Gerber et al., 2003).

The magnitude of NBP and its response to forcing depends on the ecosystem turnover rate (hereafter defined as k eco) of carbon pools in the ecosystems, and on the sensitivity of D and NPP to climate, the latter depending among other things on the biophysical coupling between soil moisture and photosynthesis (Fung et al., 2005; Luo et al., 2017; Luo & Weng, 2011; Taylor & Lloyd, 1992). The turnover time (the reciprocal of k eco) represents the average time a carbon atom spent in the ecosystem, from input through NPP until exit via D or RH. Hence, the turnover time controls the efficiency of an ecosystem to sequester carbon for a given NPP. Averaged over large spatial scales, the turnover time can vary from hundred years (e.g., in boreal forests) and to a few years (e.g., in grasslands, some tropical rainforests, or semiarid regions; Bloom et al., 2016; Carvalhais et al., 2014; Lu et al., 2018). Values of turnover rate (k eco) can also vary in time in response to climate variability over decades to longer periods, which may potentially impact NBP and cause a longer‐term adjustment to the land carbon cycle. To analyze the influences of those external factors on terrestrial carbon uptake on different time scales, here we used the theoretical framework as developed by Zhang et al. (2018) to estimate the contributions from NPP input and k eco to NBP at different timescales. During the PILM period, the carbon‐climate feedback system experienced natural climate internal variability and external forcings; focusing on this period allows us to quantify better the influences of some decadal variations of some climate modes and volcanic eruptions on the terrestrial carbon cycle (e.g., the role of NPP and k eco fluctuations on carbon storage).

In this study, we analyzed an ensemble of simulations from a coupled carbon‐climate model (earth system model [ESM]), the Community Earth System Model Last Millennium Ensemble (CESM‐LME) over the period 850‐2005, considering all external forcings together or each of them separately (Otto‐Bliesner et al., 2016). The goal is to gain a better understanding of the response of the carbon cycle in this model to internal climate variability (e.g., temperature, precipitation, solar radiation, and climate modes) and the different external forcings. The CESM‐LME ensemble allows a longer‐term perspective for detection and attribution of the effect of anthropogenic or natural external forcings (Otto‐Bliesner et al., 2016). Using the theoretical framework developed by Zhang et al. (2018), we attributed variability in NBP and RH on different timescales to NPP and k eco as simulated by the CESM‐LME. Additionally, we compared contributions of the individual external forcings (e.g., greenhouse gases and volcanic eruptions) to terrestrial carbon fluxes on different timescales. Finally, we examined the details of the land carbon cycle between the Medieval Climate Anomaly (950‐1250) and the Little Ice Age (1450‐1849).

Materials and Methods

The CESM‐LME Ensemble Simulations

The simulations analyzed come from the CESM‐LME project (Otto‐Bliesner et al., 2016), made available via the Earth System Grid (http://www.earthsystemgrid.org). The CESM version 1.1 (Hurrell et al., 2013) was integrated with the land carbon‐atmosphere system coupled at a spatial resolution of roughly 2°×2°, and the ocean and sea ice system at approximately 1°×1° (the CESM grid is encoded as 1.9x2.5_gx1v6). The land component of the ESM is the Community Land Model version 4 coupled with carbon‐nitrogen cycle (CLM4CN) including key biophysical, biogeochemical processes, and human‐made disturbances on terrestrial ecosystems (Hurrell et al., 2013; Lawrence et al., 2011; Oleson et al., 2013; Thornton et al., 2007). The CESM‐LME simulations were shown to reproduce interannual variability and decadal mean values of land and ocean carbon sinks comparable to other ESMs of the Coupled Model Inter‐Comparison Project Phase 5 (CMIP5) when compared to observational estimates during the late twentieth century (Hoffman et al., 2014; Hurrell et al., 2013; Lehner et al., 2015). The CMIP5 forcing (Schmidt et al., 2011) prescribed to the CESM includes solar irradiance, changes from Earth's orbital parameters, solar variability with an 11 years periodicity, volcanic eruptions, LULCCs, and atmospheric concentrations of long‐lived and well‐mixed greenhouse gases (GHGs; namely, CO2, CH4 and N2O). The temporal evolution of these major external forcings has been shown in Fig. 2 in Otto‐Bliesner et al. (2016). The variation in this record of GHGs (derived from ice cores in Antarctica) relates not only to natural feedback in the carbon and nitrogen cycles to changes in climate but also to some early anthropogenic perturbations, such as from land use change due to agriculture (Gasser & Ciais, 2013; Gerber et al., 2003; Schmidt et al., 2011). The imposed LULCC is the only forcing applied to land areas; this was derived from a historical mapping of agricultural areas and population‐based proxies for agricultural activity. The resulting LULCC affects the surface albedo, surface roughness, soil characteristics, water cycling, and carbon turnover. The CO2 emissions due to land use change are not considered in LULCC, since they are implicitly incorporated in the reconstructed GHG record. We noticed that there is a small step‐change in the LULCC data at the year 1500, which is an artifact from the merger of two different land‐cover datasets (Otto‐Bliesner et al., 2016), scaling the Pongratz reconstruction (Pongratz et al., 2008) to match the Hurtt dataset (Hurtt et al., 2011) at the year 1500 at each land grid cell. Our calculation shows that this artificial change at 1500 of the LULCC forcing has led to a large emission of the CESM‐LME modeled D and significant influence on NBP at that year (Figure S2 in the supporting information). To avoid its potential effect to our results, we replaced the values at the year 1500 of all variables (NBP, NPP, RH, D, and all carbon pools) from all simulations with an average of 5 years (1495‐1499) before and 5 years (1501‐1505) after 1500.

The CESM‐LME ensemble provides 13 simulations with all forcings (hereafter termed the “full‐forcing” run) and smaller ensembles with each transient forcing separately (hereafter referred to as the “single‐forcing” runs) with five simulations for volcanic eruptions, four for the solar forcing, and three for the other forcings from 850 to 2005. All ensemble simulations used the model state in the final year of an “850 control simulation” as the initial conditions; this control simulation ran from year 651 to 2005 with all external forcings being kept at their levels for the year 850 over the entire simulation period (Figure S1). In the single‐forcing runs, all other forcings were also kept at their levels for the year 850 over the entire period. Ensemble spreads were generated using small random round‐off differences in the initial atmospheric state (Otto‐Bliesner et al., 2016). Additional information about the forcing datasets, model variables, and diagnostics can be found at the CESM‐LME website (http://www.cesm.ucar.edu/projects/community‐projects/LME/). The ensemble outputs used here were the annual terrestrial NBP, NPP, RH, D, all carbon pools (in the vegetation and soil), and land surface climate variables. Ensemble means for all the full‐forcing or single‐forcing simulations were calculated in order to reduce uncertainties arising from climate internal variability at long‐term timescales, when focusing on the effects of external forcings.

Theory to Attribute Variability in NBP to NPP and k eco

It can be shown (Zhang et al., 2018; see Text S1 in the supporting information) that for a fixed timescale, the amplitude of temporal variability in NBP, NPP, and the sum RH and D can be related as follows: where A NPP, A RH+D, and A NBP denote the amplitude of NPP, RH plus D, and NBP, respectively, at a given frequency (units of year‐1, or angular frequency ω). The turnover rate (k eco) is defined as the ratio of the output (RH + D) to the total C pool within an ecosystem (C eco):

Equations (1) and (2) show that on a given timescale (i.e., for fixed ω), the magnitude of variability in RH and NBP is driven solely by changes in NPP input and variations in k eco. The turnover time (τ=1/k eco) can vary from years to centuries due to internal climate variability and changes of external forcing. In this study, we applied this theory to investigate the temporal variation in k eco and how it affects NBP variability on different timescales.

To isolate the relative contributions of NPP and k eco to the variability in NBP and RH at different timescales, we built a three‐C‐pool box model (henceforth termed the “Box‐Model”) that reproduces the response of the complex ESM given levels of NPP and D (Zhang et al., 2018). The Box‐Model simplifies the C pool structure of the CESM/CLM4CN, which includes C pools of leaf, grain, live and dead stems, fine root, live and dead coarse roots, three litter pools, coarse woody debris, and four soil organic matter pools; all of these are aggregated into three pools: a plant biomass pool, a litter pool, and a soil organic matter pool (see Figure S2). The Box‐Model is described as where is a vector of carbon pool sizes; NPP and D are the annual output from each of the CESM‐LME full‐forcing simulations with year t from 850 to 1849. A is a transfer matrix, representing the fractions of carbon transferred from one pool to the others. Following Zhang et al. (2018), we have

The Box‐Model's parameters (a 21, a 22, a 32, and a 33) were calibrated separately for each terrestrial grid‐cell with the output NBP, NPP, RH, D, and all C pools from the CESM‐LME over the PILM. Once calibrated, the Box‐Model simulates annual RH and NBP given a time series of NPP, D, and the fixed turnover rates. As the D shows relative small contribution to the variance of NBP on decadal‐to‐centennial timescales as shown in Figures 1d and 2d. In this study we focused on how NPP and turnover rate impact NBP and RH. The NBP and the RH from the ESM output can then be analyzed into components driven by NPP versus k eco via the framework outlined in equations (1) and (2) combined with the Box‐Model simulations. The NBP and the RH from the Box‐Model define the degree to which these fluxes are driven by variation in NPP only. The NBP and the RH due to variation in k eco only are estimated as NBPCESM – NBPBox‐Model and RHCESM – RHBox‐Model, respectively, given that k eco is held fixed in the Box‐Model.

Figure 1

Time‐series of the ensemble mean of the global annual net biome production (NBP), net primary production (NPP), heterotrophic respiration (RH), and disturbance (D) fluxes, based on the 13 Community Earth System Model Last Millennium Ensemble (CESM‐LME) full‐forcing simulations over the period 850‐2005. The industrial period (1850 onwards) is highlighted in light grey. The shaded‐areas are the mean ± σ, which represent for uncertainties; σ was calculated as the ensemble standard deviation for each variable.

Figure 2

Power spectra for the time‐series global annual net biome production (NBP), net primary production (NPP), heterotrophic respiration (RH), and disturbance (D) fluxes from the 13 Community Earth System Model Last Millennium Ensemble (CESM‐LME) full‐forcing simulations over the period 850‐1849. The shaded areas span the spectra of individual ensemble members and the ensemble mean of the individual spectra is overlaid.

Power Spectral Analysis for Temporal Variability in Carbon Fluxes

From equations (1) and (2), it can be seen that for a given timescale (or a time frequency; Zhang et al., 2018),

The variability of NBP, NPP, and RH plus D on a given timescale can be described as the variances of these carbon fluxes, given as the square of these amplitudes (A, A, and A). As there is relatively small variation in global annual D, we assumed that . Equation (8) then was used to attribute the contributions of NPP and RH to variability in NBP. To estimate the amplitude of variability at different timescales (i.e., to estimate A NPP), we analyzed the power spectrum of the time series of annual NBP, NPP, RH, and D from both the CESM‐LME and the Box‐Model simulations; this allowed us to quantify the variability of these fluxes on timescales of 2‐1,000 years.

Results

Figure 1 shows the time series of global annual NBP, NPP, RH, and D over the period 850‐2005, displayed as the ensemble mean of the 13 CESM‐LME full‐forcing simulations. It was found that the operation of averaging over the individual ensemble members largely dampens the magnitude of fluctuations in the C fluxes; thus, the “mean spectrum of the ensemble” refers to the mean of the individual spectra rather than the spectrum of the ensemble mean that was driven by climate internal variability (Figures 2 and 3). This is because the timing and magnitude of climate‐internal variability (e.g., El Niño) in the CESM‐LME are uncorrelated in each simulation with different initial conditions (the detrended time‐series showed correlations of R 2 < 0.1). Therefore, the ensemble mean can help to identify main signals of externally forced variability and trends (Figure 1).

Figure 3

Spectra for the global annual (a, c, and e) net biome production (NBP) and (b d, and f) heterotrophic respiration (RH) over the period 850‐1849 simulated by Community Earth System Model Last Millennium Ensemble (CESM‐LME; top row), the Box‐Model (middle row), and the differences between the two models (bottom row). The NBP (or RH) simulated by the CESM‐LME full‐forcing runs was driven by both NPP and varying turnover rate (k eco), while the NBP (or RH) simulated by the Box‐Model (using invariant k eco as parameters) was driven by NPP input only. The NBP (or RH) difference between CESM‐LME and the Box‐Model was interpreted as the contribution to variation in NBP (or RH) due to variability in k eco.

In the full‐forcing run, we found no significant trend in modeled global NBP for the PILM (850‐1849), but large anomalies (see section 3.3 for more analysis) were associated with major volcanic eruptions (e.g., Samalas in Indonesia in 1258, Kuwae of Vanuatu in 1452, and Tambora in Indonesia in 1815; Figures 1a–1c). Over the industrial period (1850‐2005), the modeled NBP, NPP, and RH show significant increasing trends, accelerating after the 1950s (Figure 1).

Variation in NPP, RH, and NBP at Different Timescales

In the context of discussing variability at different temporal scales, we adopt the following terms for specific timescales: interannual (with a period of 2‐10 years), decadal (10‐100 years), and centennial (100‐1,000 years). To quantify the contributions of NPP, RH, and D to the variability of NBP at different timescales, we calculated the power spectra of the time‐series of global NBP for each of the 13 CESM‐LME full‐forcing simulations over the PILM. The mean of the ensemble of individual power spectra for NBP displays two peaks with periods of around 4 and 50 years, which fall into the interannual and multidecadal timescales, respectively (Figure 2a). The spectrum of the NBP time series from the CESM‐LME members displays the strongest amplitude of variation on interannual timescales, namely, at the period of 4 years with a peak value of about 16 (PgC)2/year‐2 and the period of 50 years with a peak value of about 3 (PgC)2/year‐2. As seen in the power spectrum for NPP (Figure 2b), there is a similar magnitude and location (corresponding to ENSO‐like time‐scales). By contrast, the variances of RH and D are smaller than that of either NPP or NBP (Figure 2), as the sum of variance of RH and D over timescales of 2‐10 year is less than 10% of that of NBP. The peak in the spectrum for D also corresponds to an ENSO‐like timescale (Figures 2c and 2d). This may be driven by variation in biomass burning in semiarid ecosystems associated with major climate modes (O'Donnell et al., 2011).

On centennial timescales (100‐1,000 years), the spectrum of global NBP does not show obvious fluctuations, with a variance of only around 0.5 (Pg C)2/year‐2 across periods from 100 to 1,000 years. In contrast, the spectra for NPP and RH show significant increasing variances from approximately 1 to 5 (Pg C)2/year‐2 with periods increasing from 100 to 1,000 years (Figures 2b and 2c); similar variation is not seen for the NBP, as a result of near cancelation between the variations in NPP and RH (e.g., similar magnitudes and opposite contributions to NBP) on centennial timescales. Regional analysis shows that the variability in the CESM‐LME modeled global NBP is attributable almost exclusively to fluxes in the tropical regions (20°S‐20°N) on interannual to multidecadal timescales, but on centennial timescales is jointly controlled by the tropics (55%) and the area north of 20°N (45%; Figs. S4 and S5).

Attributing Variability in NBP and RH to NPP and Turnover Rate

To attribute how variations in NBP and RH are respectively driven by internal dynamics (k eco) or variance in carbon input (NPP), we calculated the power spectrum of the output of the Box‐Model (see Methods). The global NBP driven by NPP with the fitted constant k eco only (i.e., from the Box‐Model) nearly reproduces the power spectrum of the CESM‐LME full‐forcing NBP spectrum (driven by NPP and k eco). This implies that variability in NPP is the key determinant for variation in NBP from interannual to multidecadal timescales over the PILM (Figures 3a and 3c). The role of variability in k eco to variation in NBP contributes a small amount of (~13%) variance over timescales of 2‐100 years. However, over centennial timescales the contribution variability in k eco to variance in NBP increases to 32% (Figure 3e).

At regional scales, across tropical rainforests and southern hemispheric semiarid regions, NPP‐driven variability in NBP accounts for more than 80% of the variance of NBP over all timescales from inter‐annual to centennial (Figures 4a, 4c, and 4e). But for boreal forests and alpine grasslands, the variance in NBP is split nearly evenly between variability associated with NPP (50‐60%) and k eco (40‐50%) on interannual timescales (Figures 4b, 4d, and 4df).

Figure 4

Spatial patterns of the relative contributions (expressed as a fraction) of the variance in (a, c, and e) net biome production (NBP) driven solely by variation in net primary production (NPP) and (b, d, and f) NBP driven by varying turnover rate (k eco) only on timescales of 2‐10 years (a and b), 10‐100 years (c and d), and 100‐1000 years (e and f), respectively. The variance in the time‐series of NBP for a particular timescale (e.g., 2‐10 years) is calculated by integrating the power spectrum of NBP across the range of frequencies. The NPP estimates come from the Community Earth System Model Last Millennium Ensemble (CESM‐LME) full‐forcing simulations over 850‐1849, the NBP driven by NPP only was simulated by the Box‐Model, and the NBP driven by k eco was taken to be the difference between estimates of NBP from the CESM‐LME and the Box‐Model.

Variation in modeled RH is dominated by variability in the turnover rate, k eco, on timescales of 2‐100 years (97%), and by NPP on timescales of 100‐1,000 years (76%; Figures 3b, 3d, and 3f). Spatial patterns show that the variations of regional RH in almost all ecosystems are mainly (>80%) contributed by k eco on interannual and decadal timescales (Figures 5a‐5d). However, on centennial timescales, the variance in RH is dominated by NPP (>70%), except in boreal forests where variation of k eco in response to external forcings accounts for 40‐60% (Figures 5e and 5f). The NPP‐driven RH variations on centennial timescales over the southern hemispheric semiarid regions contribute more than 90%, which may be due to the relatively thin soil C stock in these water‐limited ecosystems. While, the larger fraction of k eco‐driven RH resulted from the response of the deep and thick soil C storage over the boreal forests to long‐term climate variability and external forcings.

Figure 5

The same as Figure 4 but for (a, c, and e) variance in heterotrophic respiration (RH) driven by net primary production (NPP) only or by (b, d, and f) variation in the turnover rate (k eco) only (right column).

Contributions of Different External Forcings to Global NBP Across Different Timescales

It should be noted that the NBP is not only driven by the main components of terrestrial ecosystem (i.e., NPP, RH+D, or k eco) but also affected by forcings external to the terrestrial ecosystem in the ESM, such as temperature; precipitation; radiation, which impact upon plant photosynthesis and respiration (i.e., NPP); and soil respiration (i.e., RH), and these effects may also vary across timescales. We investigate in this section the contributions of key external forcings (GHGs, LULCC, orbital, solar, and volcanic) to variability in NBP on different timescales over the PILM (Figure 6). Compared to the full‐forcing and the 850‐year control runs, the power spectra of the single‐forcing simulations of global NBP do not show significant differences across all timescales (Figure 6a). However, the variances at low frequency (timescales of 100 to 1,000 years) of global NPP and RH display clearer differences between the single‐forcing and full forcing simulations over the PILM of 850‐1849. These differences generally increase with decreasing frequency (i.e., increasing in timescale), varying from roughly 0.5 to 5.0 (Pg C)2/year‐2 across timescales as compared to the 850‐year control simulation (Figures 6b and 6c). Variance associated with variation in GHG concentrations dominates the differences (51%) on the timescales of 100 to 1,000 years, followed by volcanic eruptions (25% for NPP and 31% for RH) and LULCC (9% for NPP and 5% for RH). The other forcings (orbital, solar, and volcanic) have much smaller contributions (15% for NPP and 13% for RH) (Figures 6b–6d). Unsurprisingly, the LULCC forcing is the largest contributor to centennial‐scale variability in D (Figure 6d). Overall variation in GHG concentrations has the strongest impact upon variability in global NBP across timescales of >200 years (Figure 6b and 6c), while volcanic eruptions over the PILM are the most important external factors on interannual to multidecadal timescales, which amplify the decadal‐to‐multidecadal variability of NPP and RH resulting from GHG forcing in the CESM‐LME (Figures 1 and 7c).

Figure 6

Power spectra for time‐series of global annual net biome production (NBP), net primary production (NPP), heterotrophic respiration (RH), and disturbance (D), as simulated by the Community Earth System Model Last Millennium Ensemble (CESM‐LME) in the single‐forcing runs (greenhouse gases [GHGs], land use and land cover change [LULCC], orbital, solar, and volcanic), as well as the 13 full‐forcing ensemble members and the 850‐year control simulation.

Figure 7

The responses in the global annual net biome production (NBP), net primary production (NPP), heterotrophic respiration (RH), and disturbance (D), and mean land surface temperature (T), precipitation (P), and shortwave downwelling radiation (R) to 18 large volcanic eruptions as simulated in the Community Earth System Model Last Millennium Ensemble (CESM‐LME) full‐forcing simulations. The left panel shows the anomalies in global C fluxes and in mean T during the periods of five years before and after the volcanic eruption event year. The colored curves are shown for four significant eruptions, and the black curve with gray region shows the mean across the 18 eruptions ± one standard deviation The right panel shows scatter‐plots of the anomalies changes during the overall volcanic eruption events (summed from event year to three years after the event year) against the strength of volcanic eruptions (measured in estimated aerosol emissions).

The modeled global NBP shows significant responses to the series of large volcanic eruptions during the PILM (Figure 1a). We selected the largest 18 volcanic eruptions, of which the total reconstructed stratospheric volcanic sulfate aerosol injections (hereafter defined as the strength of volcanic eruptions) are more than 30 Tg (the amount of the 1991 Mount Pinatubo eruption; Gao et al., 2008). At interannual timescales, we found that the duration of the fast responses in land C fluxes to these volcanic eruptions is approximately limited to 3 years (Figure 7). On the year of the volcanic eruption (event year), the modeled NBP presents a positive anomaly (i.e., a larger sink; Figure 7a) mainly because of negative RH anomalies due to decreased surface/soil temperature (Figures 7c and 7e); because volcanic aerosols in the stratosphere typically reduce the direct solar radiation reaching the lower atmosphere and thus cool the land surface. Neither the modeled global NPP nor D show notable anomalies on the year of the event but show large anomalies in the next 2 years. In the first year after these eruptions (Figures 7b–7d), the NPP and RH show a negative response compared to a positive response in D, and the magnitudes of these responses increase with the strength of volcanic eruptions (Figure 7. shows Tambora in Indonesia: 109.7 Tg in 1815, Kuwae of Vanuatu: 137.5 Tg in 1452, and Samalas in Indonesia: 258 Tg in 1258). The largest eruption in 1258 (the eruption event year) resulted in decreases in the next year (1259) after the eruption event year of about 3.5 and 2.5 Pg C/year‐1 in NPP and RH, respectively, and an increase of about 0.8 Pg C/year‐1 in D, in response to decrease of 1.6 °C in land surface temperature in 1259. The modeled global NBP consequently displays a clear negative anomaly (‐1.8 Pg C/year‐1) in 1259 (Figure 7a). The increased D (Figure 7d) in response to large volcanic eruptions may be because such large volcanic aerosol forcings can change atmospheric circulation patterns, which resulted in decreased precipitation in tropical semi‐arid regions and East Asia (see Figure 8e for the Little Ice Age period) and promoted drought and fires frequencies (e.g., Zhuo et al., 2014).

Figure 8

Differences between the Little Ice Age (LIA; 1450‐1849) and the Medieval Climate Anomaly (MCA; 950‐1250). Relative changes (%) in (a) total ecosystem carbon pool, (b) net primary production (NPP), and (c) ecosystem turnover rate, and changes in (d) mean annual surface temperature, (e) annual precipitation, and (f) shortwave downwelling radiation as simulated by the Community Earth System Model Last Millennium Ensemble (CESM‐LME) full‐forcing simulations.

We evaluated accumulated anomalies from volcanic eruption event years compared to the subsequent 2 years (Figures 7a‐7d on the right panel) and found that of the 18 eruptions, five of them show significant positive anomalies in NBP (0.4 to 1.2 Pg C), while nine of them including three major eruptions present significant negative anomalies in NBP (‐0.4 to ‐1.6 Pg C). Furthermore, there was no significant correlation between the NBP response and the strength of volcanic eruptions (R 2 = 0.27, p = 0.28). Negative correlations were found between the strength of volcanic eruptions and NPP and RH (R 2 = 0.6, p<0.001), and between volcanic source strength and the corresponding anomalies in temperature (R 2 ≈ 0.56, p<0.001), precipitation (R 2 ≈ 0.66, p<0.001), and shortwave radiation (R 2 ≈ 0.68, p<0.001; Figures 7e‐7g). Further, we calculated the sensitivity of C fluxes to volcano‐induced cooling in the 3 years following the 18 eruptions (by calculating ΔNBP/ΔT, ΔNPP/ΔT, and ΔRH/ΔT for each eruption). These values are on average all negative for NBP (‐0.62 Pg C · year‐1 · °C‐1), NPP (‐0.66 Pg C · year‐1 · °C‐1), and RH (‐0.77 Pg C · year‐1 · °C‐1). The NBP shows comparable temperature sensitivity as seen for NPP and RH, suggesting that temperature variations are the main driver of short‐term changes in terrestrial carbon fluxes by volcanic eruptions.

Differences in the Land Carbon Cycle Between the Medieval Climate Anomaly and the Little Ice Age

Over the last millennium, the Earth's climate experienced a long period of cooler temperatures from 1450 to 1849, often referred as the Little Ice Age (LIA), which was preceded by a period of warmer temperatures from 950‐1250, the Medieval Climate Anomaly (MCA). Causes of the LIA are not well identified, possibly being related to volcanic eruptions (Crowley, 2000), decreased solar activity (Crowley, 2000), orbital parameter changes (Kaufman et al., 2009), altered ocean circulation (Wanamaker et al., 2012), and LULCC from early deforestation (Ellis et al., 2013). Nevertheless, analyzing differences between the warm MCA and the cool LIA in the ESM provides some clues to help understand the warming effects of climate change on the future carbon cycle.

By examining the mean states of the land carbon cycle and climate, comparing the warm MCA (on average of 300 years) and the cool LIA (on average of 400 years), we can estimate long‐term sensitivities of land carbon fluxes to climate and external forcing (e.g., volcanic aerosols, LULCC). The CESM‐LME was found to provide a reasonable representation of the cooling period of the LIA compared with historical proxy‐based reconstructions (Otto‐Bliesner et al., 2016). From the full‐forcing simulations, we found that during the cool LIA, the lower surface temperature mostly occurred over the Northern Hemisphere (NH; north of 30°N), with about ‐0.15 to ‐0.45 °C smaller than the MCA at regional scales (Figure 8d). Total land carbon storage (C eco) over the tropics was not significantly (0.3%) different between the two periods: the LIA and the MCA, as both tropical NPP and k eco varied only slightly during the LIA (Figure 8a‐8c). However, over the NH, the C eco during the LIA is 1.6% larger than during the MCA (Figure 8a). Regionally, over eastern European and Canadian forests, C eco in the LIA is 4.8% larger than that in the MCA, which resulted from both increased vegetation and soil C pools (Figure S6), because during the LIA, the 400‐year averaged NPP has increased (Figure 8b) and relative small changed turnover rate (Figure 8c). But over western Europe, which is mainly covered by croplands at current climate, both C eco during the LIA is much (‐5.5%) smaller than that in the MCA resulting from much smaller averaged NPP and higher turnover rate (larger soil C decomposition; Figures 8a‐8c and S6). For the tundra in Northern Russia (across the Arctic Circle), surface temperature during the LIA is about ‐0.5°C lower than the MCA (Figure 8d), resulting in depressed plant growth, which led to decreased NPP and C eco (by 10%) and a smaller ecosystem turnover rate k eco (Figure 8c).

The CESM‐LME single‐forcing experiments can help identify which external forcings dominate regional changes in climate variables and terrestrial carbon cycle between the LIA and MCA. Results show that changes in temperature between the MCA and the LIA are mostly attributable to the cumulative effect of several huge volcanic eruptions during the LIA (Figure S7f), followed by orbital and land‐use changes (Figure S7). The increased C eco and NPP in the LIA over the eastern European and Canadian forests may result from smaller plant and soil respiration due to lower temperature. Over western Europe, the modeled smaller C eco and NPP in the LIA (Figures 8a and 8b) may be due to reduced plant productivity from deforestation during the LIA (Pongratz et al., 2008; Pongratz et al., 2009). A decrease in C eco was accompanied with an increased turnover rate k eco (RH/Ceco) over western Europe (Figure 8c). The large changes in shortwave radiation over western Europe (Figure 8f) may be explained by large‐scale land‐atmosphere interactions and vegetation geophysical feedback resulting from major deforestation (Figure S9c); this caused increased land surface albedo, cooled regional climate, and altered general atmospheric circulation (e.g., increased cloud cover) during the LIA 400‐year period (Pongratz et al., 2009; Pongratz et al., 2010). The stratospheric aerosols from volcanoes may also contribute the reduction in shortwave radiation (Figure S9f). Overall, using the CESM‐LME, remarkable changes were found in the land carbon cycle over the NH during the LIA, especially across the boreal forests and tundra regions, due to large decreased surface temperature caused by huge volcanic eruptions during the LIA.

Discussion

Contributions of Carbon Fluxes to Variability in NBP in Terrestrial Ecosystems

Over the PILM (850‐1849), we found that the global terrestrial ecosystem carbon balance remained close to an equilibrium, with global NBP being close to zero at centennial timescales, which is mostly in agreement with estimates based on atmospheric 13C in CO2, which suggest a very small increase of only about 30 Pg C/ka between AD 755 and 1850 in land pools (Bauska et al., 2015). This equilibrium was preserved mainly because the carbon turnover rate of the global terrestrial ecosystem was not significantly reduced or enhanced from natural climate variability and external forcings for the PILM (Figure 3). Based on our theoretical framework for the variance of key components of the C cycle at different time‐scales (equations (1) and (2); Figure 9a; Zhang et al., 2018), the CESM‐LME modeling results show the preindustrial dynamics of terrestrial carbon cycle; in particular, it was seen that variation in the terrestrial net C balance (i.e., NBP) is dominated by fluctuations in C input (i.e., NPP) on the interannual to decadal timescales, while the influence from output fluxes (i.e., RH+D) tends to increase with increasing period and becomes dominant at centennial timescales (100‐1,000 years). The D flux contributes relatively little to variability in NBP on centennial timescales, which mostly arose from early anthropogenic land use change. When the ecosystem C turnover remains steady on centennial timescales, variability in RH is mostly driven by variation in NPP via the contribution of plant litter carbon to the soil C pool.

Figure 9

The relative variance in net primary production (NPP), net biome production (NBP), and heterotrophic respiration (RH)+disturbance (D) across different timescales and turnover rates. (a) Spectra for the theoretical estimates of amplitude ratios A NBP/A NPP and A RH+D/A NPP across a wide range of periods (1 to 1,000 years) for selected turnover rates (k eco = 0.01, 0.1, 1.0 year‐1) of land carbon cycle. Variance ratios of (b) A NBP/A NPP and (c) A RH+D/A NPP as a function of the period (x axis) and residence times (y axis). The variance ratios were calculated following equations (1) and (2).

Theoretically, as is shown in Figure 9, when an ecosystem is shifted to a new equilibrium state with a higher turnover rate (or shorter carbon residence time, e.g., k eco increasing from 0.01 to 0.1 year‐1 in Figure 9), as the fluctuations in NPP on a given timescale will be more offset by variation in RH of opposite sign, with minimal impact on the variability of NBP on this timescale, it implies that the rate of C sequestration would decline during the transition period (Luo & Weng, 2011); thereby, the newly input carbon into ecosystem carbon pool is retained for a shorter timescale (Figure 9). Factors affecting the change of C turnover rate include changes to ambient temperature or soil moisture. Climate changes impacting upon RH exist on various timescales, including such as the decomposition of the older soil carbon substrates that were accumulated on centennial to longer timescales. For example, the warm MCA, as modeled by the CESM‐LME, shows higher turnover rate and smaller carbon stock over boreal temperate forests than the cool LIA (Figure 8). Because of the larger climatological temperature sensitivity of soil organic carbon and larger soil organic carbon storage in cold‐climate areas (Koven et al., 2017), changes in the carbon turnover rate due to climate change may become more important over cold and temperate latitudes than the tropical regions on centennial timescales.

Contributions of Climate Variability and External Forcings to Variability in NBP

We find that climate‐internal variability dominates the variability in NBP on interannual to multidecadal timescales, which is consistent to previous studies (Bauska et al., 2015; Sitch et al., 2015; Zhang et al., 2018). The spectral analysis of NBP over 850‐1849 shows strong periodic peaks of around 4 years and 50 years. A similar result was found for interannual variability in another ensemble of land carbon models (TRENDY), which were forced by observational climate data over 1901‐2010, albeit over a much shorter period (Zhang et al., 2018). The spectrum for the TRENDY ensemble demonstrates periodic peaks at about 2.5, 4, and 30 years with much smaller peak values (e.g., about 6 (PgC)2/year‐2 at the period of 4 years) compared to the CESM‐LME NBP spectrum (about 16 (Pg C)2/year‐2 at period of 4 years). These peaks were attributed to ENSO and Pacific Decadal Oscillation climate modes controlling land temperature and precipitation, respectively, suggesting that climate modes modulate the short‐term variance in NBP by temperature and precipitation (Zeng et al., 2005; Zhang et al., 2018); the TRENDY ensemble spans 1901‐2010; thus, the differences between the TRENDY and CESM‐LME spectra may reflect differences in forcings (e.g., from GHGs), sampling error (given the short period of TRENDY), or the more accurate constraints on forcings imposed during the late twentieth century. However, the overestimated ENSO‐like periodic peak values in the power spectrum of NBP apparent in the CESM‐LME output can be related to the fact that the amplitude of the ENSO was overestimated in the CESM compared with observations (Otto‐Bliesner et al., 2016), indicating that large uncertainties remaining in modeling ocean‐atmosphere interactions could translate into large uncertainties in the interannual variability in regional/global terrestrial carbon cycle.

Let us consider how the climate system's external forcings influence the terrestrial carbon cycle under the background climate before the 1850s. To our surprise, the analysis shows that even during the PILM, GHG forcings (subject to both natural feedback and early anthropogenic activity) are identified as the dominant driver of terrestrial carbon cycle on centennial timescales; this forcing explained about 50% of variations in NPP and RH on centennial timescales in the CESM‐LME. While in the real world, the causality operates in both directions (the C pools of the terrestrial biosphere and the atmospheric CO2 levels influencing one another), the ensemble used specified GHGs and thus variation in the biosphere did not impact upon the assumed CO2 levels. Although the GHGs and the LULCC are two largest contributors to the long‐term land carbon variability during the industrial period (Le Quéré et al., 2009;, 2018; Schimel et al., 2015), the human‐induced land cover change (only forcing from change in area) contributed relatively little to variations or the trend in the land C sink before the industrial period on centennial timescales, as estimated by previous studies (Kaplan et al., 2010).

Other studies have shown that volcanic eruptions like the 1991 Mount Pinatubo eruption significantly increased annual global land C uptake during 1 to 2 years after the eruptions (Ciais et al., 1995; Jones & Cox, 2001; Le Quéré et al., 2014; Mercado et al., 2009), while the CESM‐LME results suggest that eruptions larger than the 1991 eruption during the PILM caused a sharp increase in land C uptake on the year of the event, followed by a reduction in NBP over the next 2 years. The volcanic eruptions appeared to lower both NPP and RH due to cooling of the land surface, but an increase in D (e.g., modeled fires) most likely due to increased drought and decreased precipitation from altered atmospheric circulations. For example, it has been shown that monsoonal rainfall over China following volcanic eruptions was reduced in response to volcanic aerosols over the past 700 years (Zhuo et al., 2014). This short‐term effect on land C uptake generally lasts 3 years, and the accumulated impacts over this period from huge eruptions are still very uncertain (Figure 7). However, the effect of large‐scale massive volcanic eruptions, as propagated by slower climate‐carbon feedback represented in the CESM over the last millennium, are found to strongly increase terrestrial C sink (such as increased carbon stock over the NH during LIA caused by the cumulative effect of several huge volcanic eruptions). Our modeling analysis suggests that the cooling during the LIA was largely contributed by a series of huge volcanic eruptions over the LIA period. The cooling‐induced terrestrial C uptake was consider as the cause of the observed lower ice‐core CO2 level record during the LIA (Rubino et al., 2016). These long‐term impacts were evident over timescales of around few hundred years until the land C cycle returns to equilibrium and contributed 25% of variance in NPP and RH on centennial timescales (Figure 6). This long‐term (short‐term) volcanic effect on terrestrial ecosystem can drive atmospheric CO2 changes of more than 15 years (within 3 years) in another ESM (MPI‐ESM; Giorgetta et al., 2013; Jungclaus et al., 2010), which was mostly due to changes in C pools across the tropical and subtropical regions, where RH was most significantly reduced in response to the surface cooling after eruptions (Brovkin et al., 2010; Lehner et al., 2015). However, the simulated impacts of volcanic eruptions on the C cycle are still model‐dependent (Lehner et al., 2015).

This study also has a bearing on various geoengineering techniques proposed to mitigate some of the effects of anthropogenic climate change; such proposals require a high degree of scrutiny in terms of their likely impact upon all elements of the earth system. Injection of sulphate aerosols into the stratosphere (mimicking the effect of large volcanoes) has been canvassed as a means of solar radiation management (Vaughan & Lenton, 2011). The large volcanic eruptions during the PILM offer a natural laboratory for the study of such interventions. Our findings suggest that significant volcanic eruptions cause negative anomalies in both temperature and NBP over 2‐3 years, the latter of which may lead to an increase in ambient CO2 concentration and potentially increased warming at decadal to centennial timescales, although stronger conclusions could be drawn if the impact of increased diffuse radiation were incorporated into these simulations.

Our findings also have a bearing upon the likely impacts of a reduction in atmospheric CO2 by afforestation (via sequestration of soil C). We found that the global terrestrial C pool remained close to equilibrium at centennial timescales. While the afforestation impact may be significant for regional NBP, our results support taking a conservative view to the likely global effects over very long timescales, especially when one considers the challenge of maintaining such reserves over centuries.

Uncertainties in Modeling the Terrestrial Carbon Cycle in an ESM

In this study, although ensemble simulations from one model were performed to reduce uncertainties from model initialization, uncertainties represented by intermodel comparisons remain important; we now list some of the key uncertainties and how they may impact upon the findings. Firstly, the land surface model of the CESM (CLM4CN) was reported to underestimate the total soil C pool size compared to the observations and the other ESMs and thus overestimate its turnover rate especially over the NH (He et al., 2016). This would exacerbate uncertainties in our analysis of spatial patterns of C fluxes. Secondly, the CESM‐LME lacks the negative forcing from indirect aerosol effects, resulting in a overestimation of anthropogenic climate warming (Lehner et al., 2015). Thirdly, we should also note that the land model (CLM4CN) in CESM‐LME lacks the processes to describe the high latitude permafrost carbon reservoir, which contains a large amount of organic carbon matter frozen in the soil. Large areas of the permafrost are also frozen peatlands over such as Canada and Siberia. With temperatures increasing during the last 150 years, these areas may become a large source of CO2 and CH4, as the thawing of permafrost could result in microbial decomposition of frozen organic carbon (Schuur et al., 2008). However, there remains uncertainty over the current decomposition rate of frozen organic carbon (Dutta et al., 2006; Schuur et al., 2008). This positive feedback of soil carbon respiration to global warming is important on centennial to millennial timescales (Schuur et al., 2015). The sink of atmosphere CO2 from peatlands is estimated to be significant (Spahni et al., 2012), but the decomposition of peatland/permafrost organic carbon was not parameterized in the CESM‐LME. Additionally, model structural uncertainty remains an important issue, which determines the influence of the simulated internal climate variability and external forcings over the last millennium (Lehner et al., 2015). For instance, the CLM4‐CN in CESM‐LME is simplified with only one soil C layer that induces uncertainties and/or inaccuracies in the results of this work (Oleson et al., 2013). Bringing in more soil C layers in future CLM version will increase the heterogeneity in response times as, for example, temperature and soil moisture varies with soil depth, which would present more details of soil respiration. Volcanoes have strong effects on global temperature, precipitation, and the carbon cycle, but large differences between Earth system models exist due to different representations of the response of terrestrial vegetation to climate (Friedlingstein et al., 2014; Lehner et al., 2015).

Therefore, it would be useful if our results were tested further by additional analysis including more models (e.g., from the PMIP or the CMIP), especially those that include more biochemical or biophysical processes such as the peatland and permafrost. This must be contrasted by the additional difficulties in conducting such ensemble simulations (e.g., ensuring all key model inputs are identical) and in its interpretation (given the different assumptions, parameterizations, and components of the individual models). The last millennium simulations from the PMIP or the CMIP, for instance, could be used as climate forcing for DGVMs (e.g., from TRNEDY project). This would give a spread of uncertainties that catches both climate model uncertainty, and DGVM uncertainty.

Overall, although large uncertainties remain in the CESM to characterize the timescale‐dependent variability of terrestrial C cycle over the PILM, this study highlights the importance of terrestrial ecosystem C turnover on balancing the long‐term carbon storage in response to changing climate and external forcings, especially GHGs and volcanoes over centennial timescales.

Supporting Information S1

Click here for additional data file.