Climate and carbon cycle dynamics in a CESM simulation from 850 to 2100 CE

Abstract

Under the protocols of phase 3 of the Paleoclimate Modelling Intercomparison Project, a number of simulations were produced that provide a range of potential climate evolutions from the last millennium to the end of the current century. Here, we present the first simulation with the Community Earth System Model (CESM), which includes an interactive carbon cycle, that covers the last millennium. The simulation is continued to the end of the twenty-first century. Besides state-of-the-art forcing reconstructions, we apply a modified reconstruction of total solar irradiance to shed light on the issue of forcing uncertainty in the context of the last millennium. Nevertheless, we find that structural uncertainties between different models can still dominate over forcing uncertainty for quantities such as hemispheric temperatures or the land and ocean carbon cycle response. Compared to other model simulations, we find forced decadal-scale variability to occur mainly after volcanic eruptions, while during other periods internal variability masks potentially forced signals and calls for larger ensembles in paleoclimate modeling studies. At the same time, we were not able to attribute millennial temperature trends to orbital forcing, as has been suggested recently. The climate–carbon-cycle sensitivity in CESM during the last millennium is estimated to be between 1.0 and 2.1 ppm ${}^{\circ }$C${}^{-\mathrm{1}}$. However, the dependence of this sensitivity on the exact time period and scale illustrates the prevailing challenge of deriving robust constraints on this quantity from paleoclimate proxies. In particular, the response of the land carbon cycle to volcanic forcing shows fundamental differences between different models. In CESM the tropical land dictates the response to volcanoes, with a distinct behavior for large and moderate eruptions. Under anthropogenic emissions, global land and ocean carbon uptake rates emerge from the envelope of interannual natural variability by about year 1947 and 1877, respectively, as simulated for the last millennium.

Introduction

The last about 1000 years constitute the best opportunity previous to the instrumental period to study the transient interaction of external forcing and internal variability in climate, atmospheric CO${}_{\mathrm{2}}$, and the carbon cycle on interannual to multi-decadal timescales. In fact, the instrumental record is often too short to draw strong conclusions on multi-decadal variability. The relatively stable climate together with the abundance of high-resolution climate proxy and ice core data makes the last millennium an interesting target and test bed for modeling studies. However, the large and sometimes controversial body of literature on the magnitude and impact of solar and volcanic forcing on interannual to multi-decadal climate variability illustrates the challenges inherent in establishing a robust understanding of a period that is characterized by a small signal-to-noise ratio in many quantities and for which uncertainties in the external forcing remain see, e.g.,. In addition, a process-based quantitative explanation of the reconstructed preindustrial variability in atmospheric CO${}_{\mathrm{2}}$ and carbon fluxes is largely missing.

Compared to glacial–interglacial climate change, the last millennium experienced little climate variability, yet there is evidence for distinct climate states during that period e.g.,. Within the last millennium, the Medieval Climate Anomaly (MCA, $\sim$ AD 950–1250) and the Little Ice Age (LIA, $\sim$ AD 1400–1700) are two key periods of documented regional or global temperature excursions suggested to be driven by a combination of variations in solar irradiance and volcanic activity see, e.g., . Despite great efforts in reconstructing e.g., and simulating e.g., the transition from the MCA to the LIA, substantial uncertainties remain with respect to the mechanisms at play. Recent studies point towards solar insolation playing a minor role with regard to climate over the last millennium , while regional feedback processes in response to volcanic eruptions and solar variability need to be considered to explain decadal-scale climate variability e.g.,. At high northern latitudes, the importance of millennial-scale orbital forcing is another debated issue e.g.,.

Further, the last millennium offers the possibility to study the natural variability in the carbon cycle and its response to external forcing. Models with a carbon cycle module are extensively tested against present-day observations and widely used for emission-driven future projections see, e.g.,. However, there are only few last-millennium simulations including a carbon cycle see, e.g.,. The sensitivity of the carbon cycle to climate has been shown to be mostly positive, i.e., with warming, additional CO${}_{\mathrm{2}}$ is released to the atmosphere . However, the magnitude of this feedback remains poorly constrained by observations and models e.g.. In particular, determining the role of land in past and future carbon cycle variability and trends is still challenging. In both idealized and scenario-guided multi-model studies , land constitutes the largest relative uncertainty in terms of intermediate- to long-term carbon uptake.

As for physical climate quantities, explosive volcanic eruptions constitute an important forcing for the carbon cycle. The sensitivity of the carbon cycle to such eruptions has been investigated by , , , and using different Earth system models. For this short-lived forcing, the land response appears to be the driver of most post-eruption carbon cycle changes, with a range of magnitudes and time horizons associated with the different models. Further, pointed out that the magnitude of the carbon cycle response to volcanoes depends critically on the climate state during the eruption.

The third Paleoclimate Modelling Intercomparison Project PMIP3; and fifth Coupled Model Intercomparison Project CMIP5; represent joint efforts, in which different modeling groups perform identical experiments, allowing for a systematic comparison of the models e.g.,. Here we contribute to the existing set of simulations an integration from 850–2100 CE with the Community Earth System Model, including a carbon cycle module. The aims of this study are (i) to detect coherent large-scale features of forced variability in temperature and carbon cycle quantities, in particular in response to volcanic eruptions, (ii) to investigate the relative role of forcing uncertainty and model structural uncertainty, and (iii) to provide an estimate of preindustrial variability and the time of emergence from it under anthropogenic climate change. The setup chosen here is unique in a number of ways and tailored to address the abovementioned aims: first, the carbon cycle is fully interactive with the other model components with the exception of the radiation code, which is fed by reconstructed CO${}_{\mathrm{2}}$. This allows us to study the isolated effect of climate on the carbon cycle, while guaranteeing an external forcing consistent with existing reconstructions. Second, the orbital parameters are held constant to study their importance relative to simulations with transient orbital parameters. Third, the solar forcing incorporated into the simulation has a larger amplitude than the majority of PMIP3 simulations and hence enables us to investigate whether the results are sensitive to this amplitude.

This paper is structured as follows: a description of the model and experimental setup is presented in Sect. 2. In Sects. 3 and 4 we address the general simulated climate and carbon cycle evolution and investigate forced and unforced variability in the simulated climate by comparing models to reconstructions and models to models. Section 5 focuses on the response of the climate and carbon cycle to volcanic forcing. Section 6 deals with estimating the climate–carbon-cycle sensitivity in the Community Earth System Model (CESM). A discussion and conclusions follow in Sect. 7.

Data and methods

Model description

The Community Earth System Model CESM; is a fully coupled state-of-the-art Earth system model developed by the National Center for Atmospheric Research (NCAR) and was released in 2010. In terms of physics, CESM relies on the fourth version of the Community Climate System Model CCSM4;. Additionally, a carbon cycle module is included in CESM's atmosphere, land, and ocean components. The CESM version used here is release 1.0.1 in the 1${}^{\circ }$ version and includes components for the atmosphere, land, ocean, and sea ice, all coupled by a flux coupler.

The atmospheric component of CESM 1.0.1 is the Community Atmosphere Model version 4 CAM4;, which has a finite volume core with a uniform horizontal resolution of 1.25${}^{\circ }$ $\times$ 0.9${}^{\circ }$ at 26 vertical levels. The land component is the Community Land Model version 4 CLM4;, which operates on the same horizontal grid as CAM4 and includes a prognostic carbon–nitrogen cycle that calculates vegetation, litter, soil carbon, vegetation phenology, and nitrogen states. Further, it includes a prognostic fire module, which is governed by near-surface soil moisture conditions and fuel availability.

The ocean component is the Parallel Ocean Program version 2 POP2; with an nominal 1${}^{\circ }$ horizontal resolution and 60 depth levels. The horizontal resolution varies and is higher around Greenland, to where the North Pole is displaced, as well as around the equator. Embedded in POP2 is the Biogeochemical Elemental Cycle model BEC; that builds on a nutrient–phytoplankton–zooplankton–detritus food web model and distinguishes three phytoplankton functional types . Carbon export and remineralization are parameterized according to . Alkalinity, pH, partial pressure of CO${}_{\mathrm{2}}$, concentrations of bicarbonate, and carbonate ions are diagnosed from prognostic dissolved inorganic carbon, alkalinity, and temperature- and salinity-dependent equilibrium coefficients. Biogenic calcification is implemented as being proportional to a fraction of small-phytoplankton production, which is temperature-dependent. An exponential curve is prescribed to simulate the dissolution of sinking CaCO${}_{\mathrm{3}}$ . There is no dependence of calcification–dissolution rates on saturation state. Organic material reaching the ocean floor is remineralized instantaneously, i.e., no sediment module is included. River discharge from CLM4 does not carry dissolved tracers but nitrogen deposition to the ocean surface has been prescribed. The sea ice component is the Community Ice Code version 4 (CICE4) from the Los Alamos National Laboratories , including elastic viscous plastic dynamics, energy-conserving thermodynamics, and a subgrid-scale ice thickness distribution. It operates on the same horizontal resolution as POP2.

