[1] Tropical explosive volcanism is one of the most important natural factors that significantly impact the climate system and the carbon cycle on annual to multi‐decadal time scales. The three largest explosive eruptions in the last 50 years—Agung, El Chichón, and Pinatubo—occurred in spring/summer in conjunction with El Niño events and left distinct negative signals in the observational temperature and CO 2 records. However, confounding factors such as seasonal variability and El Niño‐Southern Oscillation (ENSO) may obscure the forcing‐response relationship. We determine for the first time the extent to which initial conditions, i.e., season and phase of the ENSO, and internal variability influence the coupled climate and carbon cycle response to volcanic forcing and how this affects estimates of the terrestrial and oceanic carbon sinks. Ensemble simulations with the Earth System Model (Climate System Model 1.4‐carbon) predict that the atmospheric CO2 response is ˜60% larger when a volcanic eruption occurs during El Niño and in winter than during La Niña conditions. Our simulations suggest that the Pinatubo eruption contributed 11 ± 6% to the 25 Pg terrestrial carbon sink inferred over the decade 1990–1999 and −2 ± 1% to the 22 Pg oceanic carbon sink. In contrast to recent claims, trends in the airborne fraction of anthropogenic carbon cannot be detected when accounting for the decadal‐scale influence of explosive volcanism and related uncertainties. Our results highlight the importance of considering the role of natural variability in the carbon cycle for interpretation of observations and for data‐model intercomparison.

1 Introduction [2] The atmospheric CO 2 record measured at Mauna Loa [Keeling et al., 2001] shows considerable variability on interannual to decadal time scales associated with major volcanic eruptions [Jones and Cox, 2001; Frölicher et al., 2011; Brovkin et al., 2010] and ENSO events [Bacastow et al., 1980; Feely et al., 1999]. The origin of this variability and underlying mechanisms, however, are less well understood and quantified than the centennial increase due to human‐induced carbon emissions from fossil fuel and land‐use change. A reason is the lack of or the brevity of spatially and temporarily resolved observational records. In addition, Earth System Modelers have only recently started to employ ensemble simulations for the coupled carbon‐climate system to investigate the role of internal variability [e.g. Frölicher et al., 2009] and related analyses are largely missing in the peer reviewed literature, thereby lagging behind the physical climate modeling community. A quantitative understanding of natural internal and externally forced variability is, however, a prerequisite to detect decadal‐to‐centennial‐scale climate and carbon cycle trends, to properly attribute such trends to anthropogenic and natural forcings from the often short instrumental records, and to facilitate a meaningful evaluation of Earth System Models with observations. [3] Internal variability is the natural variability of the climate system that occurs in the absence of external forcing. The internal variability, or so‐called chaotic component of the climate system, includes processes intrinsic to the atmosphere, the ocean, and the coupled atmosphere‐ocean system [Deser et al., 2012]. An ensemble approach provides a powerful tool for quantifying internal variability. Ensemble simulations are multiple realizations that are forced by identical boundary conditions (e.g., volcanic eruptions) but started from different initial conditions. The central idea is that with a sufficient number of ensemble simulations, with the modes of variability randomized, the ensemble mean can be used to evaluate the deterministic response of the climate system to volcanic eruptions. In a complementary way, the spread of the ensemble distribution can provide an estimate of the envelope of internal variability that is superimposed on the deterministic response to volcanic eruptions. [4] Approximately 40% of the uncertainty related to warming projected over the 21st century stems from the uncertainties associated with the carbon cycle [Meehl et al., 2007]. A key question is how much the land and ocean will continue to absorb anthropogenic CO 2 . A diagnosed positive trend in the airborne fraction, the ratio of CO 2 accumulating in the atmosphere to the total anthropogenic CO 2 emissions (fossil fuel and cement, plus land emissions), has been used by authors of the Global Carbon Project to argue that the sink efficiency of the land and ocean carbon sinks at absorbing excess CO 2 has decreased [Canadell et al., 2007; Le Quéré et al., 2009]. Recent studies, however, challenge the existence of an increasing trend in the airborne fraction [Knorr, 2009], and others question whether it can be interpreted as indicating a decrease in the efficiency of the carbon sinks [Gloor et al., 2010; Sarmiento et al., 2010; Rafelski et al., 2009]. Gloor et al. [2010] and Knorr [2009] show that the noise in the observation data and the uncertainties in land use emission estimates are too large to detect trends in the airborne fraction. In addition, several factors such as non‐exponential growths of carbon emissions, changes in the response of the carbon sinks to rising CO 2 and climate, and natural variability (e.g., El Niño‐Southern Oscillation or volcanic eruptions) influence the variability of the diagnosed airborne fraction [Gloor et al., 2010; Sarmiento et al., 2010; Rafelski et al., 2009]. Thus, the airborne fraction is a relatively poor diagnostic for analyzing changes in the carbon sink efficiency. A quantitative attribution analysis of airborne fraction trends to different mechanisms, however, is missing in the peer‐reviewed literature. [5] The goal of this study is to quantify the role of volcanic forcing in shaping interannual‐to‐decadal variability in atmospheric CO 2 and carbon sinks and to highlight interdependencies with the ENSO, the seasonal cycle, and the role of internal variability. By properly accounting for the influence of volcanic eruptions on atmospheric CO 2 trends, we contribute to the ongoing debate over whether a slowdown of the land and ocean in sequestering anthropogenic carbon is already detectable. [6] Volcanoes contribute to the interannual to multi‐decadal variability of the carbon cycle and the climate system [Robock, 2000]. In the period of available direct atmospheric CO 2 measurement data, after 1958, there were three large tropical volcanic eruptions (Agung in March 1963 at 8°S, El Chichón in April 1982 at 17°N, and Pinatubo in June 1991 at 15°N). All three eruptions occurred in conjunction with a strong El Niño event in boreal spring/summer, and observations show a drop in global mean surface air‐temperature and atmospheric CO 2 . In general, the land biosphere increases carbon uptake during cool climate volcanic events [Jones and Cox, 2001; Frölicher et al., 2011]. However, the extent to which initial conditions and internal variability are important for the carbon cycle response to volcanic eruptions is currently not known. [7] In this study, we use a fully coupled Earth System Model to assess the impact of volcanic eruptions on the physical climate system and the carbon cycle. Sensitivity studies are carried out by imposing a volcanic eruption either in boreal summer or winter, or in association with El Niño or La Niña. We show that the carbon system response critically depends on these initial states, particularly in the first few years after the eruption. We investigate how CO 2 would have evolved over the past 50 years in the absence of eruptions and attribute trends in the airborne fraction to volcanic eruptions and related uncertainties. We show that natural variability related to explosive volcanic eruptions and internal climate modes must be considered for a reliable interpretation of observed trends and to arrive at robust conclusions regarding changes in carbon sinks. [8] The rest of this paper is organized as follows: section 2 presents the model used in this study and the experimental design. The results are described in section 3. A comparison with observationally based carbon fluxes after the Pinatubo eruption and a discussion about the implications for the global carbon budget is presented in section 4. The conclusions are presented in section 5. Details of the calculation of the volcanic signal in the observed temperature and CO 2 records, and the calculation of the airborne fraction are presented in Appendix A and Appendix B.

