Sensitivity of a coupled climate-carbon cycle model to large volcanic eruptions during the last millennium

The sensitivity of the climate–biogeochemistry system to volcanic eruptions is investigated using the comprehensive Earth System Model developed at the Max Planck Institute for Meteorology. The model includes an interactive carbon cycle with modules for terrestrial biosphere as well as ocean biogeochemistry. The volcanic forcing is based on a recent reconstruction for the last 1200 yr. An ensemble of five simulations is performed and the averaged response of the system is analysed in particular for the largest eruption of the last millennium in the year 1258. After this eruption, the global annual mean temperature drops by 1 K and recovers slowly during 10 yr. Atmospheric CO2 concentration declines during 4 yr after the eruption by ca. 2 ppmv to its minimum value and then starts to increase towards the pre-eruption level. This CO2 decrease is explained mainly by reduced heterotrophic respiration on land in response to the surface cooling, which leads to increased carbon storage in soils, mostly in tropical and subtropical regions. The ocean acts as a weak carbon sink, which is primarily due to temperature-induced solubility. This sink saturates 2 yr after the eruption, earlier than the land uptake.


Introduction
Volcanic eruptions have a pronounced effect on climate through injection of ample amounts of sulphate aerosols into the stratosphere. Scattering sunlight back to space by the aerosol particles leads to a reduction of solar irradiation approaching the Earth surface and to a consequent cooling of the lower troposphere (e.g. Zielinski, 1998;Robock, 2000). Coupled atmosphereocean models are capable to reproduce the direct climatic effects of volcanic eruptions quite well, although they show a considerable range of dynamic responses (e.g. Stenchikov et al., 2006). The Mount Pinatubo eruption in June 1991, resulted in a maximum global temperature decrease in the lower troposphere by 0.7 K (Santer et al., 2001), while a smaller global mean surface cooling of nearly 0.4 K in the year 1992 (Thompson et al., 2009). This volcanic eruption also had an effect on the atmospheric concentration of the three primary greenhouse gases CO 2, CH 4 and N 2 O (Schauffler and Daniel, 1994 and references therein). In 1992, 1 yr after the 1991 eruption of Mount Pinatubo, the CO 2 growth rate was the lowest for the period 1983-2003, indicating an additional carbon sink of ∼2 PgC in response to the volcanic eruption (Bousquet et al., 2000). This additional carbon sink must be located in the land biosphere or the ocean or a combination thereof. Several biogeochemical mechanisms could be responsible for it. For example, CO 2 solubility in seawater increases with decreasing temperature. The thermodynamic sensitivity of surface water pCO 2 to the temperature change is about 10 ppmv K −1 (e.g. Broecker and Peng, 1982;Archer et al., 2004) but this is valid mainly on time scales longer than decadal. The land response to the climatic change after volcanic eruption is a delicate balance between reduced productivity in response to reduced direct light, cooler and drier conditions and reduced respiration, both autotrophic and heterotrophic (Denman et al., 2007). Usually, the effect of reduced respiration prevails over the lowered productivity and land serves as a net sink of carbon during the cold episodes after volcanic eruptions, although changes in precipitation patterns and nutrient cycling make the land response more uncertain. Another factor that adds to uncertainty is an increase in diffusive radiation after an eruption, which may have created an enhanced terrestrial carbon sink (e.g. Mercado et al., 2009). Interannual variability patterns such as the El Niño-Southern Oscillation phenomenon (ENSO) can cause reduced upwelling of CO 2 rich waters in the equatorial Pacific and thus anomalously low fluxes of CO 2 to the atmosphere during the warm, El Niño phase and stronger than normal CO 2 release during the cold, La Niña phase. As a result, the atmospheric CO 2 can vary with an amplitude on the order of 1-2 ppmv (e.g. D. Keeling, personal communication in Sarmiento, 1993;Winguth et al., 1994), which is similar in size to the observed additional carbon sink. The exact role of biogeochemical mechanisms discussed above is not easy to constrain from observational data and hypotheses to explain it include a fertilization of marine biota due to deposition of volcanic ashes in iron-limited regions (Sarmiento, 1993;Watson, 1997;Duggen et al., 2007).
Coupled climate-carbon cycle models are tools capable to address these issues in a comprehensive manner. Jones and Cox (2001) investigated a response of atmospheric CO 2 to the 1991 Pinatubo eruption using the HadCM3LC model. They showed that simulated CO 2 concentration declines after the eruption, as suggested by observations, although the magnitude of simulated changes was slightly smaller than observed. Enhanced land carbon uptake due to surface cooling after the eruption is the main mechanism that draws down the atmospheric CO 2 concentration in their model. Here, we do not perform an isolated simulation for the Pinatubo eruption but investigate the effect of volcanic eruptions on the carbon cycle using millennium-scale simulations of an Earth System Model (ESM) developed at the Max Planck Institute for Meteorology (MPI-M) (Jungclaus et al., 2010). In these simulations, the model was driven by radiative forcing reconstructed for the years 800-2005. We focus on the period around the very large volcanic eruption in the year 1258, which has the strongest forcing in the reconstructed record (Crowley et al., 2008). Historical chronicles report massive rainfall anomalies, severe crop damage, pestilence and famine across Europe in 1258 (Stothers, 2000). Tree ring data show a cooling of 1.2 K over the Alpine region in the first year after the eruption (Büntgen et al., 2006). However, the decadal-scale cooling signal of the AD 1258 eruption in temperature reconstructions of the last millennium is quite small (Mann et al., 2008). Sensitivity experiments with the MPI-M ESM for the 1258 eruption show that only aerosol particles substantially larger than observed after the Pinatubo eruption yield temperature changes consistent with terrestrial Northern Hemisphere summer temperature reconstruction (Timmreck et al., 2009). The exact geographical location of the 1258 eruption is uncertain. Very large sulphate signals in both Arctic and Antarctic ice cores suggest a low latitude eruption with an initial injection of 260 ± 60 Tg SO 2 (Oppenheimer, 2003). One candidate is the eruption of Quilotoa in the Ecuadorian Andes (0.5 • S, 78.5 • W) (Barberi et al., 1995). Another plausible scenario is that multiple eruptions occurred at northern and southern high latitudes (Kurbatov et al., 2006;Schneider et al., 2009). To explore the uncertainty related to internal variability, an ensemble of five members is performed (MIL) and the averaged response of the system is analysed and compared with the control simulation (CTL).