List of simulations conducted for this study. See text for details regarding the forcing. TSI: total solar irradiance; GHGs: greenhouse gases; $E_{{\mathrm{CO}}_{\mathrm{2}}}$: anthropogenic CO${}_{\mathrm{2}}$ emissions from fossil fuel burning and cement production; LULUC; land use and land use change.

	Control simulation (CTRL)	Transient simulation (CESM)
Forcing	850 CE (500 years)	850–2099 CE
TSI	1360.228 W m${}^{-\mathrm{2}}$	adjusted and
Volcanic	none
GHGs	CO${}_{\mathrm{2}}$ (279.3 ppm)
	CH${}_{\mathrm{4}}$ (674.5 ppb)
	N${}_{\mathrm{2}}$O (266.9 ppb)
$E_{{\mathrm{CO}}_{\mathrm{2}}}$	none	and
Aerosol	1850f CE from
Orbital	1990 CE after	1990 CE after
LULUC	850 CE from	and

Experimental setup

Table  provides an overview of the simulations conducted for this study. First, a 500-year control simulation with perpetual 850 CE forcing (hereafter CTRL) was branched off from an 1850 CE control simulation with CCSM4 . However, restart files for the land component were taken from an 850 CE control simulation, kindly provided by the NCAR, in which the land use maps by were applied. This procedure has the advantage that the slowly reacting soil and ecosystem carbon stocks are closer to 850 CE conditions than in the 1850 CE control simulation. A transient simulation covering the period 850–2099 CE was then branched off from year 258 of CTRL. Despite the shortness of CTRL leading up to the start of the transient simulation, most quantities of the surface climate, such as air temperature, sea ice, or upper-ocean temperature, can be considered reasonably equilibrated at the start of the transient simulation, as the forcing levels due to total solar irradiance (TSI) and most greenhouse gases are similar between 1850 and 850 CE . However, weak trends in CTRL are still detectable in slowly reacting quantities, such as deep-ocean temperature (below 2000 m; $\sim -$0.04 ${}^{\circ }$C 100 yr ${}^{-\mathrm{1}}$), Atlantic meridional overturning circulation ($\sim -$0.22 Sv 100 yr ${}^{-\mathrm{1}}$), Antarctic Circumpolar Current ($\sim$ 0.70 Sv 100 yr ${}^{-\mathrm{1}}$), dissolved inorganic carbon in the ocean ($\sim$ $-$0.01 % 100 yr ${}^{-\mathrm{1}}$), or soil carbon storage ($\sim$ 0.2 % C 100 yr ${}^{-\mathrm{1}}$). The Antarctic Bottom Water formation rate shows no drift.

Forcings used in the last millennium simulation with CESM. (a) TSI in comparison with the different TSI reconstructions proposed by PMIP3. (b) Volcanic forcing as total volcanic aerosol mass. (c) Radiative forcing RF, calculated according to from the greenhouse gases CO${}_{\mathrm{2}}$, CH${}_{\mathrm{4}}$, and N${}_{\mathrm{2}}$O. (d) Major changes in land cover (as fraction of global land area). See text for details.

[Figure omitted. See PDF]

The applied transient forcing largely follows the PMIP3 protocols and the Coupled Model Intercomparison Project 5 CMIP5;, consisting of TSI, greenhouse gases (GHGs), volcanic and anthropogenic aerosols, and land use changes (Fig. ). Here, the TSI reconstruction by TSI${}_{\mathrm{VS}10}$ is used, to which a synthetic 11-year solar cycle is added . In light of the recently enlarged envelope of reconstructed TSI amplitude , we scale TSI by a factor of 2.2635 to have an amplitude of 0.2 % between present day (1961–1990 CE) and the late Maunder Minimum (1675–1704 CE), which is about twice as large as the 0.1 % used in most PMIP3 simulations: $$\mathrm{TSI}=2.2635\cdot ({\mathrm{TSI}}_{\mathrm{VS}10}-\stackrel{\mathrm{‾}}{{\mathrm{TSI}}_{\mathrm{VS}10}})+\stackrel{\mathrm{‾}}{{\mathrm{TSI}}_{\mathrm{VS}10}}.$$ Figure a shows that the TSI used here lies in between the large-amplitude reconstruction by and the bulk of small-amplitude reconstructions of the original PMIP3 protocol . Note that a recent detection and attribution study indicates that small-amplitude TSI reconstructions agree better with temperature reconstructions over the last millennium than large-amplitude reconstructions , in agreement with . For the twenty-first century, the last three solar cycles of the data set are repeated continuously. The insolation due to Earth's orbital configuration is calculated according to with the orbital parameters held constant at 1990 CE values.

The volcanic forcing follows from 850 to 2001 CE. It provides estimates of the stratospheric sulfate aerosol loadings from volcanic eruptions as a function of latitude, altitude, and month and is implemented in CESM as a fixed single-size distribution in the three layers in the lower stratosphere . Post-2001 CE volcanic forcing remains zero.

Land use and land use changes (LULUC) are based on from 850 to 1500 CE, when this data set is splined into , a synthesis data set that extends into the future. The two data sets do not join smoothly but exhibit a small stepwise change in the distribution of crop land and pasture at the year 1500 CE. Up until about 1850 CE global anthropogenic LULUC are small; however, they can be significant regionally . Approaching the industrial era, LULUC accelerate, dominated by the expansion of crop land and pasture. Here, only net changes in land use area are considered. The impact of shifting cultivation and wood harvest on carbon emissions from land use is neglected; these processes are estimated to have contributed around 30 % to the total carbon emissions from land use .

The temporal evolution of long-lived greenhouse gases (GHGs; CO${}_{\mathrm{2}}$, CH${}_{\mathrm{4}}$, and N${}_{\mathrm{2}}$O) is prescribed based on estimates from high-resolution Antarctic ice cores that are joined with measurements in the mid-twentieth century and references therein. While the carbon cycle module of CESM interactively calculates the CO${}_{\mathrm{2}}$ concentration originating in land use changes, fossil-fuel emissions post-1750 CE, following, and carbon-cycle–climate feedbacks, it is radiatively inactive. Instead, the ice core and measured data are prescribed in the radiative code, keeping the physical model as close to reality as possible. As a result, the impact of the interactive coupling of the carbon cycle module is minor for simulated climate and limited to changes in surface conditions due to changing vegetation. For the extension of the simulation post-2005 CE the Representative Concentration Pathway 8.5 (RCP 8.5) is used, representing the unmitigated “business-as-usual” emission scenario, corresponding to a forcing of approximately 8.5 W m${}^{-\mathrm{2}}$ in the year 2100 .

Aerosols such as sulfate, black and organic carbon, dust, and sea salt are implemented as non-time-varying up to 1850 CE, perpetually inducing the spatial distributions of the 1850 CE control simulation during this time . Post-1850 CE, the time-varying aerosol data sets provided by are used, although CESM only includes a representation of direct aerosol effects. Similarly, nitrogen (NH${}_x$ and NO${}_y$) input is held constant until it starts to be time-varying from 1850 CE onwards, also following . Iron fluxes from sediments are held fixed .

Other model simulations

In addition to comparing CESM results to output from current Model Intercomparison Projects, we also compare them to results from a similar simulation with CCSM4 and IPSL-CM5A-LR (Institut Pierre Simon Laplace–climate model 5A–low resolution) , two simulations without interactive carbon cycle. The solar and volcanic forcing reconstructions applied to CCSM4 and IPSL-CM5A-LR are identical to ours with the exception of the scaling of TSI that we applied to CESM. The goal here is to investigate the question of whether different solar forcing amplitudes applied to the same physical model (CESM vs. CCSM4) have a larger effect than applying the same solar forcing to two different physical models (IPSL-CM5A-LR vs. CCSM4).