2 Methods 2.1 Model [9] Simulations were performed with the coupled climate‐carbon cycle model of the National Centre for Atmospheric Research (NCAR) Climate System Model (CSM) 1.4‐carbon [Fung et al., 2005; Doney et al., 2006; Frölicher et al., 2009; Frölicher and Joos, 2010]. The core of the model consists of ocean, atmosphere, land, and sea ice physical components. CSM1.4‐carbon includes a modified version of the Carnegie‐Ames‐Stanford Approach terrestrial biogeochemistry model, termed CASA′, and a derivative of the Ocean Carbon Intercomparison Project 2 oceanic biogeochemistry model. CASA′ incorporates three live vegetation pools and nine dead carbon pools. Autotrophic respiration is not explicitly modeled and is assumed to be 50% of gross primary production. Heterotrophic respiration and carbon flow in litter and soil organic matter pools vary with soil temperature and moisture. Turnover times range from 20 days for the metabolic soil pool to 500 years for the passive pool. The average turnover time scale of less than 5 years of 60% of the soil carbon pools in the NCAR CSM1.4‐carbon model is consistent with flux‐weighted times derived from radiocarbon measurements [Fung et al., 2005; Trumbore, 2000]. [10] Details on strengths and weaknesses of the model can be found in the peer‐reviewed literature [Doney et al., 2006; Schneider et al., 2008; Frölicher et al., 2009; Frölicher and Joos, 2010; Steinacher et al., 2010] where modeled physical and biogeochemical fields are discussed and compared with observations. The carbon cycle response of the model to very large volcanic eruptions is described in Frölicher et al. [2011]. They show that strong reduction in heterotrophic respiration leads to enhanced carbon uptake in tropical regions after volcanic eruptions. This finding is in line with previous studies [Jones and Cox, 2001; Brovkin et al., 2010]. We extend the existing evaluation by analyzing simulated interannual variability in atmospheric CO 2 of the NCAR CSM1.4‐carbon model (Figure 1). Overall, the model captures the observed interannual variability reasonably well. The data over the volcano‐free period 1967–1981 show a linear regression coefficient of 0.38 ppm yr−1 per multivariate ENSO index (MEI) change (r2 = 0.38; green crosses in Figure 1). The linear regression coefficient is slightly smaller when including all data from 1960 to 2010 (all crosses in Figure 1). The model shows a regression coefficient of 0.58 ppm yr−1 per Niño3.4‐index change (r2 = 0.48; grey dots in Figure 1) obtained from a 680 year pre‐industrial control simulation. Large carbon releases during El Niño conditions are simulated in tropical land regions, particularly in the northern part of South America, South East Asia, India, and tropical Africa. This is qualitatively in good agreement with observationally based estimates [e.g., Rödenbeck et al., 2003]. The ocean carbon cycle response to ENSO is marginal in the NCAR CSM1.4‐carbon model. In contrast to observationally based estimates [Feely et al., 1999], largest variability in air‐sea fluxes is simulated in the Southern Ocean and North Atlantic and seems not to be associated with ENSO. This is caused by overly strong iron limitation in the tropical Pacific Ocean in the model and therefore low surface biological uptake [Doney et al., 2006; Schneider et al., 2008; Steinacher et al., 2010]. Feely et al. [1999] showed that the normally prevailing strong outgassing of CO 2 in the equatorial Pacific during an El Niño is reduced due to increased sea surface temperature and decreased supply of nutrient and carbon‐rich water by upwelling. Figure 1 Open in figure viewer PowerPoint Scatterplot of modeled annual mean atmospheric CO 2 growth rates anomalies against sea surface temperature anomalies in the Niño3.4 index region from a 680 year control simulation of the NCAR CSM1.4‐carbon model. Observed annual mean atmospheric CO 2 growth rates anomalies at Mauna Loa against the MEI index are also included. The annual mean atmospheric CO 2 growth rate has been corrected first for the anthropogenic influence. 2.2 Experimental Design [11] Four sets of six‐member ensemble perturbation simulations with identical volcanic forcing targeted to match the aerosol optical depth changes caused by the Pinatubo eruption were performed with the NCAR CSM1.4‐carbon model (Table 1). Six ensemble members are sufficient to represent internal variability in global mean atmospheric CO 2 . A larger number of ensembles do not yield a bigger uncertainty range. The ensemble spread, measured as one standard deviation in simulated atmospheric CO 2 averaged over the first 10 years in the El Niño summer cases, is 0.47 ppm and 0.48 ppm when using four and six ensemble members, respectively. The simulations are carried out by imposing a volcanic eruption either in boreal summer (1 July) or winter (1 January), or in association with El Niño or La Niña. Initial conditions are chosen based upon the phase of the model's own season and ENSO phase from selected points in a quasi‐stable preindustrial control simulation. El Niño (La Niña) initial conditions are randomly chosen when the Niño3.4 index [Trenberth, 1997] exceeds ±0.4°C for at least 9 months in the control simulation. The initial conditions are separated by at least 11 years, which leads to substantially different climate states for each ensemble member. Each ensemble member is perturbed by the similar strength and pattern as the Mount Pinatubo eruption in 1991. The volcanic forcing, which is identical in all simulations, is introduced to the model as prescribed evolution of monthly zonal mean stratospheric aerosol optical depth following Ammann et al. [2003] (Figure 1a in Frölicher et al. [2011]). Thus, potential differences in the atmospheric spread of volcanic aerosols due to season or internal variations of the stratospheric circulation are neglected. The simulated response in climate results from changes in aerosol optical depth and related changes in energy fluxes. The CO 2 response is a consequence of the altered physical climate as direct CO 2 emissions by recent volcanic eruptions are negligible and accounted for less than 0.15% of today's CO 2 emissions [Gerlach, 2011]. Potential effects of tropospheric volcanic aerosols on cloud formation are not included. Table 1. Overview of All Simulations Conducted for This Study Experiment Ensembles Pinatubo forcing Initial Conditions Years El Niño winter 6 1X El Niño and winter 10 6 ‐ El Niño and winter 10 La Niña winter 6 1X La Niña and winter 10 6 ‐ La Niña and winter 10 El Niño summer (standard run) 6 1X El Niño and summer 10 6 ‐ El Niño and summer 10 La Niña summer 6 1X La Niña and summer 10 6 ‐ La Niña and summer 10 [12] The baselines of the sensitivity simulations are delivered by a control simulation for perpetual 1820 pre‐industrial conditions [Doney et al., 2006]. From each ensemble member, the corresponding 10 year period from the control simulation is subtracted. This procedure guarantees that the “pure” volcanic‐induced response of the climate and the carbon cycle is extracted. Then, the ensemble mean and one standard deviation of the ensemble are calculated. The degree of internal variability was measured as the spread (one standard deviation) between the individual members of the ensemble during the integration period.

