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 °C−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.

includes an interactive carbon cycle, that continuously covers the last millennium, the historical period, and 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 10 forcing uncertainty for quantities such as hemispheric temperatures or the land and ocean carbon cycle response. Comparing with 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 fail to attribute millennial 15 temperature trends to orbital forcing, as has been suggested recently. The climatecarbon cycle sensitivity in CESM during the last millennium is estimated to be about 1.3 ppm • C −1 . However, the dependence of this sensitivity on the exact time period and scale illustrates the prevailing challenge of deriving robust constrains on this quantity from paleoclimate proxies. In particular, the response of the land carbon cycle to vol-

Introduction
The last about 1000 years constitute the best opportunity prior to the instrumental period to study the transient interaction of external forcing and internal variability in climate, atmospheric CO 2 , and the carbon cycle on interannual to multi-decadal time scales. In fact, the instrumental record is often too short to draw strong conclusions 5 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 testbed for modeling studies. Yet, 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 10 extracting a robust understanding from a period that is characterized by a small signalto-noise ratio in many quantities and for which uncertainties in the external forcing remain (e.g., Wanner et al., 2008;PAGES 2k network, 2013;Schurer et al., 2013). In addition, a process-based quantitative explanation of the reconstructed preindustrial variability in atmospheric CO 2 and carbon fluxes is largely missing. 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 (e.g., Hoffman et al., 2014). Yet, there are only few last millennium simulations including a carbon cycle (e.g., Gerber et al., 2003;Jungclaus et al., 2010;Brovkin et al., 2010;Friedrich et al., 2012). The sensitivity of the carbon cycle to climate has been shown to be mostly positive, i.e., with warming additional CO 2 is released to the atmosphere (Ciais et al., 2013). However, the magnitude of this feedback remains poorly constrained by observations (Frank et al., 2010) and models (e.g. Friedlingstein et al., 2006). In particular, determining the role of the land in past and future carbon cycle variability and trends is still challenging. In both ideal- 15 ized (Doney et al., 2006;Joos et al., 2013) and scenario-guided multi-model studies (Friedlingstein et al., 2014) the 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 Introduction 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 (Hunke and Lipscomb, 2010), including elastic-viscous-plastic dynamics, energy-conserving thermodynamics, and a subgridscale ice thickness distribution. It operates on the same horizontal resolution as POP2.  have an amplitude of 0.2 % between present-day (1961( -1990 and the late Maunder Minimum (1675-1704, which is about twice as large as the 0.1 % used in most PMIP3 simulations: TSI = 2.2635 · (TSI VS09 − TSI VS09 ) + TSI VS09 .

Experimental setup
(1) Figure 1a shows that the TSI used here lies in between the large-amplitude recon-5 struction by Shapiro et al. (2011) and the bulk of small-amplitude reconstructions of the original PMIP3 protocol (Schmidt et al., 2011). Note, that a recent detection and attribution study indicates small amplitude TSI reconstructions to agree better with temperature reconstructions over the last millennium than large amplitude reconstructions (Schurer et al., 2014), in agreement with Ammann et al. (2007). For the twenty-first 10 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 Berger (1978)  order of 30 % to the total carbon emissions from land use (Shevliakova et al., 2009;Houghton, 2010;Stocker et al., 2014). The temporal evolution of long-lived greenhouse gases (GHGs: CO 2 , CH 4 , and N 2 O) is prescribed based on estimates from high-resolution Antarctic ice cores that are joined with measurements at mid-twentieth century (Schmidt et al., 2011, and 5 references therein). While the carbon cycle module of CESM interactively calculates the CO 2 concentration originating from land use changes, fossil fuel emissions (post-1750CE, following Andres et al., 2012, and carbon cycle-climate feedbacks, it is radiatively inactive. Instead, 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 im-10 pact 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 −2 at the year 2100 (Moss 15 et al., 2010). 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 (Landrum et al., 2013). Post-1850 CE, the time-varying aerosol datasets provided by Lamarque et al. (2010Lamarque et al. ( , 2011 are used, whereby CESM only includes a representation of direct aerosol effects. Similarly, nitrogen (NH x and NO y ) input to the ocean is held constant until it starts to be time-varying from 1850 CE onwards, also following Lamarque et al. (2010,2011). Iron fluxes from sediments are held fixed (Moore and Braucher, 2008).

25
Besides to output from current Model Intercomparison Projects, we compare CESM results to those from similar simulation with CCSM4 (Landrum et al., 2013) and IPSL-CM5A-LR (Sicre et al., 2013) 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 ours, apply transient orbital forcing. 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 10 cycle, i.e., the prognostic CO 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 coarser resolved than CESM in both ocean and atmosphere and applies the A1B scenario for the twenty-first century (IPCC, 2000), which corresponds roughly to the current intermediate scenario as compared to the high scenario RCP8.5 used in 15 CESM.
3 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 Cli-20 mate Anomaly (MCA, ∼ 950-1250 CE), a transition into the colder Little Ice Age (LIA, ∼ 1400-1700 CE), followed by the anthropogenically driven warming of the nineteenth and twentieth century (Fig. 2a) 1.23 ± 0.15 • C, while observations report only 0.71 ± 0.13 • C (Cowtan and Way, 2014). 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 (Meehl et al., 2012). The late twentieth century being the warmest period in the NH in the past millennium is consistent with reconstructions (e.g., PAGES 2k network, 2013). 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 10 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 (PAGES 2k network, 2013;Neukom et al., 2014).
The last millennium simulation with CCSM4 shows a largely coherent behavior with 15 CESM in terms of amplitude and decadal variability of NH SAT (850-1850 CE correlation of 5 year filtered annual means r = 0.88, p < 0.001 the last millennium and its climatic impact was likely enhanced through the cumulative effect of three smaller eruptions following shortly after (Gao et al., 2008;Crowley et al., 2008;Lehner et al., 2013). 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 cool-5 ing 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 (Timmreck et al., 2010) and geographic location of past volcanic eruptions (Schneider et al., 2009), 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 currently debates whether the trees' response to volcanic eruptions resembles the true magnitude of the eruption (Mann et al., 2012;Anchukaitis et al., 2012;Tingley et al., 2014). Disagreement among the models exists on the relative amplitude of the MCA, where 15 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 how 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 20 of 2, which in turn complicates a proper detection and attribution of solar forcing during the last millennium (Servonnat et al., 2010;Schurer et al., 2014). 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 Meehl et al., 2012). The simulated mean SAT of the Southern Hemisphere (SH) generally shows a similar evolution as for the NH with the signature of the MCA and LIA superimposed on a weak millennial cooling trend (Fig. 2b). 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 5 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 (Neukom et al., 2014). However, the uncertainties in the early period of the reconstructions prohibits to robustly answer the question whether the models are too global in their response to external forcing. Similar to the NH, the industrial warming in the SH from 1851-1880 to 1981-2010 CE (0.53 ± 0.07 • C) is overestimated by CESM (0.71 ± 0.13 • 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  20 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. 3). While both models experience a negative long term trend in global TSI until about 1850 CE (Fig. 1), the difference arising from the different orbital setup can be seen best in Arctic summer land insola-25 tion (Fig. 3). Hence, Arctic summer land SAT has been proposed as a quantity to be affected by orbital forcing already on time scales of centuries to millennia (Kaufman et al., 2009). However, we find no detectable difference between the two simulations in 363 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | the trend of Arctic summer land SAT (Fig. 3b). 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. 3c shows nonoverlapping 100 and 200 year mean SAT anomalies plotted against the corresponding mean solar insolation. The results from CCSM4 suggest a clear relation of the two 5 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 Kaufman et al. (2009). This is confirmed in new simulations with decomposed forcing, again comparing simulations with fixed and

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 15 land cover changes start to have a large-scale impact on the carbon cycle. The ocean is a carbon source at the beginning and becomes a sink in the second half of the last millennium (not shown). The residual of these fluxes represents changes in the atmospheric reservoir of carbon, illustrated in Fig. 2c by the prognostic CO 2 concentration. The amplitude of the simulated concentration does not resemble the one reconstructed 20 from ice cores (i.e., imposed on the radiative code of CESM), in particular the prominent CO 2 drop in the seventeenth century is not captured by CESM. This raises the question whether the sensitivity of the carbon cycle to external forcing is too weak in CESM, whether the imposed land use changes are too modest (Kaplan et al., 2011;Pongratz et al., 2011), whether major changes in ocean circulation are not captured 25 by models (Neukom et al., 2014), or whether the ice core records are affected by uncertainties due to in-situ production of CO 2 (Tschumi and Stauffer, 2000 et al., 2010). Further, Earth System Models of Intermediate Complexity or vegetation models driven by GCM output do not reproduce the uptake of carbon by either ocean or land needed to explain the reconstructed amplitudes (Stocker et al., 2011;Gerber et al., 2003). The rise in atmospheric CO 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 multimodel mean (Hoffman et al., 2014). From the observational estimates one can diagnose that the discrepancy arises primarily from overestimated carbon release from land (  The twenty-first century sees substantial emissions from fossil fuel burning under RCP 8.5 (Fig. 2c). In addition, LULUC is associated with a positive flux into the atmosphere, particularly until around 2050 CE (not shown). After accounting for LULUC (which constitutes a carbon loss for the land) the net land sink increases to about 7 Pg C yr −1 at the end of the twenty-first century (Fig. 4a). The rate of ocean uptake, on 25 the other hand, peaks around 2070 at about 5 Pg C yr −1 , despite that atmospheric CO 2 continues to rise (Fig. 2c). This decoupling of the trends in atmospheric CO 2 growth and ocean uptake flux is linked to non-linearities in the carbon chemistry. The change in dissolved inorganic carbon per unit change in the partial pressure of CO 2 decreases with increasing CO 2 , and thus the uptake capacity of the ocean. Additionally, differences in the ventilation time scales of the upper and the deep ocean likely play a role. While the surface ocean and the thermocline exchanges carbon on annual-to-multidecadal time scales with the atmosphere, it takes century to ventilate the deep ocean as evidenced by chlorofluorocarbon and radiocarbon data (Key et al., 2004). 5 The prognostic atmospheric CO 2 increases to 1156 ppm by 2100 CE. This would imply a forcing of 7.6 W m −2 from CO 2 relative to 850 CE, significantly more than the approximately 6.5 W m −2 that are imposed by the radiative code (see Fig. 1c). This propagation of the twentieth century bias is consistent with the CMIP5 multi-model mean (Friedlingstein et al., 2014) and has motivated attempts to reduce such biases 10 by using observational constraints for ocean ventilation (Matsumoto et al., 2004), the tropical land carbon storage sensitivity to temperature variations (Cox et al., 2013;Wenzel et al., 2014;Wang et al., 2014) and for the oceanic and terrestrial carbon fluxes (Steinacher et al., 2013). CESM with CLM4, however, shows very little sensitivity in tropical land carbon, in part due to the inclusion of an interactive nitrogen cycle, 15 which -through enhanced photosynthetic uptake due to nitrogen fertilization -tends to counteract accelerated soil decomposition from warming (Lawrence et al., 2012;Wenzel et al., 2014). Together with the underestimated oceanic uptake this leads to the roughly 20 % larger airborne fraction in CESM as compared to the RCP 8.5. Figure 4 puts the current and projected changes into perspective of preindustrial vari-20 ability. Estimated interannual variability prior to 1750 CE is ±0.94 Pg C yr −1 (1 SD) for the net atmosphere-land and ±0.42 Pg C yr −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., Ciais et al., 2013). Large volcanic eruptions, as they have occurred in the last millennium, cause anomalously high uptake 25 rates that for a short period of time are on par with current uptake rates ( Fig. 4a and b, full range). We estimate when the anthropogenically forced, global-mean land and ocean uptake fluxes leave the bound of preindustrial natural variability (Hawkins and Sutton, 2012;Keller et al., 2014). As a threshold criteria, it is required that the decadal- smoothed uptake fluxes are larger than the upper bound of 2 SD of the annual fluxes prior to 1750 CE. Then, the simulated global-mean land and ocean uptake fluxes have emerged from natural interannual variability by 1947 and by 1877 CE, respectively. The prognostic atmospheric CO 2 concentration emerges already in 1755 CE, while the simulated global-mean temperature does not emerge until 1966 CE.

4 Model-model coherence
A classical approach to assess the robustness of model results is to rely on the multimodel mean response to a given forcing (IPCC, 2013). However, as there are only very few last millennium simulations with comprehensive Earth System Models to date, this approach is not feasible to investigate the decadal-scale climate-carbon cycle re-10 sponses 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 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 century are dominated by anthropogenic 15 trends, which are non-trivial to remove for a proper correlation analysis. In addition, regression analysis is used. Figure 5a 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 20 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 (Holland and Bitz, 2003;Lehner et al., 2013). The twentieth and twenty-first century then see the strong anthropogenic warming, although this occurs earlier in CESM due to missing negative forcings from indirect aerosol effects (Sect. 2). Superimposed on Introduction  , 2012;Marsland et al., 2003). 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 10 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 (Abram et al., 2014). 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 15 latitudes, except Antarctica, which shows no forced variability on these time scales (Fig. 5c). 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 20 focuses on 100 year windows, suggesting multi-decadal low-frequency forcing, such as centennial TSI trends, or internal feedback mechanisms to be responsible for the missing correlation. A regression analysis between the 5 year filtered annual TSI and SAT at each gridpoint (different filter lengths of up to 50 years have been tested as well without changing the results) reveals a clear link of the two quantities at high latitudes.

25
In CESM this seems to be driven primarily by a displacement of the sea ice edge (Arctic) and Southern Ocean heat uptake (Fig. 6a). As the sea ice response has not been detected in an earlier model version (Ammann et al., 2007, their Fig. 4), it warrants the questions whether the regression of SAT on TSI might be biased by imprints of volca-  Lehner et al., 2013), even when the timeseries are filtered, especially in a model like CESM that has a very strong volcanic imprint. Forthcoming simulations with solaronly forcing will be able to answer that question. MPI-ESM, on the other hand, shows a similar polar amplification signal from solar forcing, but not as clearly linked to sea ice (Fig. 6b). MPI-ESM also displays a stronger land-ocean contrast than CESM. 5 In addition to the comparison with MPI-ESM, Fig. 5d 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 unexpected, there are generally more robust signals of forced variability as compared to CESM vs. MPI-ESM (Fig. 5c), very likely due to the identical physical model components in CESM and 10 CCSM4. Similarly, global mean SAT shows generally stronger phasing between CESM and CCSM4 (Fig. 5e). However, the latitudinal and temporal pattern of the CESM vs. CCSM4 analysis agrees well with the one arising from CESM vs. MPI-ESM ( Fig. 5c; with exception of the much stronger phasing in CESM and CCSM4 during the volcanic eruptions in the 1450s) and suggest the physical mechanism behind periods of phasing 15 to be robust across the two models.
Applied to ocean temperature, the above approach enables us to investigate the penetration depth of a forced signal seen at the surface (Fig. 7). Indeed, most of the surface signals also show up as significant correlations down to depths of about 150-200 m, whereby their timing suggests again volcanic forcing as the origin. The Atlantic 20 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. 7d). The correlation between CESM and CCSM4 at that time is even higher and points to a significant imprint of the volcanic forcing on ocean circula-25 tion (Otterå et al., 2010;Swingedouw et al., 2013). 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 dis-5 tinctly different responses to external forcing in the two models, as will be illustrated in Sect. 5. 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. 8). While there appear to exist spurious trends at depth in both models, there are periods of coherent carbon 10 draw-down 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. 5. 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. 5). Within the tropical oceans, the models show different characteristics: CESM shows a larger 20 variability in DIC than MPI-ESM and, when influenced by anthropogenic emissions in the twentieth and twenty-first century, 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

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 following strongest seven eruptions (top10), by measure of optical depth anomaly, over the period 850-1850 CE are 5 calculated for the CESM and MPI-ESM (Fig. 9). The time series are calculated as deseasonalized monthly anomalies to 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 2-2.5 ( Fig. 9a, b, f, g). Consequently, CESM tem-10 perature and precipitation take longer (∼ 15 years) to relax back to pre-eruption values than MPI-ESM (∼ 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. 15 In the top10 case values return to normal after about 16 years, while in the top3 case they tend to return already after about six years, and overshoot. This overshoot is not straightforward to understand and did not seem to occur in earlier versions of the model (Frölicher et al., 2011). 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 20 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 (Brovkin et al., 2010).
Partitioning these atmospheric carbon changes into land and ocean changes indicates that the land is primarily responsible for the differing response behavior of the 25 two models, confirming the findings in the previous section. While in both models the land drives the atmospheric change by taking up carbon initially, it is released back to the atmosphere within about 3 years in CESM, but kept in the land for at least 15 years  Brovkin et al., 2010). 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.
A closer look at CESM reveals a distinct response to the top3 and the top10 volcanoes. The response to top3 must be understood as an interplay of a number of 5 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 also large changes in carbon storage occur (see Fig. 11a-c for the spatial pattern). Figure 10 and the following analysis therefore focuses on tropical land. Direct solar radiation decreases, indirect radiation increases, with a net decrease (Fig. 10d). 10 These unfavorable conditions cause a reduction in net primary productivity and a strong decrease of vegetation (−8 Pg C; Fig. 10a and e). At the same time, decomposition of dead biomass becomes less efficient due to reduced temperature (similar to, e.g., Frölicher et al., 2011). Despite the simultaneous decrease in net primary production this results in a build-up of dead biomass of about 5 Pg C (Fig. 10b). Although carbon 15 loss due to fire increases, it cannot get rid of the large amount of dead biomass immediately (Fig. 10f). While vegetation decrease and dead biomass buildup balance each other, the soil takes up about 2 Pg C (Fig. 10c), stores it for at least 16 years, and is therefore responsible for the initial net land uptake seen in Fig 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 Maher et al. (2014) (Fig. 11d). 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, compen-15 sating the persistent positive anomaly in the land inventory (imposed on the ocean via atmospheric CO 2 concentration Brovkin et al., 2010), a feature not present in CESM, in which the land does not store the anomalous carbon 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. 11d), while 20 the response to TSI variations of comparable radiative forcing is up to an 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.
In an attempt to validate the two models, one is restrained to the well-observed 25 eruption of Pinatubo in 1991 CE, as the CO 2 records from ice cores do not adequately resolve short-term variations induced by volcanoes over the last millennium. Figure 12 shows the global temperature and atmospheric carbon response to Pinatubo as extracted from observations, CESM, and the 3-member ensemble of MPI-ESM. Note that around 1995 CE, seemingly related to a phasing of ENSO variability in response to the eruption (see also Zanchettin et al., 2012). Further, the magnitude of atmospheric carbon response matches better in CESM, although the overshoot of the observationbased estimate is not captured. CESM's response also falls within the range of the earlier model version (Frölicher et al., 2013). It remains unclear whether this mismatch reflects a model-deficiency or is due to uncertainties arising from removing the ENSO signal from the CO 2 observations. However, the mechanisms described above that lead to an atmospheric CO 2 overshoot for large eruptions in CESM offer an opportunity for reconciliation of this dispcrepancy. Further, the precipitation response (and therewith the cloud and surface short-wave response) to volcanic eruptions is not well 15 constrained due to the small number of observed eruptions (Trenberth and Dai, 2007). 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 20 millennium represents a testbed to estimate the climate-carbon cycle feedback sensitivity γ (ppm • C −1 ) and can thus potentially help to constrain this quantity (e.g., Wood- SAT and global CO 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 2 (as seen above) and thus may hinder the detection of a low-frequency influence of temperature on CO 2 . For each filter length we determine the highest lag correlation of the two time series, considering lags of up to 100 years. By design of our simulation 5 we expect NH SAT to lead CO 2 , which is confirmed by all lag correlations indicating positive lags for NH SAT (peak of lag correlation at 80.5 ± 3.4 years). We regress the lagged time series and find a median estimate of 1.3 ppm • C −1 with a range from 1.0 to 1.8 ppm • C −1 , depending on the filter length. This is barely within the reconstruction-  (Fig. 13). The period 1300-1500 CE even 15 shows negative γ, which seems to be related to the different time scales with which SAT and CO 2 relax back to the pre-eruption conditions after perturbations from large volcanic eruptions (Fig. 9a and c): atmospheric CO 2 decreases from having overshoot while SAT increases after the initial cooling, leading to a negative correlation of the two quantities.

20
This illustrates the time-variant character of γ, which substantially complicates any attempt to constrain it by last millennium data and warrants caution when making inferences from past to future sensitivities. Besides for idealized +1 %-CO 2 year −1 simulations with CESM (11.9 ppm • C −1 ), for which a dependence on the background state, the scenario, and even the method is reported (Plattner et al., 2008;Arora et al., 2013 Applying the identical analysis to CTRL reveals other time scales of climate-carbon cycle feedback, suggesting maximum lags of less than 10 years and a γ of 2.3 (1.4-2.9) ppm • C −1 . A later peak in the lag correlation of CTRL clusters at 73.3 ± 1.1 years, i.e., close to where the forced simulation shows its highest lag correlation, but these lag correlations are much weaker (r ∼ 0.4 compared to r ∼ 0.7 in the forced simulation).

5
This is generally consistent with the finding by Jungclaus et al. (2010) 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, 10 including carbon cycle, and investigates the imprint of external forcing on different climate and carbon cycle diagnostics. For comparison we draw on a number of PMIP3 simulations, particularly, comparable 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 15 forcing are still considerable. Comparing 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. On 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 (Neukom et al., of the time (Neukom et al., 2014); but these arguments are weakened by the uncertainty in external forcing. We show here that implementing the same TSI forcing in two Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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 illustrated here offers an opportunity to reconcile the discrepancies (e.g., re-5 garding 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 one realization by chance not encompassed by a multi-model ensemble (Deser et al., 2012;Lehner et al., 2012a;Bothe et al., 2013;Neukom et al., 2014). 10 Further, we compare simulations with and without orbital forcing and fail 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 (Kaufman et al., 2009;Esper et al., 2012). Instead, the decreasing trend in annual TSI -as opposed 15 to seasonal and regional insolation -together with local feedbacks are 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 2 variability as suggested by ice cores. Notably, the large drop of CO 2 in the sev-20 enteenth century is not reproduced, similar as in earlier studies (Gerber et al., 2003;Stocker et al., 2011;Jungclaus et al., 2010). Neukom et al. (2014) 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 dy-25 namics or carbon uptake and do not show any marked CO 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 (Otterå et al., 2010;ESDD 6, 351-406, 2015 Climate and carbon cycle dynamics in a CESM simulation from 850-2100 CE  Swingedouw et al., 2013). Under 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 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 observa-5 tional data (Keller et al., 2015). 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. 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 suggest that small volcanic eruptions, which are typically not well-resolved in reconstructions of volcanic activity, exerct a significant cumulative effect on global temperature and climate (Ridley et al., 2014).
Volcanoes trigger a coherent global response in SAT and precipitation that is qualita-15 tively in line with earlier studies on the volcanic influence on climate and carbon cycle (e.g., Jones and Cox, 2001;Brovkin et al., 2010;Frölicher et al., 2011Frölicher et al., , 2013. 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 20 to which such structural uncertainties matter is illustrated by the large spread in the airborne fraction of CO 2 between these two (and other) models in the twenty-first century (see also Friedlingstein et al., 2014). In particular, known biases in CESM's carbon uptake in response to anthropogenic emissions in the twentieth and twenty-first century lead to a 20 % overestimation of the atmospheric CO 2 concentration and the corre-25 sponding 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 about 1.3 ppm • C −1 , with a depen- dency on the filtering and the exact time period considered. Generally, the sensitivity of the carbon cycle to temperature variations in CESM is comparably small (Frank et al., 2010) 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 5 superposition of perturbations and responses. In addition to the classic climate-carbon cycle sensitivity experiments (e.g., Arora et al., 2013) 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 (Gerber et al., 2003). Despite the challenges that paleoclimate modelling 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 vs. externally-forced variability -which remains particularly critical for a period of relatively weak external forcing, such as the last millennium -larger simulation ensembles and ensembles with decomposed forcing should become a standard procedure in paleoclimate modelling. At the same time, 15 uncertainties in forcings and reconstruction 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 constrains remain the sensitivity of temperature to solar and volcanic forcing and the climate-carbon cycle sensitivity. Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Copyright of Earth System Dynamics Discussions is the property of Copernicus
Gesellschaft mbH and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.