Selected forcing details of simulations used in comparisons with CESM. TSI: total solar irradiance; LULUC: land use and land use change.

	CESM	CCSM4	IPSL-CM5A-LR	MPI-ESM
Forcing
TSI	adjusted
	and	and	and
Volcanic
Orbital	1990 CE,	transient;	transient;	transient;
LULUC			non-transient
	and	and

Further, we compare CESM to MPI-ESM (Max-Planck-Institut Earth System Model) ECHAM5/MPIOM;, a model that includes an interactive carbon cycle, to assess the robustness of the simulated climate and carbon cycle variations in response to external forcing. MPI-ESM uses as TSI forcing and as volcanic forcing. These differ from the CESM forcing in amplitude much more than in timing and therefore allow for a comparison of the forced response. All simulations, except for ours, apply transient orbital forcing and are summarized in Table . If not using the full ensemble of MPI-ESM, we focus on the member “mil0021”. Another difference in terms of experimental setup between CESM and MPI-ESM is that MPI-ESM was run with a fully interactive carbon cycle, i.e., the prognostic CO${}_{\mathrm{2}}$ interacts with the radiation and through that again influences climate, while in our setup this is a one-directional interaction only. Further, MPI-ESM is more coarsely resolved than CESM in both ocean and atmosphere and applies the A1B scenario for the twenty-first century , which corresponds roughly to the current intermediate scenario as compared to the high-scenario RCP8.5 used in CESM.

(a) Northern Hemisphere and (b) Southern Hemisphere temperature anomalies in model simulations and reconstructions. The anomalies are with reference to 1500–1899 CE (left panels) and 1850–1899 CE (right panels). Gray shading in (a) indicates the reconstruction overlap ; in (b), it indicates the reconstruction by . The 5–95 % range of the simulations from the third Paleoclimate Modelling Intercomparison Project (PMIP3) and the fifth Coupled Model Intercomparison Project (CMIP5; applying the RCP8.5) are given in green and red shading, respectively. Note that MPI-ESM applies the A1B scenario , which has a weaker forcing than RCP8.5. Hemispheric means from observations are shown as thick black line . All time series have been smoothed by a local regression filter which suppresses variability higher than 30 years. The Medieval Climate Anomaly (MCA) and the Little Ice Age (LIA) are indicated as defined in . (c) Evolution of atmospheric CO${}_{\mathrm{2}}$ in CESM (black), MPI-ESM (grey; ensemble range), from ice cores (red), from measurements (orange), and from RCP8.5 used to force the radiative code in CESM (magenta). The small inset in the middle panel shows the observed annual cycle at Mauna Loa, Hawaii, and a 2${}^{\circ }$ $\times$ 2${}^{\circ }$ average over Hawaii from CESM, both derived from the period 1958–2012.

[Figure omitted. See PDF]

(a) Mean June–August (JJA) Arctic ($>$ 60${}^{\circ }$ N land) solar insolation in CCSM4 with time-varying orbital parameters and CESM with fixed orbital parameters. (b) Arctic JJA temperature difference between CCSM4 and CESM. The least-squares linear trend of this temperature difference is given in red. (c) Arctic JJA temperature anomalies (from their 850 to 1850 AD mean) versus solar insolation as 100- and 200-year averages (10 and 5 circles, respectively) from CCSM4 and CESM (red and blue, respectively). The least-squares linear trend for each cloud of 100- and 200-year averages is given in the respective color. The shading envelops the range of temperature versus solar insolation for each cloud of means.

[Figure omitted. See PDF]

General climate and carbon cycle evolution

Temperature

The simulated annual mean Northern Hemisphere (NH) surface air temperature (SAT) follows the general evolution of proxy reconstructions: a warm Medieval Climate Anomaly (MCA, $\sim$ 950–1250 CE), a transition into the colder Little Ice Age (LIA, $\sim$ 1400–1700 CE), followed by the anthropogenically driven warming of the nineteenth and twentieth centuries (Fig. ). The NH MCA-to-LIA cooling amounts to 0.26 $\pm$ 0.18 ${}^{\circ }$C taking the time periods defined above, which are as in, placing it in the lower half of reconstructed amplitudes that range from about 0.1 to 0.7 ${}^{\circ }$C . The subsequent warming between 1851–1880 CE and 1981–2010 CE amounts to 1.23 $\pm$ 0.15 ${}^{\circ }$C, while observations report only 0.71 $\pm$ 0.13 ${}^{\circ }$C . This overestimation by CESM takes place almost entirely after 1960 and arises largely from missing negative forcing from the indirect aerosol effect, which is not implemented in CAM4 . The late twentieth century being the warmest period in the NH in the past millennium is consistent with reconstructions e.g.,.

In CESM, the inception of the NH LIA occurs in concert with decreasing TSI and a sequence of strong volcanic eruptions during the thirteenth century. Reconstructions differ substantially in this matter and start to cool as early as 1100 CE or as late as 1400 CE. Further, new regional multi-proxy reconstructions of temperature provide no support for a hemispherical or globally synchronous MCA or LIA but show a clear tendency towards colder temperatures and exceptionally cold decades over most continents in the second half of the millennium .

The last millennium simulation with CCSM4 shows a largely coherent behavior with CESM in terms of amplitude and decadal variability in NH SAT (850–1850 CE correlation of 5-year filtered annual means $r=0.88$, $p<0.001$). The difference in NH SAT due to the different TSI amplitudes in CESM and CCSM4 scales roughly with the regression slope of NH SAT vs. TSI of both CESM and CCSM4 ($\sim$0.13 ${}^{\circ }$C per W m${}^{-\mathrm{2}}$), although internal variability can easily mask this effect at times. For example, the Maunder Minimum (1675–1704 CE), the 30-year period with the lowest TSI values and – when using TSI scaling – with the largest difference between CESM and CCSM4 (1.5 W m${}^{-\mathrm{2}}$), is only 0.14 ${}^{\circ }$C cooler than in CCSM4 and not 0.20 ${}^{\circ }$C as expected from the regression.

The NH temperature evolution of additional PMIP3 and CMIP5 simulations shows that the multi-model range is within the range of the reconstructions and encompasses the instrument-based observations (Fig. ). Disagreement between models and reconstructions exists in particular on the magnitude of response to the eruptions at 1258 CE and around 1350 CE. The 1258 CE eruption is the largest volcanic event recorded for the last millennium, and its climatic impact was likely enhanced through the cumulative effect of three smaller eruptions following shortly after . However, the pronounced cooling that is simulated by the models for this cluster of eruptions is largely absent in temperature reconstructions. Conversely, around 1350 CE temperature reconstructions show a decadal-scale cooling presumably due to volcanoes that is absent in the models, as the reconstructed volcanic forcing shows only two relatively small eruptions around that time. Part of this incoherent picture may arise from the unknown aerosol size distribution , the geographic location of past volcanic eruptions , and differences in reconstruction methods. As many proxy reconstructions of temperature rely heavily on tree ring data, it is worth noting that the dendrochronology community is currently debating whether the trees' response to volcanic eruptions resembles the true magnitude of the eruption .

Disagreement among the models exists on the relative amplitude of the MCA, where most models show colder conditions than CESM and CCSM4. Remarkably, the simulation by IPSL-CM5A-LR applied the same TSI and volcanic forcing as CCSM4, yet it comes to lie at the lower end of the PMIP3 model range during the MCA. In other words, the way in which models respond to variations in TSI and other forcings can still make a larger difference in the simulated amplitude than the scaling of TSI by a factor of 2, which in turn complicates a proper detection and attribution of solar forcing during the last millennium . Further disagreement among the models exists on the response to volcanic eruptions, where CESM and CCSM4 are among the more sensitive models an oversensitivity of CCSM4 to volcanoes based on twentieth-century simulations was reported by. Turning to the century-scale change over the industrial era, CESM and CCSM4 are at the upper end of the CMIP5 range and show an overestimation of the observed warming.