3 Results [13] The simulations starting with El Niño and boreal summer initial conditions (standard runs) mimic the observed overlap of the Agung, El Chichón, and Pinatubo eruption with El Niño and summer conditions. We will focus our discussion first on these standard runs. The ensemble mean of the standard runs yields a decrease in global mean atmospheric surface air temperature and atmospheric CO 2 of 0.42 ± 0.37°C and 1.64 ± 0.26 ppm, respectively (green lines in Figures 2a and 2b). The temperature decrease in the standard runs is consistent with the observation‐based estimate [Thompson et al., 2009], whereas the global mean CO 2 decrease at the peak underestimates the observation‐based estimate by about 19% (dashed black lines in Figures 2a and 2b). It is very difficult, however, to estimate amplitude, timing, and recovery of the atmospheric CO 2 anomaly from observations as it is based on several assumptions (see Appendix A for more details). Thus the observational‐based estimate must be viewed as a tentative attempt. The underestimation may also be related to the neglect of changes in the fraction of diffusive radiation in the primary productivity algorithm of the NCAR CSM1.4‐carbon model; an increase in the diffusive radiation as seen after volcanic eruptions may boost productivity in the first few years and thus carbon uptake by alleviating light limitation [Gu et al., 2003; Mercado et al., 2009]. Alternatively, the underestimation may point to a low carbon cycle‐climate sensitivity on the multi‐year time scale of the NCAR CSM1.4‐carbon model, which has been identified in transient historical simulations for this model (4.3 ppm/°C) [Frölicher et al., 2011]. This suggests that our simulated CO 2 response to volcanic forcing and internal variability estimates is rather at the lower bound. A larger number of ensembles may not yield a bigger uncertainty range, as the ensemble spread does not significantly decrease when using four instead of six ensemble members. Figure 2 Open in figure viewer PowerPoint 2 , (c) land carbon inventory, (d) ocean carbon inventory, (e) vegetation carbon inventory, and (f) soil carbon inventory between six simulations with and without volcanic eruptions is shown. The eruptions started in El Niño winter (black), La Niña winter (red), El Niño summer (green), and La Niña summer (blue) conditions. Shadings indicate one standard deviation confidence interval of the ensemble simulations. Dashed black lines indicate observation‐derived temperature and atmospheric CO 2 changes after the Mount Pinatubo eruption. Note that the calculation of the observation‐derived atmospheric CO 2 changes must be viewed as a tentative attempt ( Simulated global mean changes after the Pinatubo eruption. Temporal evolution of global ensemble monthly mean differences in (a) atmospheric surface temperature, (b) atmospheric CO, (c) land carbon inventory, (d) ocean carbon inventory, (e) vegetation carbon inventory, and (f) soil carbon inventory between six simulations with and without volcanic eruptions is shown. The eruptions started in El Niño winter (black), La Niña winter (red), El Niño summer (green), and La Niña summer (blue) conditions. Shadings indicate one standard deviation confidence interval of the ensemble simulations. Dashed black lines indicate observation‐derived temperature and atmospheric COchanges after the Mount Pinatubo eruption. Note that the calculation of the observation‐derived atmospheric COchanges must be viewed as a tentative attempt ( Appendix A ). [14] Temperature anomalies in all cases last up to 7 years (Figure 2a)—much longer than the pure volcanic‐induced total surface energy flux changes. The perturbations in CO 2 are multi‐decadal (Figures 2b–2f). Nine years after the perturbation, the land carbon pools are perturbed by 2.6 Gt C in the El Niño summer case (Figure 2c). Thus, the land carbon pool changes account for about 130% of the total −2.0 Gt C atmospheric CO 2 perturbation after 9 years. Most of the carbon is taken up by the soil carbon pools while the vegetation carbon pools show small changes (Figures 2e and 2f). In comparison to the land, the ocean inventory changes are much smaller (Figure 2d). The ocean initially absorbs carbon mainly in response to cooling, but turns into a source after a couple of years in response to the lowered atmospheric CO 2 [Brovkin et al., 2010; Frölicher et al., 2011]. Thus, the ocean compensates somewhat the long‐term decrease in atmospheric CO 2 due to land carbon inventory increases. [15] Next, we analyze the aspects of the terrestrial model that leads to the long‐tail response in the land carbon inventory in the standard case (Figure 3). In general, a small increase in net primary production (NPP) in the tropics (Figure 3a) leads to an enhanced carbon flux from the vegetation carbon pools to the soil pools. The related increase in soil carbon inventory is amplified by the slower decay rate of soil carbon due to temperature and soil moisture changes (Figure 3b). Thus, tropical soil carbon pools with a slowly decomposing time scale of several years remain significantly perturbed by the end of the simulation (Figure 3b). Figure 3 Open in figure viewer PowerPoint Simulated land carbon inventory changes after the Pinatubo eruption in the El Niño summer case. (a) Temporal evolution of zonally integrated ensemble monthly mean differences in vegetation carbon between six simulations with and without volcanic eruptions started in El Niño summer conditions. (b) Same as Figure 3 a, but for soil carbon. [16] Next, we compare all sensitivity simulations. In contrast to the surface energy flux, temperature (Figure 2a), and precipitation, the amplitude and timing of the atmospheric CO 2 response are strongly sensitive to the initial ENSO phase and season (Figure 2b). The CO 2 response is larger (−2.11 ± 0.39 ppm) with a maximum after 2–3 years with El Niño and boreal winter initial conditions, whereas the response is smaller (−1.31 ± 0.48 ppm) with a maximum after 4–5 years with La Niña and boreal summer starting conditions. The impact of initial conditions diminishes with time, as the four different cases converge to similar states after 4–5 years. [17] Differences in atmospheric CO 2 between the cases can be attributed to differences in global mean terrestrial vegetation and soil carbon inventory changes (black and blue lines in Figure 2e and 2f). These changes are in turn driven by changes in NPP and heterotrophic respiration (R h ) in the first four years after the volcanic eruption (Figure 4). Globally, R h is suppressed in all cases, whereas the sign of NPP changes varies between the cases. The differences in the atmospheric CO 2 response between the cases are thus caused by different NPP responses in the NCAR CSM1.4‐carbon model. For example, NPP increases during the first 2 years in the El Niño winter case which, together with a decrease in R h , leads to a large atmospheric CO 2 response. In the La Niña summer case, however, both NPP and R h decrease and the simulated atmospheric CO 2 anomaly is small. Averaged over the first 4 years, El Niño cases tend to have enhanced NPP after a volcanic eruption (El Niño winter: 0.46 Gt C/yr; El Niño summer: 0.17 Gt C/yr; Figures 4a and 4b) and reduced R h (El Niño winter: −0.44 Gt C/yr; El Niño summer: −0.46 Gt C/yr; Figures 4a and 4b), whereas La Niña cases tend to have both reduced NPP (La Niña winter: −0.23 Gt C/yr; La Niña summer: −0.43 Gt C/yr; Figures 4c and 4d) and reduced R h (La Ninã winter: −0.80 Gt C/yr; La Niña summer: −1.00 Gt C/yr; Figures 4a and 4b). Winter cases also tend to have larger NPP than summer cases. Figure 4 Open in figure viewer PowerPoint Temporal evolution of simulated global integrated ensemble monthly mean net primary production and heterotrophic respiration changes between six simulations with and without volcanic eruptions. Pinatubo eruption started in (a) El Niño winter, (b) El Niño summer, (c) La Niña winter, and (d) La Niña summer conditions. Shadings indicate one standard deviation confidence interval of the ensemble simulations. Numbers in all panels show NPP and R h changes in Gt C/yr averaged over the first 4 years after the volcanic eruption. [18] On a regional scale, the differences in carbon uptake between the El Niño winter and La Niña summer case during the first 4 years after a volcanic eruption are located in the tropics and parts of North America and Europe (Figures 5a–5c). Figure 5 Open in figure viewer PowerPoint Simulated regional changes in (a–c) total carbon inventory, (d–f) net primary production, (g–i) soil respiration, (j–l) precipitation, and (m–o) surface temperature averaged over the first 4 years after the Pinatubo eruption. Shown are the differences between the ensemble mean of all 24 ensemble members and the corresponding control simulation (Figures 5 a, 5 d, 5 g, 5 j, and 5 m), the ensemble mean of the anomalies of the six El Niño winter simulations and the ensemble mean of the anomalies of all 24 ensemble members (Figures 5 b, 5 e, 5 h, 5 k, and 5 n), and the ensemble mean of the anomalies of the six La Niña summer cases and the ensemble mean of the anomalies of all 24 ensemble members (Figures 5 c, 5 f, 5 i, 5 l, and 5 o). [19] The NPP (Figures 5d–5f) changes are mainly a result of changes in precipitation (Figures 5j–5l) and soil moisture in the NCAR CSM1.4‐carbon model. Frölicher et al. [2011] show a significant correlation between NPP and soil moisture in these regions after volcanic eruptions (see Figure 7 in Frölicher et al. [2011] for more details). Whereas drying (Figure 5l) and soil moisture decrease lead to a decrease in NPP (Figure 5f) in the La Niña summer case in parts of Africa and northeast South America, precipitation (Figure 5k) and soil moisture increase lead to a net increase in NPP (Figure 5e) in these regions in the El Niño winter cases. R h changes are driven by temperature and soil moisture changes in the NCAR CSM1.4‐carbon model [Figure 7 in Frölicher et al., 2011]. As volcanic‐induced cooling occurs in all cases (Figures 5m–5o), R h decreases in all cases and in most regions (Figures 5g–5i). However, internal variability in NPP and R h is large on a regional scale and complicates a clean attribution of mechanism. Our results show that regional differences in climate response patterns associated with the state of ENSO and the seasonal cycle contribute globally to widespread land carbon cycle responses to volcanic eruptions.

4 Discussion 4.1 Comparison With Observationally Based Carbon Fluxes After the Pinatubo Eruption [20] We compare simulated land and ocean carbon flux anomalies averaged over the first 2 years after the Pinatubo eruption with atmospheric inversion estimates based on atmospheric CO 2 measurements from Bousquet et al. [2000], Rödenbeck et al. [2003], and Baker et al. [2006] (Figure 6). Both data and models show that the land is mainly responsible for the large carbon uptake after Pinatubo. The inversion‐based estimates indicate an anomalous global air‐land carbon flux of 0.88–1.11 Gt C yr−1 and an anomalous air‐sea carbon flux of 0.09–0.25 Gt C yr−1. This is comparable to the simulated carbon flux in the El Niño summer case of 1.38 ± 0.19 Gt C yr−1 for the land and 0.02 ± 0.06 Gt C yr−1 for the ocean, when accounting for the large ensemble spread. The small mismatch may be related to the inclusion of the counterbalancing El Niño signal in the inversion‐based estimates. On regional scales, however, partly conflicting results from the inversions make a quantitative evaluation of the model performance somewhat more difficult. Uncertainties in the inversion results arise from the underconstrained network of CO 2 observations, the choice of the atmospheric transport model, and the sampling network strategy in the inversion studies. Overall, our ensemble‐mean results agree with the multi‐model inversion estimates of Baker et al. [2006] with a large tropical land uptake and relatively small fluxes elsewhere. Individual ensemble members, however, are also consistent with the small tropical land fluxes inferred by Rödenbeck et al. [2003] and Bousquet et al. [2000]. The Bousquet et al. [2000] inversion shows a dipole structure with carbon uptake in North America and a large release in Eurasia, a feature that is not supported by the other two inversions and the model ensemble. Air‐sea carbon flux anomalies are an order of magnitude smaller both in the inversion and in the ensemble simulations. Figure 6 Open in figure viewer PowerPoint Bousquet et al. [ 2000 Rödenbeck et al. [ 2003 Baker et al. [ 2006 Simulated and observationally based carbon fluxes after the Pinatubo eruption are shown. (a) Air‐land carbon fluxes for different regions after the Pinatubo eruption (positive: uptake of carbon through land biosphere; negative: release of carbon). (b) Same as Figure 6 a, but for air‐sea carbon fluxes. Shown are atmospheric inversion estimates (blue:. []; orange: run s90 with 19 sites from. []; and green:. []) and results from this study (white). Results for the NCAR CSM1.4‐carbon model are calculated as the ensemble mean changes from the El Niño summer cases averaged over the first 24 months. Solid and dashed error bars indicate one standard deviation of the six El Niño summer cases and the 24 ensemble cases, respectively. Anomalies for the inversion results are calculated as differences between period June 1991 to May 1993 and June 1990 to May 1991. For each inversion, a 12 month running mean was first calculated to remove the seasonal cycle from the anomaly fluxes. [21] The spread in regional fluxes from the ensemble simulations demonstrates that it is absolutely essential to consider the role of internal climate variability when comparing model with observation‐based evidence to gauge model performance and to avoid premature conclusions. Here, each of the six El Niño summer simulations represents an equally likely evolution of the climate‐carbon cycle system after the Pinatubo eruption. The ensemble of these six El Niño summer simulations as well as the ensemble of the combined 24 El Niño winter, El Niño summer, La Niña winter, and La Niña summer simulations suggests that a wide range of global and regional fluxes responses to a Pinatubo‐like volcanic eruption is plausible. 4.2 Implications for the Global Carbon Budget [22] The multi‐decadal and nonlinear impact of volcanic eruptions on the global carbon cycle has important implications for the detection and attribution that are required in the context of the global carbon budget. Our results from the El Niño summer case suggest that the Pinatubo eruption contributed 2.7 ± 1.6 Pg C (11 ± 6%) to the 25 Pg land sink [Pan et al., 2011] and −0.4 ± 0.2 Pg C (−2 ± 1%) to the 22 Pg ocean carbon sink [Le Quéré et al., 2009] in the 1990s. The internal variability, however, is large for both the land and the ocean. [23] To detect trends in the airborne fraction, the variability induced by volcanoes and ENSO needs to be removed from the atmospheric growth rate. Current studies [Canadell et al., 2007; Knorr, 2009; Le Quéré et al., 2009; Gloor et al., 2010; Sarmiento et al., 2010] assume that the volcanic signal in the atmospheric CO 2 record has time scales of a few years, related to short‐term changes in aerosol optical depth changes, and that the atmospheric CO 2 concentration responds linearly to volcanic eruptions; thus no uncertainties are assigned to the responses. We quantify the uncertainty in the airborne fraction and its trend arising from decadal‐scale nonlinear volcanic CO 2 signals by “correcting” the actual airborne fraction (month by month) for the influence of volcanic eruptions using our ensemble results in a probabilistic Monte Carlo approach (Figure 7). Two hundred sixteen volcano‐free atmospheric CO 2 time series are generated by correcting the actual observed atmospheric CO 2 data at Mauna Loa with the CO 2 time series from our volcanic ensemble simulations (Figures 7a and 7b). The trends in the airborne fraction and the carbon sink fluxes are then calculated following the methodology outlined by Raupach et al. [2008]. A detailed description of the airborne fraction calculation and the Monte Carlo approach is given in Appendix B. Figure 7 Open in figure viewer PowerPoint 2 anomalies resulting from volcanic eruptions, (b) atmospheric CO 2 growth rate as observed and after correcting for the influence of volcanic eruptions, (d) global carbon sink flux to the ocean and land, and (e) airborne fraction over the period 1960–2009 as summarized and described in Table Le Quéré et al. [ 2009 Stocker et al. [ 2011 2 record without correcting for the influence of volcanic eruptions. The multivariate ENSO index (MEI) is also shown in Figure 7b. The three major volcanic eruptions Agung, El Chichón, and Pinatubo are identified by the grey vertical band in all panels. Components of the global carbon budget. Time series of (a) model‐based atmospheric COanomalies resulting from volcanic eruptions, (b) atmospheric COgrowth rate as observed and after correcting for the influence of volcanic eruptions, (d) global carbon sink flux to the ocean and land, and (e) airborne fraction over the period 1960–2009 as summarized and described in Table 2 . The monthly time series in Figures 7 b, 7 d, and 7 e are filtered with a 15 month running mean. (c) Carbon emissions from fossil fuel combustion and cement production (black) and from land use change as given by. [] (brown) and. [] (purple). The shaded areas in Figure 7 a show one standard deviation confidence interval of the ensemble simulations. Thick black lines in Figures 7b, 7 d, and 7 e indicate time series that are calculated from the Mauna Loa COrecord without correcting for the influence of volcanic eruptions. The multivariate ENSO index (MEI) is also shown in Figure 7b. The three major volcanic eruptions Agung, El Chichón, and Pinatubo are identified by the grey vertical band in all panels. Table 2. Observation‐Based Estimates of the Airborne Fraction of Anthropogenic CO 2 Emissions, Its Long‐Term Trend, and Sensitivity of Results to Assumptions Regarding the Influence of Explosive Volcanic Eruptions and CO 2 Emissions From Land Use Airborne Fraction Trend in Airborne Fraction Trend in Global Carbon Sink Flux 1960–1969 1970–1980 1980–1990 1990–1999 2000 to Jun 2009 Jan 1960 to Jun 2009 Jan 1960 to Jun 2009 With volcano correction Standardb 0.46 ± 0.01f 0.45 ± 0.01 0.46 ± 0.02 0.37 ± 0.02 0.45 ± 0.02 −0.17 ± 0.07%/yr 1.62 ± 0.05%/yr Small CO 2 response for Agung and El Chichónc 0.44 ± 0.01 0.46 ± 0.01 0.46 ± 0.01 0.38 ± 0.01 0.46 ± 0.02 −0.05 ± 0.05%/yr 1.55 ± 0.04%/yr Large CO 2 response to volcanoesd 0.47 ± 0.02 0.45 ± 0.02 0.46 ± 0.03 0.38 ± 0.03 0.45 ± 0.02 −0.23 ± 0.09%/yr 1.66 ± 0.06%/yr Different land use emission datae 0.48 ± 0.01 0.47 ± 0.01 0.47 ± 0.02 0.41 ± 0.02 0.47 ± 0.02 −0.12 ± 0.07%/yr 1.53 ± 0.05%/yr Without volcanic correction 0.41 0.47 0.45 0.37 0.48 0.15%/yr 1.40%/yr Le Quéré et al. [ 2009 0.39 ± 0.07 0.45 ± 0.06 0.48 ± 0.05 0.40 ± 0.04 0.45 ± 0.04g 0.3 ± 0.2%/yrg [24] We show that the trend in the “anthropogenic component of the airborne fraction” becomes negative when accounting for volcanoes by using our El Niño summer ensemble simulations (red line in Figure 7 and row 1 in column 6 of Table 2). The airborne fraction trend is slightly positive when no volcanic corrections are applied (row 5 in column 6 of Table 2). The magnitude of the airborne fraction trend depends on the scaling of the smaller Agung and El Chichón eruptions relative to Pinatubo and the correct magnitude of the CO 2 drop after a volcanic eruption. Smaller scaling factors for the Agung and El Chichón eruption resulted in smaller negative airborne fraction trends (row 2 in column 6 of Table 2 and green line in Figure 7). Larger CO 2 scaling for all volcanoes, by assuming that the NCAR CSM1.4‐carbon model underestimates the true response, results in larger negative airborne fraction trends (row 3 in column 6 of Table 2 and blue line in Figure 7). The airborne fraction trends for all cases, however, are not significantly different from zero using a two sided t test (p‐value < 0.05). [25] The uptake rates of the combined land and ocean carbon sink have increased over the last 50 years with an average trend of 1.40%/yr (Figure 7d and row 5 in column 7 of Table 2). When correcting for volcanoes, we show that trends in the uptake rates are even larger (1.62 ± 0.05%/yr in our standard case; row 1 in column 7 of Table 2). [26] It is well known that uncertainties in land use emission data are significant [e.g., Houghton, 2010]. Several papers pointed out that the ability to detect trends in the airborne fraction is heavily impacted by uncertainties in land use emission data [e.g., Le Quéré et al., 2009]. To assess the implications of uncertainty in land use emission data for trends in airborne fraction and carbon sinks, we replaced the land use emission bookkeeping data from Houghton [2003] in our standard case with estimates from Stocker et al. [2011]. We bring the Stocker et al. [2011] land use emission estimates into discussions here, as these estimates include, in addition to the bookkeeping flux considered in the Houghton data, also land use feedback fluxes associated with changes in climate and atmospheric CO 2 [Strassmann et al., 2008]. The Stocker et al. estimates are based on spatially explicit land use maps including cropland, pasture, and built areas, and were not included in earlier studies [e.g., Le Quéré et al., 2009]. The trend in airborne fraction becomes only slightly less negative when using land use change estimates from Stocker et al. [2011] in comparison to the computed airborne fraction using the Houghton [2003] data (row 4 in Table 1 and yellow line in Figure 7). Thus, volcanic eruptions may have a larger impact on airborne fraction trend estimates than uncertainties in land use emission data. [27] Our results are not directly comparable with the results of Le Quéré et al. [2009] (row 6 in Table 2) as they use monthly global volcanic aerosol index data as a measure of volcanically induced aerosol optical depth changes and as a proxy for the impact of volcanic eruption on the global carbon cycle. In earlier studies [Le Quéré et al., 2009; Canadell et al., 2007], a linear assumption between aerosol index data and atmospheric CO 2 growth rate is usually applied. Here, we show that this assumption may not be adequate as we found a much longer impact of volcanic eruption on the carbon cycle than 1–2 years, typically seen when using volcanic aerosol indices, and a large sensitivity on the pre‐existing initial state of the atmospheric CO 2 response to volcanic eruptions.

5 Conclusions [28] In this study, sensitivity simulations are carried out by imposing a volcanic eruption either in the boreal summer or northern winter, or in association with El Niño or La Niña events to study the impact on the physical climate system and the carbon cycle. Implications for the global carbon budget are analyzed. [29] Our analysis shows that the atmospheric CO 2 response to volcanic eruptions depends sensitively on the pre‐existing climate state of the Earth system at the point of the eruption and the evolution of internal climate variability in the first years after the eruption. Furthermore, volcanic eruption produces perturbations to land carbon fluxes that extend for a decade or more. [30] To our knowledge, volcanic sensitivity simulations with a fully coupled Earth System Model have not been used before to assess how atmospheric CO 2 would have evolved over recent decades in the absence of volcanic eruptions. This separation provides a first step toward an attribution analysis of causes of variations in ocean and land carbon sinks and the airborne fraction. Our analysis suggests that the decadal‐scale volcanic effects are underestimated in recent studies describing airborne fraction trends [Le Quéré et al., 2009; Canadell et al., 2007]. Furthermore, we show that the omission of volcanoes is sufficient to explain the observed trend in the airborne fraction [Gloor et al., 2010]. However, analysis of trends in the airborne fraction, as described in this study and elsewhere [Le Quéré et al., 2009; Canadell et al., 2007], is insufficient to make firm conclusions about whether the land and ocean uptake of anthropogenic carbon has slowed down due to recent climate change. As pointed out by Gloor et al. [2010], trends in the airborne fraction are not identical with changes in the ocean and land carbon sinks. [31] Further volcano modeling studies will be needed to constrain our results in the quantitative details. Ensemble simulations are routinely used in the climate modeling community to quantify the amplitudes of the target signal (e.g., change in atmospheric CO 2 after volcanic eruptions) in comparison with the noise (internal variability), but ensemble simulations are not yet standard in the carbon cycle community. Our study may motivate analyses of observations for the coupled atmosphere‐ocean‐land system and model‐observation comparison studies that take full account of the role of natural variability. The quantification of such variability is an essential step in distinguishing natural from human‐induced perturbation of the carbon cycle‐climate system. Future explosive volcanic eruptions have the potential to shed more light on Earth system response relationships and to narrow uncertainties in our quantitative understanding. Volcanic eruptions are of sufficient importance for process understanding so that observing systems on land and ocean should be able to resolve the fingerprints of the Earth system response to volcanic eruptions.

Acknowledgments [51] TLF and JLS are supported by the Carbon Mitigation Initiative project at Princeton University, sponsored by BP, and TLF is also supported by the Swiss National Science Foundation. CCR and FJ are supported by the National Centre of Competence in Research (NCCR) Climate, by the Swiss National Science Foundation, and by the European Project CARBOCHANGE (grant 264879). We thank C. Rödenbeck, P. Bousquet, and B. Stocker for making data from their studies available. Helpful comments from W. Buermann, M. Gloor, and D. Medvigy are acknowledged. The simulations were performed at the Swiss National Supercomputing Centre (CSCS).

Appendix A [32] In this appendix, we provide details about the calculation of the volcanic signal in observed temperature and CO 2 records. Temperature [33] We used three different time series from Thompson et al. 2009 T obs , (ii) the ENSO‐induced variability in global mean monthly mean surface temperature T obs,ENSO , and (iii) the dynamically induced variability (e.g., through North Atlantic Oscillation) in global mean monthly mean surface temperature T obs,dyn (see Figure 4 in Thompson et al. [ 2009 We used three different time series from] (available at [http://www.atmos.colostate.edu/˜davet/ThompsonWallaceJonesKennedy/) : (i) the global mean monthly mean surface temperature, (ii) the ENSO‐induced variability in global mean monthly mean surface temperature, and (iii) the dynamically induced variability (e.g., through North Atlantic Oscillation) in global mean monthly mean surface temperature(see Figure 4 in. [] for details). We calculated the pure volcanic signal in the temperature record as follows:Then, the anthropogenic trend (30 year linear trend centered on the June 1991 Pinatubo eruption) is subtracted from the residual time series. CO 2 [34] We used following three data sets to calculate the pure volcanic CO 2 signal in the CO 2 record: (i) the yearly atmospheric CO 2 data [Keeling et al., 2001, given in column 2 of ftp://ftp.cmdl.noaa.gov/ccg/trends/co2_annmean_mlo.txt ], (ii) the yearly mean fossil fuel and cement emissions (F foss ) as well as the yearly mean land use emissions (F luc ) from Le Quéré et al. [2009] (available at http://lgmacweb.env.uea.ac.uk/lequere/co2/carbon_budget.htm), and (iii) the bimonthly MEI values from Wolter and Timlin [1993,1998] (available at http://www.esrl.noaa.gov/psd/people/klaus.wolter/MEI/table.html). Yearly series of total emission (F E = F foss + F LUC ) and MEI values are constructed. [35] We calculated the pure volcanic CO 2 signal in the observations as follows: [36] The atmospheric growth rate in ppm month−1 is calculated year (t) by year using the equation: [37] We subtracted the contribution of emissions (F E in Pg C yr−1; Figure 7c) from the atmospheric growth rate by using the assumption that 43% (average over period 1959–2006) of the emissions remain airborne: [38] We determined an approximate relationship over the “volcano‐free” period 1967–1981 between the multivariate ENSO index (MEI; Figure 7b) and agr emis_cor by using a linear regression (Figure 1): [39] By subtracting the MEI variations from agr emis_cor , we get an estimate of atmospheric CO 2 without ENSO‐induced variability (agr volc ): [40] CO 2 changes are then calculated as time‐integrated changes of agr volc (t) relative to year 1990. We calculated the pure volcanic COsignal in the observations as follows: [41] We tried to get an uncertainty estimate of agr volc by (i) changing the constant airborne fraction from 43% to 40% and 46%, (ii) using a different land use estimate from Stocker et al. [2011], and (iii) using the entire time period 1960–2010 for the calculation of the linear regression coefficient between MEI and atmospheric CO 2 growth rate. The volcanic signal is marginally sensitive to the used airborne fraction (not shown), but using the Stocker et al. [2011] land use emissions and using the entire time period for linear regression leads to a slightly smaller magnitude of CO 2 peak and shorter recovery time. However, even after the correction for ENSO and emission, the time series agr emis_cor shows large variability associated with volcanoes as well as internal variability making a clear isolation of the volcanic signal difficult to estimate. For example, the atmospheric growth rate in the late 1980s is anomalous low (approximately −0.7 ppm yr−1) and has been associated with an abrupt increase in net land carbon uptake after 1988 [Beaulieu et al., 2012]. Thus, this leads to small anomalies over the period 1991–2000 relative to 1990 and to a “fast recovery” in the observed atmospheric CO 2 signal. Finally, we assumed a linear relationship between MEI and atmospheric CO 2 growth rate as was done in previous analysis [Raupach et al., 2008]. However, the strong El Niño years 1997 and 1998 (Figure 1; years 6 and 7 in Figure 2b) may lead to anomalous high CO 2 fluxes from tropical regions through enhanced biomass burning from recently exposed peat that is very sensitive to drought. The shoot‐off in the high direction in the years 1997 and 1998 may therefore be caused by our linear relationship assumption that does not take into account nonlinearity.

Appendix B [42] In this appendix, we provide detail about the calculation of the airborne fraction. We calculated the trend in the airborne fraction as follows: [43] We generated 216 volcano‐only monthly CO 2 time series over the period January 1960 to June 2010 (CO 2,volc ; Figure 7a). As we have six ensemble members (results from Figure 2b) and three eruptions (Agung, El Chichón, Pinatubo) during this time period, we got n = 216 (6 × 6 × 6) different evolutions of atmospheric CO 2 for the El Niño summer “standard run.” We assumed that a volcanic signal lasts 20 years and extrapolated the time series from our experiments linearly to zero from year 10 to year 20 after the volcanic eruption. We scaled the simulated atmospheric CO 2 anomalies by 0.66 (and by 0.33 for sensitivity tests) for Agung and El Chichón as the observed aerosol optical depth changes for these volcanoes are two thirds of the Pinatubo. [44] We subtracted the n different simulated volcano‐only monthly CO 2 times series over the period January 1960 to June 2010 from the observed smoothed deseasonalized atmospheric CO 2 record of Keeling et al. [2001] (data in column 8 of http://scrippsco2.ucsd.edu/data/in_situ_co2/monthly_mlo.csv) to get n different “volcano‐corrected” CO 2 time series: [45] We calculated the “volcano‐corrected” atmospheric growth rates in Pg C yr−1 (agr volc_cor ; Figure 7b) for each time series n as [46] To filter out the component of the atmospheric growth rate associated with ENSO, we first calculated the total surface atmospheric CO 2 sink flux (F sinks,volc_cor ; Figure 7d) for each time series n: 2 budget occur through the total (land plus ocean) surface‐air exchange flux rather than through the atmospheric growth rate itself [Raupach et al., 2008 F E (t) are fossil fuel and land use emissions in Pg C yr−1 (Figure 7c) and are obtained from Le Quéré (available at http://lgmacweb.env.uea.ac.uk/lequere/co2/carbon_budget.htm). A monthly series of total emission is constructed by linear interpolation of the annual emissions series thereby neglecting a small seasonal cycle in emissions. Uncertainties in the emissions estimates are not taken into account.The MEI‐correlated component is then removed from F sinks,volc_cor for each time series n: μ(n) is the sensitivity of MEI to the F sinks,volc_cor for each realization n. The bimonthly MEI values are obtained from Wolter and Timlin [ 1993 1998 http://www.esrl.noaa.gov/psd/people/klaus.wolter/MEI/table.html) and are related to monthly values as described on the NOAA webpage. [47] The volcano‐ and MEI‐corrected CO 2 airborne fractions (af volc_cor,mei_cor ) are then calculated for each time series n as [48] The linear trend of the airborne fraction (r(af volc_cor,mei_cor )) in units % yr−1 is then calculated for each af volc_cor,mei_cor as volc_cor,mei_cor )′ as the slope and avg(af volc_cor,mei_cor ) as the mean over the time period. The linear trend of the carbon sink fluxes has been calculated accordingly. [49] The statistical significance of the linear trend of each time series n, (af volc_cor,mei_cor )′(n), is calculated using two‐sided Student's t test (p‐value < 0.05). Autocorrelation in the time series has been neglected, as significance levels are low for all cases. In this appendix, we provide detail about the calculation of the airborne fraction. We calculated the trend in the airborne fraction as follows: [50] Additionally, we used the land use emissions estimate from Stocker et al. [2011] (not used in Le Quéré et al. [2009]) to calculate the airborne fraction. The land use emission estimates are first filtered with a spline using a cutoff period of 15 years. Note that we assume constant land emissions for the period 2004‐2009 (Figure 7c).