Methods
Simulations of the last millennium have been conducted applying the MPI-ESM (Jungclaus et al., 2010). The MPI-ESM used in the simulations consists of the general circula-tion models for the atmosphere ECHAM5 (European Centre/ HAMburg model, Roeckner et al., 2003) and for the ocean MPIOM (the Max-Planck-Institute Ocean Model, Marsland et al., 2003). ECHAM5 is run at T31 resolution (3.75 • ) with 19 levels and MPIOM applies a conformal mapping grid with a horizontal resolution ranging from 22 to 350 km. This grid set-up is a low-resolution version of the model used for the CMIP3 (Jungclaus et al., 2006), including scenario simulations for the International Panel of Climate Change and C4MIP simulations (Friedlingstein et al., 2006;Raddatz et al., 2007). Ocean and atmosphere are coupled daily without flux corrections. The carbon cycle model comprises the ocean biogeochemistry module HAMOCC5 (HAMburg Ocean Carbon Cycle, Maier-Reimer et al., 2005;Wetzel et al., 2005) and the land surface scheme JSBACH (Jena Scheme for Biosphere-Atmosphere Coupling in Hamburg, Raddatz et al., 2007). HAMOCC5 simulates inorganic carbon chemistry, phyto-and zooplankton dynamics in dependence of temperature, solar radiation and nutrients as well as formation and dissolution of sediments. JSBACH distinguishes 12 plant functional types and terrestrial organic carbon in five vegetation and soil carbon pools. The carbon fluxes between the ocean, the land and the atmosphere are determined at each coupling time step.
In the 3000-yr CTL simulation, there is no external forcing and the simulated CO 2 evolution is characterized by internal variability. The carbon dioxide concentration is calculated interactively within the model. In the CTL simulation, there is no significant trend in the atmospheric CO 2 concentration (0.001 GtC yr −1 ), which guarantees nearly constant climate conditions. In the land and ocean carbon pools a small compensating trend of 0.007 GtC yr −1 is still present. In the MIL five-member ensemble integrations, the model was driven by superposition of volcanic forcing, solar irradiance, CH 4 and N 2 O concentrations (Meure et al., 2006), land cover changes (Pongratz et al., 2008) and anthropogenic aerosol reconstructions. More details of these forcing reconstructions are discussed by Jungclaus et al. (2010).
Here we summarize the volcanic forcing reconstruction. The radiative volcanic forcing caused by the release of sulphate aerosols into the stratosphere is calculated online in the model using time series of aerosol optical depth (AOD) at 0.55 μm and of the effective radius (R eff ) (Crowley et al., 2008). The data are provided at four equal area latitude bands (30-90 • N, 0-30 • N, 0-30 • S, 30-90 • S) with a time resolution of 10 d. Due to estimates of the timing of the eruption derived from electrical conductivity measurements in Arctic (Clausen et al., 1997;Taylor et al., 1997) and Antarctic ice cores (Traufetter et al., 2004), the date of the largest eruption of the last millennia is set to mid of September 1257 with a maximum of stratospheric sulphate aerosol load in summer 1258.
AOD estimates are based on a correlation between sulphate in Antarctic ice cores and the Stratospheric Aerosol and Gas Experiment (SAGE) and the Stratospheric Aerosol Measurement satellite data for the most recent eruptions (Sato et al., 1993).
R eff growth and decay is based on SAGE satellite observations of the Pinatubo eruption in 1991 (Sato et al., 1993) and what has been attributed to particle sizes of very large eruptions by microphysical simulations (Pinto et al., 1989). In the model, the volcanic aerosol is placed over three stratospheric levels, with a maximum at 50 hPa. Optical parameters for the ECHAM5 radiation scheme are calculated from the time dependent AOD. These are normalized extinction, absorption coefficients and a normalized asymmetry factor assuming a constant SD of 1.8 for the aerosol size distribution. The optical parameters are derived for all spectral bands of the ECHAM5 radiation scheme with 6 short-wave and 16 long-wave wavelength bands between 0.185 and 100 μm. They depend on the assumed aerosol size distribution with an effective radius between 0.02 and 1 μm. Sensitivity experiments for the model response to the 1991 Pinatubo eruption yield a global average temperature drop of 0.4 K consistent with observations (Timmreck et al., 2009).