The simulated mean SAT of the Southern Hemisphere (SH) generally shows a similar evolution as in the NH, with the signature of the MCA and LIA superimposed on a weak millennial cooling trend. Models and reconstructions disagree to a larger extent in the SH than in the NH, in particular regarding cold excursions due to large volcanic eruptions, which are largely absent in the reconstructions. Similar results have been reported in a recent study on interhemispheric temperature variations that finds much less phasing of the two hemispheres in reconstructions than in models, potentially related to underestimated internal variability on the SH in models . A lingering question of climate modeling in general is whether models are too global in their response to external forcing. That is, they might show too little regional variability that is independent of the global mean response during a forced period. However, the uncertainties in the early period of the reconstructions make it impossible to robustly answer this question. Similar to the NH, the industrial warming in the SH from 1851–1880 to 1981–2010 CE (0.53 $\pm$ 0.07 ${}^{\circ }$C) is overestimated by CESM (0.71 $\pm$ 0.13 ${}^{\circ }$C).

The differential warming between the hemispheres in CESM is among the smallest among CMIP5 models (not shown). This is mainly due to the underestimated deep-water formation in the Southern Ocean, leading to a comparably strong warming of the SH and likely an underestimation of the oceanic uptake of anthropogenic carbon . With a transient climate response of 1.73 ${}^{\circ }$C and an equilibrium climate sensitivity of 3.20 ${}^{\circ }$C , CESM lies in the middle of recent estimates of 1.0 to 2.5 ${}^{\circ }$C and 1.5 to 4.5 ${}^{\circ }$C .

Orbital forcing

To detect and attribute the influence of orbital forcing on SAT trends during the last millennium, we compare our simulation with fixed orbital parameters to the CCSM4 simulation with time-varying orbital parameters (Fig. ). While both models experience a negative long-term trend in global TSI until about 1850 CE (Fig. ), the difference arising from the different orbital setup can be best seen in Arctic summer land insolation (Fig. ). Hence, Arctic summer land SAT has been proposed as a quantity that, on timescales of centuries to millennia, may be affected by orbital forcing . However, we find no detectable difference between the two simulations in the trend of Arctic summer land SAT (Fig. b). In fact, the Arctic multi-decadal to centennial summer land SAT anomalies in CESM span a very similar range as in CCSM4, despite CESM not accounting for time-varying orbital parameters: Fig. c shows non-overlapping 100- and 200-year mean SAT anomalies plotted against the corresponding mean solar insolation. The results from CCSM4 suggest a clear relationship between of the two quantities; however, the results of CESM show that nearly identical SAT anomalies are possible without orbital forcing. In other words, while we detect a long-term cooling trend in Arctic summer SAT in both CESM and CCSM4, we fail to attribute this trend to orbital forcing alone, as suggested by . This is confirmed in new simulations with decomposed forcing, again comparing simulations with fixed and time-varying orbital parameters (B. Otto-Bliesner, personal communication, 2014).

Cumulative carbon emissions in Pg C by different components over different time periods in CESM. Positive (negative) values indicate emission to (uptake from) the atmosphere.

	850–1500 CE	1501–1750 CE	1751–2011 CE	2012–2100 CE
Ocean	26.0	$-$4.0	$-$151.3	$-$413.0
Land	$-$15.0	10.3	82.5	139.4
Land (without LULUC)	$-$24.4	$-$9.3	$-$94.7	$-$436.3
Fossil Fuels	0.0	0.0	358.0	1851.5

Carbon cycle

The prognostic carbon cycle module in CESM allows us to study the response of the carbon cycle to transient external forcing. The land biosphere is a carbon sink during most of the first half of the last millennium, but becomes a source as anthropogenic land cover changes start to have a large-scale impact on the carbon cycle (Table ). The ocean is a carbon source at the beginning and becomes a sink in the second half of the last millennium. The residual of these fluxes represents changes in the atmospheric reservoir of carbon, illustrated in Fig. c by the prognostic CO${}_{\mathrm{2}}$ concentration. The amplitude of the simulated concentration does not resemble the one reconstructed from ice cores (i.e., imposed on the radiative code of CESM); in particular, the prominent CO${}_{\mathrm{2}}$ drop in the seventeenth century is not captured by CESM. This raises the question of whether the sensitivity of the carbon cycle to external forcing is too weak in CESM, whether the imposed land use changes are too modest , whether major changes in ocean circulation are not captured by models , and whether the ice core records are affected by uncertainties due to in situ production of CO${}_{\mathrm{2}}$ . Ensemble simulations with MPI-ESM also do not reproduce the reconstructed amplitudes or the drop . Further, Earth system models of intermediate complexity or vegetation models driven by GCM (general circulation model) output do not reproduce the uptake of carbon by either ocean or land needed to explain the reconstructed amplitudes .

Annual mean net carbon flux from the atmosphere to (a) land and (b) ocean. Green bars given the full and 10–90 % range from the preindustrial part of the simulation. Observational estimates are from .

[Figure omitted. See PDF]

The rise in atmospheric CO${}_{\mathrm{2}}$ due to fossil-fuel combustion is in good agreement with ice cores until about the 1940s. After that, a growing offset exists, leading to an overestimation of about 20 ppm by 2005 in CESM, qualitatively similar to the CMIP5 multi-model mean . From the observational estimates one can diagnose that the discrepancy arises primarily from overestimated carbon release from land Fig. a; see also. From 1750 to 2011 CE the cumulative total land release (including LULUC) is 83 Pg C compared to 30 $\pm$ 45 Pg C from observational estimates;, while the cumulative net land uptake is 95 Pg C 160 $\pm$ 90 Pg C;. The ocean cumulative uptake of 151 Pg C compares more favorably to current estimates of 155 $\pm$ 30 Pg C . Note, however, that given the overestimation of atmospheric CO${}_{\mathrm{2}}$, one would expect a higher ocean uptake. This bias originates largely in an underestimation of the uptake in the Southern Ocean . Along with this goes an underestimated seasonal cycle in CESM, originating from too weak a growing season net flux in CLM4 . MPI-ESM, on the other hand, underestimates atmospheric CO${}_{\mathrm{2}}$ due to weak emissions from LULUC .

The twenty-first century sees substantial emissions from fossil-fuel burning under RCP8.5 (Fig. c). In addition, LULUC is associated with a positive flux into the atmosphere, particularly until around 2050 CE (Table ). After accounting for LULUC (which constitutes a carbon loss for land) the net land sink increases to about 7 Pg C yr${}^{-\mathrm{1}}$ at the end of the twenty-first century (Fig. a). The rate of ocean uptake, on the other hand, peaks around 2070 at about 5 Pg C yr${}^{-\mathrm{1}}$, despite the fact that atmospheric CO${}_{\mathrm{2}}$ continues to rise (Fig. c). This decoupling of the trends in atmospheric CO${}_{\mathrm{2}}$ growth and ocean uptake flux is linked to nonlinearities in the carbon chemistry . The change in dissolved inorganic carbon per unit change in the partial pressure of CO${}_{\mathrm{2}}$ decreases with increasing CO${}_{\mathrm{2}}$, and thus so does the uptake capacity of the ocean. Additionally, differences in the ventilation timescales of the upper and the deep ocean likely play a role. While the surface ocean and the thermocline exchanges carbon on annual-to-multi-decadal timescales with the atmosphere, it takes a century to ventilate the deep ocean, as evidenced by chlorofluorocarbon and radiocarbon data . CESM has a documented low bias in Southern Ocean ventilation due to too shallow mixed layer depths, contributing to the underestimated carbon uptake of the ocean .

The prognostic atmospheric CO${}_{\mathrm{2}}$ increases to 1156 ppm by 2100 CE. This would imply a forcing of 7.6 W m${}^{-\mathrm{2}}$ from CO${}_{\mathrm{2}}$ relative to 850 CE, significantly more than the approximately 6.5 W m${}^{-\mathrm{2}}$ that are imposed by the radiative code calculated according tosee also Fig. c. This propagation of the twentieth-century bias is consistent with the CMIP5 multi-model mean and has motivated attempts to reduce such biases by using observational constraints for ocean ventilation , the tropical land carbon storage sensitivity to temperature variations , and the oceanic and terrestrial carbon fluxes . CESM with CLM4, however, shows very little sensitivity to tropical land carbon, in part due to the inclusion of an interactive nitrogen cycle, which – through enhanced photosynthetic uptake due to nitrogen fertilization – tends to counteract accelerated soil decomposition from warming . Together with the underestimated oceanic uptake, this leads to the roughly 20 % larger airborne fraction in CESM as compared to what is actually prescribed as atmospheric concentration in the radiative code according to the RCP8.5.

Figure  places the current and projected changes within the context of preindustrial variability. Estimated interannual variability prior to 1750 CE is $\pm$0.94 Pg C yr${}^{-\mathrm{1}}$ (1 standard deviation) for the net atmosphere–land and $\pm$0.42 Pg C yr${}^{-\mathrm{1}}$ for the net atmosphere–ocean flux. The much larger interannual variability in land than ocean flux is consistent with independent estimates and results from other models e.g.,. Large volcanic eruptions, as they have occurred in the last millennium, cause anomalously high uptake rates that, for a short period of time, are on par with current uptake rates (Fig. a and b, full range). We estimate when the anthropogenically forced, global mean land and ocean uptake fluxes exceeds the envelope of preindustrial natural variability . As a threshold criteria, it is required that the decadally smoothed uptake fluxes are larger than the upper bound of 2 standard deviations of the annual fluxes prior to 1750 CE. In that case, the simulated global mean land and ocean uptake fluxes have emerged from natural interannual variability by 1947 CE and by 1877 CE, respectively. The prognostic atmospheric CO${}_{\mathrm{2}}$ concentration emerges already in 1755 CE, while the simulated global mean temperature does not emerge until 1966 CE.

Five-year filtered zonal mean anomalies of surface air temperature (SAT), relative to the 850–1849 CE mean from (a) CESM and (b) MPI-ESM. (c) 100-year running-window correlation of zonal mean SAT from CESM and MPI-ESM. A 0.75 Tukey window has been applied to the data before correlation to weaken sharp transitions. Stippling indicates significance (5 % level), taking into account autocorrelation estimated from the entire time period. Panel (d): as (c) but for the correlation of CESM with CCSM4. Panel (e): as (d) but for global mean SAT. Small inset on top shows volcanic and solar forcing of CESM and MPI-ESM. Volcanic forcing of CESM scaled to have the same radiative forcing as MPI-ESM for Pinatubo in 1991 CE. Solar forcing relative to 1850 CE.

[Figure omitted. See PDF]

Model–model coherence

A classical approach to assess the robustness of model results is to rely on the multi-model mean response to a given forcing . However, as there are only very few last millennium simulations with comprehensive Earth system models to date, this approach is not feasible for investigating the decadal-scale climate–carbon-cycle responses to external forcing in the period before 1850 CE. Instead, we estimate periods of forced variability with a 100-year running-window correlation of CESM and MPI-ESM, indicating phasing of the two models. The time series are anomalies from their 850–1849 mean and are smoothed with a 5-year local regression filter before calculating the correlation. Thereby, we focus on the preindustrial period, as the twentieth and twenty-first centuries are dominated by anthropogenic trends, which are nontrivial to remove for a proper correlation analysis. In addition, regression analysis is used.

Temperature

Figure a and b show anomalies of zonal mean annual SAT from CESM and MPI-ESM. In both models the northern high latitudes show the strongest trend, from positive anomalies during the MCA to negative anomalies during the LIA. This is consistent with the current understanding of polar amplification during either warm or cold phases . The twentieth and twenty-first centuries then see the strong anthropogenic warming, although this occurs earlier in CESM due to missing negative forcings from indirect aerosol effects (Sect. ). Superimposed on the preindustrial long-term negative trend are volcanic cooling events. In CESM many of these are global and are able to considerably cool the SH extra-tropics around 60${}^{\circ }$ S, while in MPI-ESM the SH extra-tropics are only weakly affected. These differences are likely related to the Southern Ocean heat uptake rates in the two models arising from under- and overestimation of Southern Ocean mixed layer depths in CESM and MPI-ESM, respectively;. This is evident also in the delayed warming at these latitudes in the twenty-first century in MPI-ESM as compared to CESM. The consistent SH high-latitude positive anomalies before the thirteenth century, on the other hand, appear to be related to a positive phase of the Southern Annular Mode (SAM) in both models (not shown), a behavior common to most PMIP3 models. Note, however, that a recent reconstruction of the SAM finds the models to lack amplitude in their simulated variability, challenging the models' capabilities to represent SAM .

Regression of total solar irradiance (TSI) on surface air temperature (SAT) for the period 850–1850 CE in (a) CESM and (b) MPI-ESM. Time series at each grid point have been 5-year-filtered. Only significant regression coefficients at the 5 % level are shown. The small panel shows zonal means.

[Figure omitted. See PDF]

The phasing on interannual to decadal scales between the two models is largely restricted to periods of volcanic activity and, within those, mainly to land-dominated latitudes (except Antarctica, which shows no forced variability on these timescales; Fig. c). Despite the largest absolute temperature anomalies occurring in the Arctic, the correlations are highest in the subtropics, due to the smaller interannual variability there. Periods of centennial trends, such as the MCA or the Arctic cooling during the Maunder Minimum around 1700 CE, do not show up in the correlation analysis that focuses on 100-year windows, suggesting that multi-decadal low-frequency forcing, such as centennial TSI trends, or internal feedback mechanisms are responsible for the missing correlation. A regression analysis between the 5-year filtered annual TSI and SAT at each grid point (different filter lengths of up to 50 years have been tested as well without the results changing) reveals a clear link between the two quantities at high latitudes. In CESM this link seems to be driven primarily by a displacement of the sea ice edge (Arctic) and Southern Ocean heat uptake (Fig. a). As the sea ice response has not been detected in an earlier model version their Fig. 4, it warrants the question of whether the regression of SAT on TSI might be biased by imprints of volcanoes , even when the time series are filtered, especially in a model like CESM that has a very strong volcanic imprint. Forthcoming simulations with solar-only forcing will be able to answer that question. MPI-ESM, on the other hand, shows a similar polar amplification signal from solar forcing but is not as clearly linked to sea ice (Fig. b). MPI-ESM also displays a stronger land–ocean contrast than CESM see also.

Five-year-filtered zonal mean anomalies of horizontally averaged ocean temperature, relative to 850–1849 CE from (a) CESM and (b) MPI-ESM. Panel (c): 100-year running-window correlation of zonal mean SAT from CESM and MPI-ESM. A 0.75 Tukey window has been applied to the data before correlation to weaken sharp transitions. Stippling indicates significance at the 5 % level, taking into account autocorrelation estimated from the entire time period. Panel (d): 100-year running-window correlation of the Atlantic meridional overturning circulation (AMOC) in CESM and MPI-ESM.

[Figure omitted. See PDF]

In addition to the comparison with MPI-ESM, Fig. d shows results from the correlation analysis between CESM and CCSM4, two simulations that in terms of physics differ only in their applied TSI amplitude and orbital parameters. Not unexpectedly, there are generally more robust signals of forced variability as compared to CESM vs. MPI-ESM (Fig. c), very likely due to the identical physical model components in CESM and CCSM4. Similarly, global mean SAT shows generally stronger phasing between CESM and CCSM4 (Fig. e). However, the latitudinal and temporal pattern of the CESM vs. CCSM4 analysis agrees well with the one arising from CESM vs. MPI-ESM (Fig. c; with exception of the much stronger phasing in CESM and CCSM4 during the volcanic eruptions in the 1450s) and suggests that the physical mechanism behind periods of phasing is robust across the two models.

Five-year filtered zonal mean anomalies of horizontally integrated dissolved inorganic carbon (DIC), relative to 850–1849 CE, from (a) CESM and (b) MPI-ESM. Panel (c): 100-year running-window correlation of zonal mean SAT from CESM and MPI-ESM. A 0.75 Tukey window has been applied to the data before correlation to weaken sharp transitions. Stippling indicates significance at the 5 % level, taking into account autocorrelation estimated from the entire time period.

[Figure omitted. See PDF]

Applied to ocean temperature, the above approach enables us to investigate the penetration depth of a forced signal seen at the surface (Fig. ). Indeed, most of the surface signals also show up as significant correlations down to depths of about 150–200 m, and their timing again suggests volcanic forcing as the origin. Reduced heat loss from the tropical equatorial Pacific together with reduced heat uptake at high latitudes are responsible for ocean cooling after volcanoes (not shown). The Atlantic Meridional Overturning Circulations (AMOC) in the CESM and MPI-ESM shows no significant correlation; however, the highest correlation occurs during the thirteenth century and coincides with a phasing of the upper-ocean temperatures due to strong volcanic forcing (Fig. d). The correlation between CESM and CCSM4 at that time is even higher and points to a significant imprint of the volcanic forcing on ocean circulation . However, during the remainder of the millennium, no phasing of the AMOC is found.