Results
In response to the largest eruption in the year 1258, the globally annually-averaged near surface air temperature dropped by 1 K (see Fig. 1). The cooling is even more pronounced at the basis of monthly means, where the monthly mean temperature drop between November 1258 and November 1259 is 1.1 K compared to the long-term mean of the control run. After strong cooling over 2 yr, the temperature slowly recovers back to the pre-eruption level. The temperature decrease over land is even more pronounced and reaches a global monthly mean of 1.7 K in October 1259. The spatial distribution of the temperature difference between the year 1259, corresponding to the minimum in global temperature and the year 1257, before eruption, is shown in Fig. 2. In general, the cooling is more pronounced over land than over ocean, except for the equatorial Pacific where the evolving temperature pattern resembles a La Niña type pattern. The temperature decrease over the landmasses in northern high latitudes reaches 2.5 K. The cooling is most significant over the Northern Hemisphere and tropical areas, while regions south of 30 • S show little response. This asymmetric response is mainly explained by the unequal distribution of landmasses between the Northern and Southern Hemisphere. In the vast ocean areas south of 30 • S the large heat capacity of the ocean dampens the temperature signal. Furthermore, the prescribed optical aerosol depth in extratropics after the eruption is higher for the North than for the South motivated by the observed difference in sulphate records in Greenland and Antarctic ice cores (Crowley et al., 2008). This enhanced northern forcing contributes to the asymmetric response as well.
All carbon reservoirs (atmospheric, terrestrial and marine) respond significantly to the volcanic eruptions of the year 1258 as well as to another eruption of the year 1228 (Fig. 3) . Natural variability of all carbon pools as estimated from the control run (shown as grey bars in Fig. 3) is substantial, but the response of the atmospheric CO 2 concentration is significant and clearly beyond the ±1 SD range of the CTL simulation. After the very strong volcanic eruption in the year 1258, the atmospheric carbon pool declines by 4 Gt (corresponding to 2 ppmv, Fig. 3a). The recovery of atmospheric CO 2 concentration to the pre-eruption level takes about 20 yr, much longer than the temperature recovery presented in Fig. 1. This atmospheric CO 2 dynamics (fast decline and slow increase thereafter) is explained by (1) fast response of terrestrial and marine carbon cycle to the abrupt disturbance and (2) longer time scales of biogeochemical mechanisms in comparison with the physical Fig. 1. Anomalies in global mean near surface air temperature (K), averaged over 12 months and over five MIL ensemble members, relative to the averaged temperature of the CTL simulation for the years 1200-1400. The grey bar is for ±1 SD range, calculated over 1000-yr of the CTL simulation. The temperature changes around the 1258 eruption are shown in the inserted box. Arrows and dash lines indicate volcanic eruptions of years 1228 and 1258. processes in the atmosphere and the surface ocean. Marine carbon inventory increases immediately after eruption by 1 GtC and declines to lower than pre-eruption levels within the next 20 yr (Fig. 3b). The terrestrial biosphere takes most of the carbon from the atmosphere in a few years after the eruption (about 3 GtC, Fig. 3c). This carbon is accumulated in the soil carbon pool while the vegetation pool does not change much ( Fig. 3d and e). The land carbon content then stays at the same level over the next 15 yr before declining to the pre-eruption level in about 50 yr. Therefore, there is a difference between the response of terrestrial biosphere, which holds the carbon gains for two decades after the eruption and the ocean response with a faster relaxation towards a level below the pre-disturbance carbon inventory to compensate for the atmospheric CO 2 disturbance. As a result, the magnitude of atmospheric CO 2 concentration response to the volcanic forcing is determined mainly by the land carbon content, while its duration is determined by the marine carbon cycle.
The initial ocean carbon uptake is explained by the cooling effect of the volcanic eruptions. In response to this cooling, solubility of CO 2 in surface waters increases leading to a temporal net gain in the oceanic CO 2 inventory. In the tropical Pacific and to the south of 30 • S, a CO 2 uptake due to the solubility effect is compensated by an enhanced outgassing due to a deepening of the mixed layer. Oceanic CO 2 uptake occurs primarily in extratropical regions and areas of deep-water formation in the Nordic Seas. There is no significant change in the marine biology due to the simulated reduction of solar radiation or changes in surface water temperatures.
The increase in the land carbon content occurs mostly in the tropical regions (Fig. 4a). Regions north of 40 • N are sources of carbon for the atmosphere. A direct effect of the cooling is a reduction of plant productivity and biomass in the boreal region, where the plant growth period is limited by temperature ( Fig. 4b). A reduction in carbon influx into the soil leads to reduced soil carbon content in this region (Fig. 4c), although lower temperature results in a lower CO 2 source due to reduction in the heterotrophic (soil) respiration. In tropical and subtropical regions, where productivity is limited by moisture availability, the effect of decreasing temperature on heterotrophic respiration prevails over the effect of climate change on productivity and biomass. As a result, tropical soils serve as enhanced carbon reservoir after the volcanic eruption (Fig. 4c). The soil carbon model includes two carbon pools, fast and slow and the slow carbon pool has a time scale of several decades in the tropics. This multidecadal time scale governs dynamics of the land carbon response to the volcanic eruption.
The post-eruption carbon uptake occurs despite the common feature of terrestrial and marine carbon cycles to counteract an atmospheric CO 2 decline by decreasing the carbon pools and releasing CO 2 back to the atmosphere. This negative feedback between carbon pools and atmospheric CO 2 , sometimes called 'CO 2 -concentration-carbon feedback' (Gregory et al., 2009), operates through CO 2 fertilization, land CO 2 uptake feedback (e.g. House et al., 2006) and feedback between surface waters pCO 2 and atmospheric CO 2 . As soon as the cooling effect of the volcanic eruption fades away, the ocean tries to re-establish the pre-eruption atmospheric CO 2 concentrations by enhanced sea-air carbon fluxes. The land is able to hold the carbon gain longer than the ocean due to the slower process of decay of excessive soil carbon accumulated during the post-eruption cooling phase.