Carbon cycle

We apply the same correlation analysis to zonally integrated land and ocean carbon fluxes from the two models to detect forced variability in the carbon cycle. Compared to SAT hardly any phasing can be found between the models in atmosphere-to-land carbon fluxes (not shown), which is due to its large interannual variability and to distinctly different responses to external forcing in the two models, as will be illustrated in Sect. . Similarly, there is little model phasing in net atmosphere-to-ocean carbon fluxes (not shown). Results become somewhat clearer when considering globally integrated upper-ocean dissolved inorganic carbon (DIC; Fig. ). There appear to exist spurious trends in CESM, likely related to model drift. We repeated the analysis, but with the CESM output detrended in each grid cell by subtracting the CTRL over the corresponding period 850–1372 CE. Due to the shortness of CTRL, we cannot apply this method to the whole simulation. However, these tests showed that the correlation between the two simulations is largely insensitive to the drift in CESM. In Fig. c there are periods of coherent carbon drawdown coinciding with volcanic eruptions around 1450 and 1815 CE in response to temperature-driven solubility changes. Interestingly, MPI-ESM shows a distinct behavior for the strong eruption of 1258 CE, with a prolonged ocean carbon loss after a weak initial uptake. CESM shows a stronger and more sustained carbon uptake, leading to no correlation between the two models for this eruption. The reasons for this discrepancy are discussed in Sect. .

Superposed epoch analysis of the strongest three (top3) and subsequent strongest seven eruptions (top10) of the period 850–1850 CE in (a–e) CESM and (f–j) MPI-ESM for (a, f) global mean surface air temperature, (b, g) global mean precipitation, (c, h) atmospheric carbon given in Pg C on the left $y$ axis and in ppm CO${}_{\mathrm{2}}$ on the right $y$ axis, (d, i) ocean carbon, and (e, j) land carbon. Time series are deseasonalized and calculated as anomalies to the mean of the preceding 5 years. The shading shows the 10–90 % range.

[Figure omitted. See PDF]

Generally, the largest changes in upper-ocean carbon storage occur in response to volcanoes and take place in the tropical Pacific , with other significant changes occurring in the North and South Pacific, the subtropical Atlantic, and the Arctic (Sect. ). Within the tropical oceans, the models show different characteristics: CESM shows a larger variability in DIC than MPI-ESM and, when influenced by anthropogenic emissions in the twentieth and twenty-first centuries, takes up a larger portion of the total ocean carbon uptake than in MPI-ESM (not shown). In MPI-ESM, the Southern Ocean shows stronger variability and larger carbon uptake in the twenty-first century, illustrating the different behavior of the two models in terms of ocean carbon cycle variability and trend magnitude, closely related to the different mixed layer depths in the Southern Ocean region.

Volcanic forcing

To further isolate the response of the climate system and carbon cycle to volcanic eruptions, a superposed epoch analysis is applied to both simulations. Thereby, composite time series for the strongest three (top3) and subsequent strongest seven eruptions (top10), according to optical depth anomaly, over the period 850–1850 CE are calculated for the CESM and MPI-ESM (Fig. ). The time series are calculated as deseasonalized monthly anomalies from the 5 years preceding an eruption.

The physical parameters' global mean surface air temperature and global mean precipitation decrease in both models after volcanic eruptions, although the response of CESM is stronger by roughly a factor of 2–2.5 (Fig. a, b, f, g). Consequently, CESM temperature and precipitation take longer ($\sim$ 15 years) to relax back to pre-eruption values than MPI-ESM ($\sim$ 9 years).

The atmospheric carbon inventory, on the other hand, shows a remarkably different response in the two models. In CESM the atmosphere initially looses about 2–3 Pg C, irrespectively of the eruption strength, with the minimum occurring after about 1–2 years. In the top10 case, values return to normal after about 16 years, while in the top3 case, they tend to return already after about 6 years and overshoot. This overshoot is not straightforward to understand and did not seem to occur in earlier versions of the model . In MPI-ESM the response is a priori more straightforward and slower: in the top10 case the atmosphere looses about 2.5 Pg C, reaches a minimum after 2–4 years, and returns to pre-eruption values after 10–16 years. The top3 case reaches its minimum ($-$6 Pg C) a bit faster, but then takes about 20 years to return to pre-eruption values .

Partitioning these atmospheric carbon changes into land and ocean changes indicates that land is primarily responsible for the differing response behavior of the two models, confirming the findings in the previous section. While in both models, land drives the atmospheric change by taking up carbon initially, it is released back to the atmosphere within about 3 years in CESM but retained in land areas for at least 15 years in MPI-ESM and up to 50 years for the 1258 CE eruption;. In the top3 case of CESM the land starts to even loose carbon after about 5 years, causing the overshoot seen in the atmospheric carbon.

Superposed epoch analysis of the strongest three (top3) and subsequent strongest seven eruptions (top10) for tropical land (25${}^{\circ }$ S to 25${}^{\circ }$ N) in CESM during the period 850–1850 CE. Land carbon inventory changes split up in (a) vegetation, (b) dead biomass (litter and wooden debris), and (c) soil. Furthermore, changes in (d) solar radiation, (e) net primary production (NPP), and (f) loss of carbon through fire. Time series are deseasonalized and calculated as anomalies from the mean of the preceding 5 years. The shading shows the 10–90 % range.

[Figure omitted. See PDF]

A closer look at CESM reveals a distinct response to the top3 and the top10 volcanoes. The response to top3 must be understood as the interaction of a number of processes: the initial global cooling triggers a La Niña-like response and a corresponding cloud and precipitation reduction that is particularly pronounced over tropical land, where large changes in carbon storage also occur (see Fig. a–c for the spatial pattern). Figure  and the following analysis therefore focuses on tropical land. Direct solar radiation decreases, and indirect radiation increases, with a net decrease (Fig. d). These unfavorable conditions cause a reduction in net primary productivity and a strong decrease in vegetation ($-$8 Pg C; Fig. a and e). At the same time, decomposition of dead biomass becomes less efficient due to reduced temperature similar to, e.g., the case in. Despite the simultaneous decrease in net primary production this results in a buildup of dead biomass of about 5 Pg C (Fig. b). Due to the dry conditions and availability of dead biomass, there is increased fire activity, leading to increased carbon loss from land. However, the fire cannot get rid of the large amount of dead biomass immediately (Fig. f). While vegetation decrease and dead biomass buildup balance each other, the soil takes up about 2 Pg C (Fig. c), stores it for at least 16 years, and is therefore responsible for the initial net land uptake seen in Fig. e (see also the left section of Fig. c). After about 2 years, tropical precipitation increases again and puts a halt to the decrease in vegetation (Figs. a and b right). The vegetation does not recover fully for about another 20 years. The dead biomass, on the other hand, gets decomposed entirely within about 15 years and therefore turns land into a carbon source, causing the overshoot in CO${}_{\mathrm{2}}$. In the top10 case, the precipitation and radiation response is about half of the top3 case and so is the vegetation decrease. Consequently, vegetation recovers more quickly. The decomposition of dead biomass, however, takes about the same amount of time as in the top3 case as the decomposition rates are similar for both cases. Hence, land acts as a more sustainable carbon sink in the top10 case. In MPI-ESM it is the soil as well which acts as the main land carbon storage pool, while the vegetation decrease is significantly less than in CESM , leading to the different response behavior of the two land models, particularly striking in the top3 case. Note that there are subtle regional differences between CESM and the earlier version of the carbon-cycle-enabled NCAR model CSM1.4 carbon : tropical Africa sees a reduction in land carbon in CESM, related to a persistent increase in cloud cover and precipitation after volcanoes, while CSM1.4 carbon saw a decrease in precipitation and an increase in land carbon.

Composites of top10 post-volcanic eruption years as anomalies from the preceding 5 years, averaged over (left) the first 2 years starting with the year of the eruption and (right) the following 3 years. (a) Surface air temperature, (b) precipitation, (c) total land carbon, (d) dissolved inorganic carbon (DIC) integrated over the top 200 m. Shading or stippling indicates significance at the 5 % level. Note that for land carbon in an individual grid cell hardly any significant changes are detected due to the large interannual variability.