Summary and discussion
In response to the largest volcanic eruption in the previous millennium in 1258, the simulated atmospheric CO 2 concentration decreases by 2 ppmv. This CO 2 decline is explained mainly by Fig. 3. Dynamics of carbon pools (GtC) in atmosphere (a), ocean (b), total land (c), soil (d) and land vegetation (e); moving average over 36 months. On (a)-(c), the grey bars indicate the ±1 SD range, calculated over 1000-yr of the CTL simulation. This variability range for the soil and vegetation carbon pools, 2422 ± 1.7 and 587 ± 1.2 GtC, respectively, is not shown because of offsets in these pools between the MIL and CTL simulations due to the land cover change forcing. Arrows and dash lines indicate volcanic eruptions of years 1228 and 1258. Note that scaling of carbon storage differs among compartments. reduced heterotrophic respiration of the terrestrial biosphere in response to the surface cooling, which leads to increased soil carbon content, mostly in tropical and subtropical regions. This conclusion is in line with the study by Jones and Cox (2001) who found a strong reduction in both autotrophic and heterotrophic respiration after a volcanic eruption.
The high-resolution ice core CO 2 record from the Law Dome does not reveal a CO 2 drop after the 1258 eruption (Etheridge et al., 1996). However, the strongest CO 2 signal of the order of few ppmv in our simulation is seen on a time scale of less than a decade. This period is still not resolved by the ice core records, where the process of gas diffusion exerts a smoothing of the CO 2 concentration trapped in the air bubbles on the scale of several decades. In this respect, an absence of the strong CO 2 signal after the 1258 eruption in the available ice core records cannot invalidate our results.
In our simulations, land in northern latitudes north of 45 • N serves as a small net source of carbon to the atmosphere. This result contradicts the atmospheric CO 2 inversion by Bousquet et al. (2000) and the modelling study by Lucht et al. (2002) who found an enhanced carbon uptake in the boreal zone after the Pinatubo eruption. However, a later CO 2 inversion study by Baker et al. (2006) attributed the post-Pinatubo carbon sink to the tropical regions and not to the boreal zone. The modelling study by Jones and Cox (2001) also did not reveal an increase in carbon content in the northern high latitudes after the Pinatubo eruption, in line with our simulations. We conclude that while the latitudinal distribution of the land carbon sink is model dependent, an increase of land carbon uptake is the robust finding of both modelling and atmospheric inversion studies.
An asymmetric response of atmospheric CO 2 to the volcanic eruptions-a fast drop followed by a slow recovery-is a prominent feature of the MIL simulations. In general, the atmospheric CO 2 anomaly is about twice as long as the climate anomaly. This asymmetric response is explained by the combined dynamics of terrestrial and marine carbon pools. The magnitude of atmospheric CO 2 concentration response to the volcanic forcing is primarily determined by the land carbon storage while its duration is set up by the marine carbon cycle. The fast drop and slow recovery of the atmospheric CO 2 level is also visible in the post-Pinatubo CO 2 dynamics, although the CO 2 signal is mingled with an ENSO warm phase (Denman et al., 2007).
Our model and the experimental design have some limitations. The terrestrial carbon cycle does not include nitrogen and phosphate cycles, which may alter the land ecosystems response to the climate change. The effect of increased diffuse radiation is neglected albeit the stimulation of NPP by diffuse radiation was probably only enough to compensate for reduced direct light (Angert et al., 2004). Wildfires were not included into our model, while Angert et al. (2004) suggested that reduced biomass burning following eruption could contribute to the post-eruption carbon sink. Effects of volcanic eruption on CH 4 and N 2 O concentrations were not simulated interactively but accounted for implicitly (by using their ice core reconstructions). We did not account for possible deposition of volcanic ashes on the oceanic surface and consequent fertilization of the marine biology. Nonetheless, a qualitative agreement with the post-Pinatubo CO 2 dynamics supports our main result that the atmospheric CO 2 concentration drops after the strong volcanic eruption.
The enhanced land and ocean carbon uptake after volcanic eruptions is a prominent feature in climate-carbon cycle simulations as well as in analyses of observed post-Pinatubo CO 2 dynamics. For this reason large volcanic eruptions are often considered as an analogue for the artificial injection of sulphate aerosols into the stratosphere, which has been suggested as a geo-engineering measure to counteract the warming due to fossil fuel emissions (e.g. Crutzen, 2006;Wigley, 2006;Rasch et al., 2008). However, the volcanic cooling was also associated with a substantial decrease in precipitation over land and a record decrease in runoff and river discharge into the ocean (Trenberth and Dai, 2007) and stratospheric ozone destruction (WMO, 2007). These implications stress out the potential risks of climate engineering (e.g. Brovkin et al., 2009;Hegerl and Solomon, 2009). On longer timescale artificial sulphate injections seem to enhance the oceanic carbon uptake relatively to the experiments without sulphate injections (Brovkin et al., 2009). However, as soon as these massive artificial sulphate injections into the stratosphere would be stopped surface air temperature and the atmospheric CO 2 quickly recover to the pre-injection level. Thus, the geoengineering is not the solution to the global warming problem, in particular because of a danger of catastrophic climate change in case of technological failure as well as ongoing ocean acidification.

Acknowledgments
The work contributes to the MPI-M Millennium and SV projects. S.J.L. acknowledges support by the ENIGMA project from the Max Planck Society. C.T. has partially been support by the German science Foundation DFG grant TI-344/1-1 and Collaborative Research Center 574 'Volatiles and Fluids in Subduction Zones'. The simulations have been enabled with the extensive help of the 'Millennium working group' and carried out on the supercomputing system of the German Climate Computation Centre (DKRZ) in Hamburg. Data processing and storage was provided by the Model and Data group at the MPI-M. We thank Chris Jones and an anonymous reviewer for providing us with helpful and constructive comments.