[Figure omitted. See PDF]

The ocean, on the other hand, shows a qualitatively similar response in CESM and MPI-ESM with an uptake of carbon and a gradual relaxation back to pre-eruption values over 20 or more years. In CESM the radiative cooling leads to increased uptake in the Western Pacific, while in the Eastern Pacific, cooling is less as this region is more controlled by upwelling rather than direct radiative forcing, as suggested by (Fig. d). Two or more years after the volcano, a La Niña-like pattern settles in both surface temperature as well as carbon uptake. Some model differences exist; e.g., in the top3 case of MPI-ESM, the ocean starts to release carbon, compensating for the persistent positive anomaly in the land inventory imposed on the ocean via atmospheric CO${}_{\mathrm{2}}$ concentration;, a feature not present in CESM, in which the land does not store the anomalous carbon for as long. In CESM the tropical oceans appear to be more sensitive to volcanic forcing than to TSI variations. The equatorial Pacific shows the strongest response in DIC to volcanoes (Fig. d), while the response to TSI variations of comparable radiative forcing is up to 1 order of magnitude weaker and confined to higher latitudes (not shown). Overall, it seems, therefore, that the response of the land vegetation governs the overall different responses in the two models.

Global mean changes in response to Pinatubo. (a) Global mean surface air temperature and (b) atmospheric carbon, as (left y axis) Pg C and (right y axis) ppm CO${}_{\mathrm{2}}$ equivalent, both deseasonalized and linearly detrended over 30 years centered on June 1991; temperature observations were corrected for El Niño–Southern Oscillation and other dynamical components ; CO${}_{\mathrm{2}}$ observations were corrected for El Niño–Southern Oscillation and anthropogenic emissions .

[Figure omitted. See PDF]

In an attempt to validate the two models, one is restrained to the well-observed eruption of Pinatubo in 1991 CE, as the CO${}_{\mathrm{2}}$ records from ice cores do not adequately resolve short-term variations induced by volcanoes over the last millennium. Figure  shows the global temperature and atmospheric carbon response to Pinatubo as extracted from observations, CESM, and the three-member ensemble of MPI-ESM. Note that the effects of El Niño–Southern Oscillation and anthropogenic emissions have been removed from the CO${}_{\mathrm{2}}$ observations to obtain a tentative estimate of the actual CO${}_{\mathrm{2}}$ response to the Pinatubo eruption . The initial cooling of about $-$ 0.5 ${}^{\circ }$C and the relaxation back to initial temperatures around 1998 CE is captured well by both models. The MPI-ESM ensemble, however, shows a large and consistent variation around 1995 CE, seemingly related to a phasing of El Niño–Southern Oscillation (ENSO) variability in response to the eruption see also. Further, the magnitude of atmospheric carbon response matches better in CESM, although the overshoot of the observation-based estimate is not captured. CESM's response also falls within the range of the earlier model version . It remains unclear whether this mismatch reflects a model deficiency or is due to uncertainties arising from removing the ENSO signal from the CO${}_{\mathrm{2}}$ observations. However, the mechanisms described above that lead to an atmospheric CO${}_{\mathrm{2}}$ overshoot for large eruptions in CESM offer an opportunity to resolve this discrepancy. Further, the precipitation response (and therewith the cloud and surface shortwave response) to volcanic eruptions is not well constrained due to the small number of observed eruptions . Biases in the representation of these processes can influence a model's carbon cycle response.

Climate–carbon-cycle sensitivity

Due to the absence of large anthropogenic disturbances of the carbon cycle, the last millennium represents a test bed to estimate the climate–carbon-cycle sensitivity, expressed as ppm ${}^{\circ }$C${}^{-\mathrm{1}}$, and can thus potentially help to constrain this quantity e.g.,. Note, however, that there exist important differences between studies in how this sensitivity is calculated and what it implies. Studies using observations and fully coupled simulations estimate the sensitivity from the ratio of changes in CO${}_{\mathrm{2}}$ to changes in temperature. This quantity incorporates feedbacks into the response as the initial sensitivity of the carbon cycle to climate change modifies itself via the climate change that arises from the changed carbon cycle. This is distinct from the climate–carbon-cycle feedback sensitivity $\mathit{\gamma }$, which uses idealized simulations with atmospheric CO${}_{\mathrm{2}}$ held constant, while the climate varies naturally to isolate the feedback parameter . The sensitivity that can be derived from our CESM transient simulation is subtly different again in that the carbon cycle will respond to changes in climate and this response will feed back on the carbon stocks through increased or decreased atmospheric CO${}_{\mathrm{2}}$ concentrations, yet these changes in CO${}_{\mathrm{2}}$ are not allowed to feed back on the climate. Such a sensitivity is expected to be lower than $\mathit{\gamma }$, which we can derive from CTRL.

Here, we estimate the climate–carbon-cycle sensitivity for CESM as follows. We focus on the period before significant LULUC (850–1500 CE) and apply different low-pass filters of 20 to 120 years, taking 5-year increments, to the time series of NH SAT and global CO${}_{\mathrm{2}}$. The filtering aims at minimizing the influence of short-lived forcings such as volcanic eruptions that have a relatively direct impact on temperature and CO${}_{\mathrm{2}}$ (as seen above) and thus may hinder the detection of a low-frequency influence of temperature on CO${}_{\mathrm{2}}$. For each filter length, we determine the highest lag correlation of the two time series, considering lags of up to 100 years. Due to the design of our simulation, we expect NH SAT to lead CO${}_{\mathrm{2}}$, which is confirmed by all lag correlations indicating positive lags for NH SAT (peak of lag correlation at 80.5 $\pm$ 3.4 years). We regress the lagged time series and find a median estimate of 1.3 ppm ${}^{\circ }$C${}^{-\mathrm{1}}$ with a range from 1.0 to 1.8 ppm ${}^{\circ }$C${}^{-\mathrm{1}}$, depending on the filter length. About $-$1 ppm ${}^{\circ }$C${}^{-\mathrm{1}}$ is explained by the land carbon cycle, while the ocean shows smaller sensitivities of about $-$0.4 ppm ${}^{\circ }$C${}^{-\mathrm{1}}$. Note that we use NH SAT in order to be comparable with existing studies . Using global SAT instead of NH SAT can influence the sensitivity estimate, especially for the forced simulation: including the vast ocean area of the SH tends to dampen temperature variability induced by volcanoes and TSI variations. With temperature variability dampened, the sensitivity increases to 1.7 ppm ${}^{\circ }$C${}^{-\mathrm{1}}$ (1.4–2.1).

Temporal dependence of the climate–carbon-cycle sensitivity in CESM. Normalized probability density functions (PDFs) of climate–carbon-cycle sensitivity for 200-year windows overlapping by 50 years (color-filled) for the full period 850–1500 CE (black solid) and for the CTRL (black dashed). The spread of each PDF arises from the range of low-pass filters applied (20 to 120 years).

[Figure omitted. See PDF]

This estimate is barely within the reconstruction-constrained range of 1.7–21.4 ppm ${}^{\circ }$C${}^{-\mathrm{1}}$ and suggests a comparably low sensitivity of the carbon cycle in CESM. This low sensitivity is in agreement with, e.g., . Note that found different sensitivities for the early and late part of the last millennium with the mean for the period 1050–1549 CE being 4.3 ppm ${}^{\circ }$C${}^{-\mathrm{1}}$. Indeed, a strong temporal dependence of the climate–carbon-cycle sensitivity is also found in CESM when looking at individual 200-year windows (Fig. ). The period 1300–1500 CE even shows negative sensitivity, which seems to be related to the different timescales with which SAT and CO${}_{\mathrm{2}}$ relax back to the pre-eruption conditions after perturbations from large volcanic eruptions (Fig. a and c): atmospheric CO${}_{\mathrm{2}}$ decreases after having overshot, while SAT increases after the initial cooling, leading to a negative correlation of the two quantities.

This illustrates the time-variant character of the climate–carbon-cycle sensitivity, which substantially complicates any attempt to constrain it by last millennium data and warrants caution when making inferences from past to future sensitivities. Besides, and found the sensitivity to vary greatly in a coupled model with the timescale and magnitude of volcanic forcing considered. This issue is further highlighted by the larger sensitivity derived for idealized +1 % CO${}_{\mathrm{2}}$ year${}^{-\mathrm{1}}$ simulations with CESM (11.9 ppm ${}^{\circ }$C${}^{-\mathrm{1}}$), for which a dependence on the background state, the scenario, and even the method is reported . Further, it is worth stressing that such sensitivity estimates cannot be extrapolated easily across timescales, as different processes might be at play .

Applying the identical analysis to CTRL reveals other timescales of climate–carbon-cycle feedback, suggesting maximum lags of less than 10 years and a sensitivity of 2.3 (1.4–2.9) ppm ${}^{\circ }$C${}^{-\mathrm{1}}$. Using global SAT instead of NH SAT has no discernible effect (2.3 ppm ${}^{\circ }$C${}^{-\mathrm{1}}$), as the CTRL does not see volcanoes or TSI variations. A later peak in the lag correlation of NH SAT and CO${}_{\mathrm{2}}$ is found at 73.3 $\pm$ 1.1 years in CTRL, i.e., close to where the forced simulation shows its highest lag correlation, but these lag correlations are much weaker ($r$ $\sim$ 0.4 compared to $r$ $\sim$ 0.7 in the forced simulation). This is generally consistent with the finding by that a forced simulation exhibits increased power on lower frequencies compared to a control simulation.

Discussion and conclusions

This study presents a simulation from 850 to 2100 CE with the fully coupled CESM, including the carbon cycle, and provides an overview of the imprint of external forcing on different climate and carbon cycle diagnostics in the simulation. For comparison we draw on a number of PMIP3 simulations, particularly simulations with CCSM4 and MPI-ESM. The evolution of NH SAT during the preindustrial era in CESM is in reasonable agreement with both reconstructions and other models, albeit the uncertainties in reconstructions and forcing still being considerable. Compared to more reliable data in the twentieth century, the anthropogenic warming in CESM is overestimated due to a lack of negative forcing from indirect aerosol effects. In the SH, CESM and most other models do not capture the evolution of the mean SAT as well. The discrepancies could be explained by (i) significant model biases in SH and also interhemispheric SAT variability , (ii) spectral biases in proxies used in the reconstructions , (iii) uncertainties in the external forcing , or (iv) natural internal variability . Unfortunately, these potential explanations are neither exclusive nor independent. Arguments for model bias come from the fact that reconstructed interhemispheric SAT variability lies outside the models' range over 40 % of the time ; but these arguments are weakened by the uncertainty in external forcing. We show here that implementing the same TSI forcing in two different models results in a larger difference in simulated SAT than implementing two different TSI forcings in the same model. Hence, model structural uncertainty remains an issue in determining the role of external forcing over the last millennium.

Albeit beyond the scope of this study, detecting structural and spatial dependencies such as those illustrated here offers an opportunity to reconcile the discrepancies (e.g., regarding SH volcanic signals) between reconstructions and simulations, which might originate from sampling bias, model deficiencies, a combination of these, or the fact that reality may be the one realization that, by chance, is not encompassed by a multi-model ensemble .

Further, we compare simulations with and without orbital forcing and were not able to attribute northern high-latitude SAT trends over the last millennium to orbital forcing. This hampers, if not challenges, the validation of recent findings based on proxy archives that claim a distinct low-frequency orbital component in millennial trends . Instead, the decreasing trend in annual TSI – as opposed to seasonal and regional insolation – together with local feedbacks is able to account for a similar magnitude of trend.

When forced with emissions from LULUC, TSI variations, and volcanic eruptions over the last millennium, both CESM and MPI-ESM do not reproduce atmospheric CO${}_{\mathrm{2}}$ variability, as suggested by ice cores. Notably, the large drop of CO${}_{\mathrm{2}}$ in the seventeenth century is not reproduced, similar to the case in earlier studies . hypothesized that the unique, globally synchronous cooling during the LIA (which might be related to ocean dynamics) can serve as an explanation for this drop. While both CESM and MPI-ESM show a global cooling during the LIA, they develop no apparent phasing of ocean dynamics or carbon uptake and do not show any marked CO${}_{\mathrm{2}}$ reduction around that time, leaving this issue unresolved. The strong volcanic forcing during the thirteenth century, on the other hand, is able to synchronize the AMOC on decadal scales, confirming similar results from the Bergen Climate Model and IPSL-CM5A-LR . With anthropogenic emissions, land and ocean carbon uptake rates emerge from the envelope of natural variability as simulated for the last millennium by about 1947 and 1877 CE, respectively. Atmospheric CO${}_{\mathrm{2}}$ and global temperature emerge by 1755 and 1966 CE, suggesting that changes in carbon-cycle-related variables would be easier to detect than temperature, given sufficient observational data .

We find forced decadal-scale variability in CESM and MPI-ESM in response to major volcanic eruptions in both SAT and upper-ocean temperature, while the response in carbon cycle quantities is less coherent among models see also. Outside volcanically active periods, large parts of the decadal-scale variations cannot be attributed to external forcing, suggesting that internal variability masks external forcing influence. Note, however, that recent work suggests that small volcanic eruptions, which are typically not well-resolved in reconstructions of volcanic activity, exert a significant cumulative effect on global temperature and climate .

Volcanoes trigger a coherent global response in SAT and precipitation that is qualitatively in line with earlier studies on the volcanic influence on climate and the carbon cycle e.g.,. However, the carbon cycle response, in particular on land, shows fundamental model differences in terms of perturbation amplitude and persistence after volcanic eruptions. These differences arise from a differing land vegetation responses in the two models. The extent to which such structural uncertainties matter is illustrated by the large spread in the airborne fraction of CO${}_{\mathrm{2}}$ between these two (and other) models in the twenty-first century see also. In particular, known biases in CESM's carbon uptake in response to anthropogenic emissions in the twentieth and twenty-first centuries lead to a 20 % overestimation of the atmospheric CO${}_{\mathrm{2}}$ concentration and the corresponding prognostic radiative forcing as compared to the prescribed RCP8.5 at year 2100 CE.

The climate–carbon-cycle sensitivity of CESM as estimated from the anthropogenically unperturbed first part of the last millennium is between 1.0 and 2.1 ppm ${}^{\circ }$C${}^{-\mathrm{1}}$, depending on the filtering and the exact time period considered. Generally, the sensitivity of the carbon cycle to temperature variations in CESM is comparably small and reveals a strong component of unforced natural variability. In a transient last-millennium simulation with small temperature variations, the proper detection of a lead–lag relation between temperature and the carbon cycle is complicated by the superposition of perturbations and responses. In addition to the classic climate–carbon-cycle sensitivity experiments e.g.,, it is therefore desirable to conduct step-function-like sensitivity experiments in order to isolate the response of the carbon cycle to a particular external forcing .

Despite the challenges that paleoclimate modeling faces, a number of lessons regarding forcing and structural uncertainties can be learned from these experiments. In order to better understand the role of internal versus externally forced variability – which remains particularly critical for a period of relatively weak external forcing, such as the last millennium – larger simulation ensembles as well as ensembles with decomposed forcing should become standard in paleoclimate modeling. Since these are computationally expensive simulations, this calls for an informed discussion on the optimal usage of computing resources, to which studies like the one here can contribute valuable information. At the same time, uncertainties in forcings and reconstructions need to be further reduced to be able to better validate models in the past with the goal of constraining their future response. Key targets for such constraints are the sensitivity of temperature to solar and volcanic forcing and the climate–carbon-cycle sensitivity.

Acknowledgements

We gratefully acknowledge Axel Timmermann, Bette Otto-Bliesner, Peter Lawrence, and Rosie Fisher for valuable discussions as well as four anonymous reviewers for very helpful comments. We are grateful to the NCAR in Boulder, USA, for providing the code of the CESM, to the World Climate Research Programme's Working Group on Coupled Modelling, which is responsible for CMIP, and to the climate modeling groups for producing and making available their model output. This study is supported by the Swiss National Science Foundation (grant no. 200020 147174), and the European Commission through Seventh Framework Program (FP7) projects CARBOCHANGE (grant no. 264879) and Past4Future (grant no. 243908). J. Mignot has benefited from the support of the French Agence Nationale de la Recherche (HAMOC: ANR 13-BLAN-06-0003). The simulations for this study were performed on a CRAY XT5 and XE6 at the Swiss National Supercomputing Centre (CSCS) in Lugano. Edited by: G. Bala