Background conditions influence the decadal climate response to strong volcanic eruptions
Abstract
[1] Background conditions have the potential to influence the climate response to strong tropical volcanic eruptions. As a case study, we systematically assess the decadal climate response to the April 1815 Tambora eruption in a set of full-complexity Earth system model simulations. Three 10-member simulation ensembles are evaluated which describe the climate evolution of the early 19th century under (1) full-forcing conditions, (2) volcanic forcing–only conditions, and (3) volcanic forcing–only conditions excluding events preceding the Tambora eruption. The amplitude of the simulated radiative perturbation induced by the Tambora eruption depends only marginally on the background conditions. In contrast, simulated near-surface atmospheric and especially oceanic dynamics evolve significantly differently after the eruption under different background conditions. In particular, large inter-ensemble differences are found in the post-Tambora decadal evolution of oceanic heat transport and sea ice in the North Atlantic/Arctic Ocean. They reveal the existence of multiple response pathways that depend on background conditions. Background conditions are therefore not merely a source of additive noise for post-eruption decadal climate variability but actively influence the mechanisms involved in the post-eruption decadal evolution. Hence, background conditions should appropriately be accounted for in future ensemble-based numerical studies.
Key Points
- The background state affects the decadal climate response to volcanic eruptions
- Background conditions actively influence the climate response mechanisms
- North Atlantic/Arctic oceanic heat transport and sea ice are key factor
1 Introduction
[2] Strong volcanic eruptions (SVEs) can eject large amounts of sulfur into the atmosphere. When the ejected sulfur penetrates into the stratosphere, the consequent temporarily (1–2 years) increased concentrations of sulfate aerosols impose a major transient energy imbalance on the climate system, which is due to aerosol absorption of outgoing infrared radiation and scattering of incoming solar radiation [Robock, 2000; Cole-Dai, 2010; Timmreck, 2012]. Strong near-surface global cooling is a typical direct consequence of the altered radiative energy balance [Robock, 2000; Cole-Dai, 2010].
[3] Short-term climatic impacts of SVEs are, however, not limited to direct radiative effects, as these, particularly for tropical events, trigger dynamical alterations of the atmospheric general circulation that can result in large variations in regional near-surface responses [e.g., Shindell et al., 2004]. Numerical climate simulations produce a considerable range of dynamical short-term responses to volcanic forcing [Stenchikov et al., 2006; Driscoll et al., 2012] likely depending on various aspects of model formulation [e.g., Shindell et al., 2004; Cagnazzo and Manzini, 2009; Spangehl et al., 2010], on the simulated background internal climate variability [e.g., Bengtsson et al., 2006; Thomas et al., 2009], and also on eruption details including magnitude [Timmreck et al., 2011], latitude [Schneider et al., 2009], and season [Toohey et al., 2011]. Important aspects of the short-term dynamical response of the climate to SVEs therefore remain unclear, particularly due to paucity of observed episodes. Moreover, although robust short-term dynamical responses to SVEs are detected in European seasonal temperature and precipitation reconstructions [Fischer et al., 2007; Esper et al., 2013], uncertainty in the background climate state hinders unambiguous discrimination between the contributions of individual internal and external factors to the reconstructed variability [Hegerl et al., 2011].
[4] Considering the oceanic component of the Earth system, the SVE-driven radiative cooling and the associated changes in the atmospheric circulation result in colder and saltier, hence denser, upper ocean conditions over the convective regions, triggering modifications of the oceanic thermohaline circulation [Stenchikov et al., 2009; Otterå et al., 2010; Zhong et al., 2011; Zanchettin et al., 2012a; Ortega et al., 2012; Miller et al., 2012]. This allows the cold signal to penetrate deep into the ocean and to significantly affect the global ocean temperature and heat content for several decades or even longer [Delworth et al., 2005; Church et al., 2005; Gleckler et al., 2006; see also Sévellec and Fedorov, 2012], thereby protracting the recovery of the climate system from volcanic perturbations. Numerical simulations suggest that, similar to short-term dynamical responses, also long-term responses are only partially constrained by the initial radiative anomaly. Background climate conditions are important sources of variety and uncertainty in the simulated near-decadal oceanic response to SVEs [Mignot et al., 2011; Miller et al., 2012; Wang et al., 2012; Zanchettin et al., 2012a, 2012b]. Here “background climate conditions” refers to both the initial climate state and the presence and strength of additional external forcings. Background conditions may impact the efficiency of the various processes and feedback mechanisms involved in the post-eruption climate evolution. Accordingly, whereas systematic long-term responses are simulated at the regional scale, e.g., protracted cooling in the Arctic [Schneider et al., 2009; Zanchettin et al., 2012a], the distinct traits of the long-term climate response to individual eruptions may be potentially elusive, particularly in the noisy atmospheric and near-surface components [Zhong et al., 2011; Zanchettin et al., 2013].
[5] In this study, we systematically show the dependency of simulated volcanically forced decadal climate variability on background constraints. The eruption of Mount Tambora in April 1815 serves as a case study, as it is the largest known historical SVE [Oppenheimer, 2003; Cole-Dai, 2010], which followed the 1809 VEI-6 tropical eruption of unknown location [Crowley et al., 2008; Cole-Dai et al., 2009] and occurred during the ~1790–1830 Dalton solar minimum. Due to the combined effects of the 1809 and Tambora eruptions, the decade 1810–1819 was the coldest during the past 500 years or longer in the Northern Hemisphere and the tropics [Cole-Dai et al., 2009]. By investigating the simulated post-Tambora climate evolutions in a set of simulation ensembles using the identical aerosol forcing from the Tambora eruption but different imposed boundary and initial climate state conditions, we aim at answering the following questions: Can background conditions significantly affect the simulated dynamical decadal climate response to a strong tropical volcanic eruption? And if so, do background conditions simply add noise to the forced response or, rather, diversify it deterministically? We focus on oceanic processes and sea ice dynamics in the North Atlantic/Arctic sector, which has been identified as a critical region for long-term climate responses to volcanic forcing [Stenchikov et al., 2009; Otterå et al., 2010; Zhong et al., 2011; Zanchettin et al., 2012a; Ortega et al., 2012; Miller et al., 2012]. We discuss the implications of our findings for the dynamical interpretation of reconstructed climate variability and for decadal climate predictability. We further explore ensemble forced responses associated to a succession of perturbations induced on a simple, idealized nonlinear system.
2 Data and Methods
[6] We use the COSMOS-Mil version of the Earth system model developed at the Max-Planck-Institute for Meteorology (MPI-ESM) [Jungclaus et al., 2010], which is based on the coupled atmosphere-ocean ECHAM5-MPIOM general circulation models. The atmosphere model ECHAM5 [Roeckner et al., 2006] is run at T31 resolution (3.75° × 3.75° spatial resolution), with 19 vertical levels and a model top at 10 hPa. The ocean model MPIOM [Marsland et al., 2003; Jungclaus et al., 2006] uses a conformal mapping grid with a horizontal grid spacing of 3.0° and 40 vertical levels and includes a dynamic/thermodynamic sea ice model. The full carbon cycle is considered with modules for terrestrial biosphere (JSBACH) as well as ocean biogeochemistry (HAMOCC). Notwithstanding its rather coarse resolution and an only partial description of stratospheric dynamics, the COSMOS-Mil version of MPI-ESM represents an optimal choice for studies requiring long simulations in a largely populated ensemble thanks to its cheap computational requirements. In addition, COSMOS-Mil simulations for the last millennium [Jungclaus et al., 2010] comply with the “Paleoclimate Modelling Intercomparison Project Phase III” requirements [Schmidt et al., 2011]. Based on these simulations, Zanchettin et al. [2012a] discuss the typical interannual and decadal climate responses to tropical strong volcanic eruptions in the context of the MPI-ESM dynamical framework. COSMOS-Mil simulations for the last millennium compare well with reconstructions of multidecadal North Atlantic sea surface temperature (SST) variability under periods dominated by external forcings [Zanchettin et al., 2012b]. Their consistency with available temperature reconstructions is limited to some specific regions and periods of the last millennium [Bothe et al., 2012]. In particular, COSMOS-Mil simulations are climatologically and probabilistically consistent with reconstructed annual central European mean temperatures [Dobrovolný et al., 2010] for the last ~500 years [Bothe et al., 2012].
[7] Three 10-member ensembles are created to investigate the simulated climate evolution around the Tambora eruption under different background conditions. The three ensembles differ in the applied external forcing, and they entail 10 all-forcing (AF) simulations, including all natural (solar, volcanic, orbital) and anthropogenic (aerosols and greenhouse gas, land-cover changes) forcing, 10 volcanic forcing–only (VO1) simulations, and 10 additional volcanic forcing–only (VO2) simulations, where the forcing by the 1809 eruption is omitted. The reconstructions of natural and anthropogenic forcing are those described in Jungclaus et al. [2010]. In particular, volcanic forcing is calculated online, using 10-day average time series covering the last millennium of aerosol optical depth (AOD) at 0.55 µm and of effective particle radius (Reff) for four equal-area latitudinal bands [Crowley et al., 2008; see also Crowley and Unterman, 2012]. Volcanic aerosols are placed over three stratospheric levels, with a maximum at 50 hPa. Optical parameters for the 22 wavelength bands in the range from 0.185 to 100 µm of the ECHAM5 radiation scheme are calculated from the time-dependent AOD and Reff, assuming a constant standard deviation of 1.8 µm for the aerosol size distribution. Sensitivity experiments for the model response to the 1991 Pinatubo eruption yield a global average temperature drop comparable with observations [Timmreck et al., 2009].
[8] Individual simulations within each ensemble differ in the initial climate state at the time of the Tambora eruption. AF simulations started in year 1751 from initial conditions taken from two of the full-forcing COSMOS-Mil experiments described in Jungclaus et al. [2010]. The ensemble is then created by inducing infinitesimal perturbations in the atmospheric vertical diffusivity during the first integration years. In VO2, the AOD and Reff input values are set equal to their background values before the Tambora eruption. Initial climate states for simulations in the VO1 and VO2 ensembles are taken from a multimillennial control run describing an unperturbed preindustrial climate [Timmreck et al., 2010]. The global deep-ocean (below 2200 m) potential temperature in the control run is affected by a ~0.005 K/century trend. While this may generate biases in the upper ocean [Sévellec and Fedorov, 2012], these are included in our estimation of the natural variability range. We expect therefore the trend in the deep ocean not to affect the validity of our general inferences. Boundary conditions for the control run are set constant at year 800 values; therefore, the initial states for VO1 and VO2 simulations are generally warmer than those for AF simulations. In both VO1 and VO2, the 1815 Tambora eruption is initialized in April of the 7th integration year. In VO1, the 1809 eruption is initialized in January of the first integration year, in line with the tentative date for this eruption set on February 1809 based on ice core timing data [Cole-Dai et al., 2009]. Figure 1 provides an overview on the initial climate states of all the individual simulations. It highlights within- and inter-ensemble differences in the sampled climate state. The AF ensemble, in particular, samples from overall colder—compared to VO1 and VO2, and the control run—SST conditions over the North Atlantic (see Appendix A for a definition of ocean basins), while individual AF simulations sample the warmest SST conditions in the tropical Pacific (Figure 1a). The AF ensemble further samples from overall stronger states of the North Atlantic subpolar gyre and weaker states of the Atlantic meridional overturning circulation (AMOC), which is described here as the maximum of the mass transport stream function in the Atlantic near 30°N and below 1000 m depth (Figure 1b). We also remark that the different phasing between individual processes adds to their individual pre-eruption states in generating the complexity of simulated responses to external forcing [Zanchettin et al., 2012a].
[9] The ocean heat transports HT (with reference 0°C) are calculated as in Zanchettin et al. [2012a] based on the equation HT = ΣzΣxvTcpρdxdz, where v (in m/s) is the velocity component perpendicular to the respective section (Denmark Strait, Iceland-Scotland Ridge, different latitudes in the North Atlantic), T (in °C) is temperature, cp (in J/(kg °C)) is heat capacity, ρ (in kg/m3) is density, dz (in m) represents the vertical integral, and dx (in m) represents the horizontal (along the section) integral. For the calculation of the AMOC component of the heat transport at the different latitudes in the North Atlantic, vmTm is used in the above-mentioned equation, with m indicating the zonal mean. For the calculation of the gyre component of the heat transport at the different latitudes in the North Atlantic, v′T′ is used in the above-mentioned equation, with ′ indicating the deviation from the zonal mean.
[10] A simple assessment of the competing effects of the sea ice–albedo and water vapor feedbacks and their dependence on temperature is performed. The sea ice–albedo feedback is diagnosed through the difference between post-eruption and pre-eruption values of the albedo of ice. The water vapor feedback is diagnosed through the ratio between post-eruption and pre-eruption values of LW0/LWs↑, where LW0 is the global average top of atmosphere (TOA) thermal radiation and LWs↑ is the global average upward thermal radiation at the surface.
[11] Oceanic signals are diagnosed by investigating deseasonalized and low pass–filtered (61-months backward running mean) values. Atmospheric signals are diagnosed through superposed epoch analysis [e.g., Fischer et al., 2007] on deseasonalized and low pass–filtered (3-month backward running mean) values. In the superposed epoch analysis, post-eruption anomalies are evaluated with respect to the pre-eruption climatology, defined as the average state over the decade 1799–1808 to avoid including spurious signals from the 1809 eruption.
[12] Statistical significance is estimated based on the likelihood of a random occurrence of the signal [e.g., Graf and Zanchettin, 2012]. Specifically, ensemble-average signals of VO1, VO2, and AF are compared to those from a set of 500 random ensembles created by sampling the 10 eruption years from the whole length of the control run. If one ensemble member is investigated individually, significance is similarly estimated based on one eruption year. The distribution yielded by the years/random sequences is used as a basis to evaluate the confidence level associated to a chance occurrence of the signal. In the estimation of the significance, the autocorrelation is therefore preserved. Additionally, the Mann-Whitney U (MWU) test [Steel and Torrie, 1960] is used to assess the significance of the difference between each ensemble and the other two through the integration time. The MWU tests the null hypothesis that data in two groups are (independent) samples from continuous distributions with equal medians against the alternative that they are not. It is a nonparametric rank-based alternative to the independent t test.
[13] Compatibility between simulated and observed/reconstructed climates of the early 19th century is preliminary assessed by comparing the AF ensemble with the observation-based Berkley Earth Surface Temperature (BEST) data [Rohde et al., 2012] and the proxy-based Northern Hemisphere (NH) temperature reconstruction ensemble by Frank et al. [2010] (full raw calibration ensemble). We further compare the simulated winter (DJF) 500 hPa geopotential height (Z500) anomalies for the first two post-eruption winters with respect to the period 1799–1808 with corresponding reconstructions available for the North Atlantic/European region [Casty et al., 2007]. The proposed assessment is meant to delineate the range of uncertainty to circumscribe the representativeness of our general inferences for the specific case of the 1815 Tambora eruption. It is not meant as an evaluation of the realism of simulated climate responses to volcanic forcing in MPI-ESM [to this regard, see Timmreck et al., 2009; Zanchettin et al., 2012a, 2013] nor as an in-depth validation of the ensemble [e.g., Bothe et al., 2012]. Neither is it meant to challenge considerations about a possible discrepancy between simulated and reconstructed magnitude of volcanically induced cooling [Mann et al., 2012; Anchukaitis et al., 2012; Brohan et al., 2012; Esper et al., 2013].
3 Results
[14] In the following, first, the compatibility of simulated and observed/reconstructed temperature and atmospheric responses to the Tambora eruption is analyzed. Then, the imposed forcing is assessed (4), and key dynamical atmospheric (4) and oceanic (5) features are presented, highlighting the dualism between the inter-ensemble consistency of the imposed forcing and the inter-ensemble differences in the forced decadal climate pathways. 6 investigates potential drawbacks due to the (unavoidably) limited sampling of the assessed responses. 7 substantiates our findings by illustrating the state-dependent clustering of post-perturbation evolutions in a simple, idealized nonlinear system.
3.1 Simulation-Reconstruction Compatibility
[15] Figure 2 compares the simulated (AF ensemble only) and observed/reconstructed evolutions of the annual land-only global surface (2 m) air temperature (GST) (Figure 2a) and of the annual NH surface air temperature (Figure 2b) around the 1815 Tambora eruption. Note that no winter temperature reconstruction exists for the full NH. The two temperature data sets provide consistent pictures of simulation-observation/reconstruction compatibility, both indicating that simulations feature, in their ensemble-average value, a stronger and longer-lasting post-Tambora cooling compared to observations/reconstructions. Due to the large ensemble spread, the coldest post-Tambora GST anomalies for individual simulations exceed the coldest anomaly in the observation-based estimate by a factor between 1.3 and 2 (Figure 2a, note the large reconstruction's uncertainty). The relative magnitude of post-Tambora and post-1809 coldest anomalies inversely depends on the magnitude of post-1809 cooling (compare individual simulations in Figure 2a). Such dependency is consistently found in the VO1 ensemble (not shown). This means that the 1809 eruption is a major factor affecting the GST response to the Tambora eruption, hence indirectly the scaling between simulated and observed/reconstructed post-Tambora GST anomalies. Notably, contrasting with the post-Tambora evolutions, the range of simulated GST evolutions is more constrained and closer to the observation-based estimate after the 1809 eruption.
[16] The proxy-based annual NH temperature reconstructions provide a more constrained estimate than observations do for GST (Figure 2b). Reconstructions depict a decadal-scale negative anomaly rather than large, short-lived, post-eruption excursions. Indeed, even though consisting of annually resolved data, the raw calibration ensemble by Frank et al. [2010] includes time-filtered series and is therefore most representative of decadal and longer fluctuations. Additionally, we note that large-scale temperature signals extracted from tree ring chronologies, as employed in this reconstruction ensemble, are largely restricted to decadal and longer fluctuations [Cook et al., 2004; D'Arrigo et al., 2006]. Accordingly, the simulated temperature exceeds the reconstructions’ uncertainty range after both the 1809 and Tambora eruptions if no smoothing is applied to it (dashed black line in Figure 2b). Decadally smoothing the simulated data results in the ensemble-average trajectory as well as most individual simulation trajectories lying within the reconstructions’ uncertainty range most of the time. Some smoothed simulated trajectories even closely trace the reconstructed ensemble-average trajectory.
[17] In summary, simulations and observations/reconstructions provide overall compatible representations of global and hemispheric temperature evolution around the 1815 Tambora eruption if ensembles are treated as probabilistic ensembles (see Bothe et al. [2012] for a probabilistic assessment of last-millennium MPI-ESM simulations). In the following analyses, we highlight differences in the ensemble-average responses and treat within-ensemble variability as uncertainty. Accordingly, we compare ensembles following the common truth-centered approach, which highlights differences in the expectation values. Concentrating on ensemble-average evolutions and patterns does not affect the generality of our inferences since inter-ensemble-average differences reflect both initial-state and boundary-forcing driven contributions to background condition–related variability. At the same time, this simplifies the presentation of our results.
[18] Geopotential height at 500 hPa (Z500) anomalies depicting a positive phase of the North Atlantic Oscillation (NAO) during the first and second winters after tropical SVEs are a typical feature of reconstructed atmospheric variability [Fischer et al., 2007; Hegerl et al., 2011]. The reproduction of such a feature still represents a challenge for dynamical climate models due to a deficient representation of stratospheric processes and of stratosphere-troposphere interactions as well as to internally generated variability masking the forced signal [Stenchikov et al., 2006; Zanchettin et al., 2012a; Driscoll et al., 2012]. Figure 3 compares reconstructed and simulated (AF ensemble) post-eruption DJF Z500 anomalies for the 1809 (top panels) and Tambora (bottom panels) eruptions. Reconstructions [Casty et al., 2007] show a post-eruption dipolar anomalous pattern for both the 1809 and Tambora eruptions, with centers displaced compared to a typical positive NAO-like pattern in the first event (Figures 3a and 3d). Anomalies in the centers of actions of large-scale atmospheric circulations are largely smeared out in the ensemble-average simulated responses to both eruptions (Figures 3b and 3e). Consequently, the lowering of geopotential height anomalies over the tropics emerges as the only robust forced signature at the hemispheric scale considering both eruptions. This is consistent with the expected thermodynamical response, i.e., the thickening of the lower tropospheric layers driven by the volcanically forced surface cooling, which is initially stronger within the tropical latitudinal band [Zanchettin et al., 2012a]. Nonetheless, ensemble-mean post-Tambora eruption anomalies over the extra-tropical North Pacific consistently describe a strong deepening of the Aleutian low (Figure 3e). Within the variety of post-eruption anomalous patterns from individual simulations, some are found to be consistent with the reconstructed patterns, particularly for the Tambora event (Figures 3c and 3f). The limited significance of Z500 anomalies in individual simulations indicates that they are compatible with internally generated variability. Therefore, reconstructed post-eruption Z500 anomalies for the early 19th century are compatible with the range of simulated behaviors expressed by the AF ensemble. Their attribution as predominantly forced short-term dynamical responses is, however, questionable.
3.2 Imposed Forcing, and Atmospheric and Near-Surface Response
[19] Figure 4 summarizes the imposed forcing in the three ensembles, as it is diagnosed through the globally averaged full-sky TOA shortwave, longwave, and net radiative flux anomalies. The post-Tambora eruption evolutions feature the typical [e.g., Timmreck et al., 2009, 2010; Zanchettin et al., 2012a] twofold radiative effect of the stratospheric aerosols: (1) a reflection to space of the solar radiation, here accounting for a ~11 W m−2 peak reduction of the downward shortwave flux (Figure 4a), and (2) an entrapment of thermal radiation, here resulting in a ~6 W m−2 maximum net decrease in the outgoing longwave radiation (Figure 4b). The net TOA radiative imbalance accounts for a peak reduction of ~4.5 W m−2 about one year after the beginning of the eruption (Figure 4c). Inter-ensemble differences in the forcing are sporadic during the first four post-eruption years and are more appreciable in the thermal radiation. Also, at the peak forcing, the ~0.5 W m−2 inter-ensemble-average spread in the net radiative anomaly is nonsignificant. Whereas the shortwave radiation reverts back to pre-eruption values after approximately 5 years, the decay time of the thermal radiation anomaly is longer, leading to a temporary positive (i.e., strengthened downward) net radiation. Such behavior emerges also in the clear-sky radiation (not shown), suggesting that clouds do not play a role in this simulated feature.
[20] Figure 5 summarizes the simulated evolution of global atmospheric and near-surface anomalies around the 1815 Tambora eruption for the three ensembles. Plotted variables are global averages for the surface net radiative fluxes, the latent heat flux, GST, and the global cloud cover. The post-Tambora anomaly of total net surface energy flux (solar and thermal radiation, sensible and latent heat) tracks the TOA evolution, accounting for an all-ensemble-average peak value of ~−4.5 W m−2 followed by a weaker positive anomaly (not shown). Focusing on the surface radiative flux anomalies (Figure 5a), the different ensembles yield very similar maximum imbalances, particularly for the 1809 event (VO1 and AF ensembles only). Differently from TOA radiative and surface energy fluxes, the post-eruption evolution of surface radiative flux does not entail positive anomalies. As shown in Figure 5b, the positive anomaly is related to a positive (downward) latent heat anomaly, likely a reduced latent heat release from the cooler ocean surface.
[21] The simulated GST response to the 1809 eruption is slightly larger in AF compared to VO1 (Figure 5c). In AF, GST anomalies remain significantly colder than the range of natural variability until the Tambora eruption. GST anomalies are, however, not significantly different in the three ensembles at the time of the Tambora eruption. The post-Tambora cold anomaly is significantly smaller and shorter-lived in VO2 (i.e., when the forcing from the 1809 eruption is omitted) compared to VO1 and AF. On the other hand, the peak cold excursion in AF occurs later and is then significantly larger than in VO1 and VO2. After the maximum anomaly, interannual fluctuations become apparent in all ensemble-average evolutions. The fluctuations have different phases and amplitudes in the different ensembles, and it is not clear whether the eruption acts as a pacemaker for simulated interannual climate variability. The global cloud cover (Figure 5d) confirms the twofold effect of the eruption, with a significant reduction in the first two post-eruption years suddenly followed by a nearly same-amplitude significant increase. The negative cloud-cover anomaly is strongest in VO1. The noisier character of ensemble-average global cloud-cover anomalies compared to other variables highlights the role of the global hydrological cycle as a constraint for the simulated volcanically induced climate anomaly. In short, the volcanic eruption primarily affects the evaporative part of the hydrological cycle. The hydrological cycle feeds back on the forced radiative imbalance by modulating the amount of incoming solar radiation reaching the surface (by means of cloud cover) and the radiative heating within the atmosphere (by means of convection) (see Andrews et al. [2009] for a discussion on radiative forcing and the hydrological cycle).
[22] Figure 6 further illustrates how the volcanically induced changes in the water vapor and sea ice–albedo feedbacks depend, in the first three post-Tambora eruption years (January–December), on the background global surface temperature. Generally, a strengthening of the feedback would correspond to increased values in the difference for the sea ice–albedo feedback and to decreased values in the ratio for the water vapor feedback. For both feedbacks, the assessment is complicated by the presence of stratospheric volcanic aerosols during the first post-eruption years. The albedo of ice strengthens in post-eruption years, further showing a positive dependence on the pre-eruption global surface temperature (r = 0.72, p < 0.0001 for year 1816) (Figure 6a). Such dependence of albedo of ice on temperature contrasts with the general relationship found in an unperturbed climate, though the latter is characterized by large uncertainty (Figure 6a, inset). The ratio between surface and planetary albedo, the latter defined by the ratio between the TOA upward short-wave radiation and the incoming solar radiation, is indeed reduced by ~8% in the first post-eruption year compared to pre-eruption conditions (not shown). The effect of sea ice on the planetary albedo is therefore overwhelmed by the scattering of incoming solar radiation by volcanic aerosols. The water vapor feedback strengthens during the first two post-eruption years (Figure 6b). Individual simulations indicate that such strengthening depends on the pre-eruption global surface temperature (r = 0.75, p < 0.0001 for year 1816) (Figure 6b). Specifically, colder climates are expected to yield a comparatively strengthened water vapor feedback in post-eruption years and vice versa. This is counterintuitive, since, as also shown by MPI-ESM under unperturbed conditions (Figure 6b, inset), the Clausius-Clapeyron equation indicates a clear positive dependency of atmospheric water vapor on the surface temperature (~0.6%/K). The thermodynamical dependence may therefore be overwhelmed by dynamically induced changes (e.g., changes in the large-scale circulation). The different constraints imposed on the two feedbacks complicate our inferences about the short-term forced response of the global climate system under different background states. Overall, the direct forcing dominates the short-term response and significant inter-ensemble differences in post-Tambora global average atmospheric and near-surface conditions emerge only sporadically.
[23] Figure 7 illustrates the simulated regional responses of northern hemispheric annual near-surface air temperature immediately after the Tambora eruption (1st–3rd post-eruption years, top panels) and around one decade after the Tambora eruption (10th–12th post-eruption years, bottom panels) for the three ensembles. Consistently, among the ensembles, a heterogeneous and extensively significant cooling characterizes the immediate response (Figures 7a–7c). Large extensive negative anomalies (in the range from −1.5 to −1.25 K) are robustly detected over landmasses, particularly over North America, Western Europe, and subtropical continental Asia. Temperature anomalies over the Bering Strait/Alaska region are generally nonsignificant or positive (in VO2). Together with the extremely cold conditions over North America, this corresponds to a large-scale atmospheric circulation change identifiable as a positive phase of the Pacific North American (PNA) pattern. In VO1 and AF, but not in VO2, the hemispheric anomaly pattern entails strong cooling over the Barents Sea/Nordic Seas tracing the advance/retreat region of the sea-ice front.
[24] Approximately one decade after the eruption, most regions have reverted back to the natural variability range (Figures 7d and 7e). Consistent for all ensembles, the tropical and the eastern North Atlantic remain significantly colder than in the pre-eruption period, though only locally in VO2. Over landmasses, the cold anomaly is more extensive in AF, covering Western North America, Western Europe, Siberia, and Eastern Asia. The Sahara and the Arabian Peninsula tend to remain colder in VO1 and VO2, but not in AF. So, on the one hand, all ensemble post-eruption evolutions indicate initial strong near-surface cooling followed by a recovery to pre-eruption conditions in most regions one decade after the eruption. On the other hand, inter-ensemble differences arise in both the magnitude and the duration of local anomalies, which trace back to the efficiency of the various processes and feedback mechanisms involved in the post-eruption climate evolution.
3.3 Oceanic Response
[25] Figures 8-10 summarize the simulated evolution of the global ocean and of the (extra-tropical) North Atlantic/Arctic Ocean states around the Tambora eruption for the three ensembles (see Appendix A for definition of ocean basins). On the global scale, upper ocean temperatures remain outside the natural variability range and evolve near parallel throughout the period encompassing the Tambora eruption and the following eruptions of Babuyan Claro and Cosiguina in the early 1830s (Figure 8a). Pre-eruption conditions govern the post-eruption upper ocean temperature differences between the three ensembles. Post-Tambora global ocean heat losses are comparable for all ensembles, though they are notably stronger in VO2 in their peak values (Figure 8b). The higher upper ocean temperatures in VO2 imply a larger sea surface atmosphere temperature difference supporting a more efficient energy transfer at the air-sea boundary. The negative flux anomaly characterizing the 5th–15th post-eruption years is, within certain subperiods, significantly larger in AF than in VO1 and VO2. Global heat flux anomalies in the different ensembles seem to converge again at the onset of the 1830s eruptions (Figure 8b).
[26] Similar to the global ocean, simulated North Atlantic upper ocean temperatures drop significantly outside the natural variability range during the first post-eruption decade (Figure 8c). Peak negative anomalies are comparable in VO1 and VO2, where temperatures in the second post-eruption decade recover to within the natural variability range. In AF, the peak negative anomaly roughly doubles that in VO1 and VO2, and temperatures remain outside the natural variability range until the early 1830s eruptions. The North Atlantic variability within the ensembles is generally larger than that of global values, particularly for AF in the second post-eruption decade. This reflects both the variety of local dynamical responses initiated by the volcanic radiative perturbation and the noisier character of local and basin-scale processes. Major inter-ensemble differences emerge for the ocean heat releases in the North Atlantic Ocean (Figure 8d). The ensembles differ in the peak value of the anomaly, with significantly larger values in AF, but not in its persistence. Different from the global ocean, the North Atlantic does not experience significant below-normal heat loss during the 5th–15th post-eruption years.
[27] A ~0.7 Sv post-eruption strengthening of the AMOC emerges as a robust inter-ensemble feature (Figure 9a). The 1830s eruptions, occurring within the time frame of the typical duration of post-eruption AMOC fluctuations in this model [Zanchettin et al., 2012a], keep AMOC values outside the natural variability range [see also Zanchettin et al., 2010]. A significant increase in Northern Hemisphere sea-ice cover is consistently diagnosed in all ensembles after the eruption (Figure 9c), but the average fluctuations in the different ensembles differ in both magnitude and duration. The similarity of pre-1809 conditions and the differences in post-Tambora evolutions in VO1 and AF suggest that boundary conditions including variations in the total solar irradiance (in this case, the Dalton solar minimum) influence feedback mechanisms involving sea ice dynamics.
[28] The meridional ocean heat transport in the North Atlantic can be separated into two contributions associated to the AMOC and to the gyre circulation, respectively (see 2). In all ensembles, the AMOC strengthening allows for a significant, though temporary, increase in the northward heat transport about one decade after the eruption, especially from the subtropical region to the midlatitudes (Figure 10a). Initially, however, the AMOC acts toward a reduction of northward heat transport especially in the midlatitude to subpolar region. The North Atlantic subpolar gyre is significantly stronger in AF than in VO1 and VO2 during the first post-Tambora decade (not shown, but compare initial states in Figure 1b), suggesting a stronger recirculation of temperature anomalies within the North Atlantic basin. This is confirmed by the significant reduction of gyre-related northward heat transport detected in AF during the first five post-eruption years (Figure 10b). In AF, the combined AMOC and gyre effects produce a strong initial reduction of northward total heat transport extending through the subtropical to subpolar latitudinal band (Figure 10c). In VO1 and VO2, the post-Tambora heat transports evolve qualitatively similar to AF, but the anomalies are weaker and only occasionally significant (Figure 10c). The changes in the North Atlantic meridional heat transport are reflected in the amount of heat transported by the ocean into the Nordic Seas across the Greenland-Scotland Ridge (Figure 9b). Likely explanations to the distinct behavior of AF reside in the background conditions, including the significantly weaker, compared to VO1 and VO2, pre-eruption AMOC state (Figures 1b and 9a), the significantly colder upper North Atlantic Ocean (Figures 1b and 8c), and a southward-shifted midlatitude zonal atmospheric circulation (not shown) influencing the wind-driven oceanic circulation and hence gyre strength (Figure 1b) and structure.
[29] In VO2, the net Arctic Ocean/Nordic Sea heat loss—defined as the total heat losses to the atmosphere from the Nordic Seas and the Arctic Ocean minus the heat transport into the Nordic Seas—only sporadically exceeds the natural variability range during the first two post-eruption decades (Figure 9d). This means that the post-eruption anomalous heat transport into the Nordic Seas (Figure 9b) is almost totally balanced by anomalous ocean heat release to the atmosphere over the Nordic Seas and the Arctic Ocean (the transport through the Bering Strait is neglected). In VO1 and especially in AF, the 1809 eruption induces anomalously strong Arctic Ocean/Nordic Sea heat losses (Figure 9d). After the Tambora eruption, despite the increased insulating effect from the extended post-eruption sea ice cover (Figure 9c), the oceanic heat loss to the atmosphere in the Arctic Ocean/Nordic Seas is not compensated by a possible oceanic heat gain from strengthened northward heat transport in the North Atlantic (Figure 10). In VO1 and AF, the post-eruption Arctic Ocean/Nordic Sea heat loss remains therefore largely outside the natural variability range (Figure 9d). Pre-eruption conditions are, again, primarily responsible for this anomalous behavior. In VO1 and AF, the Tambora eruption sustains an anomalous state induced by the 1809 eruption, while in VO2, it works against the descending phase of an internal fluctuation. Inter-ensemble differences in Arctic Ocean/Nordic Sea heat transport and balance (Figures 9b and 9d) are particularly relevant during the first 5 years after the Tambora eruption. About one decade after the eruption, the culmination of the post-eruption AMOC strengthening (Figure 9a) and the partially recovered lower-troposphere temperatures (Figures 7d–7f) contribute to a rebounding fluctuation, yielding a temporarily significant reduction of Arctic Ocean/Nordic Sea heat loss (Figure 9d).
[30] Overall, background conditions seem to actively influence the mechanisms involved in the post-eruption decadal climate evolution. The three ensembles, under indistinguishable imposed forcing (Figure 4), significantly differ in the magnitude and coherence of post-Tambora decadal climate signals. Depending on background conditions, after the perturbation, the system can therefore either evolve through trajectories mostly within the natural variability range (e.g., the VO2 ensemble) or be pushed largely outside it (e.g., the AF ensemble). An admissible interpretation is that background conditions modulate the strength of the feedbacks initiated by the imposed forcing, implying both different signals and different signal-to-noise ratios.
3.4 Implications of Limited Sampling
[31] The attribution of the different evolutions illustrated in sections 4 and 5 to a clustering of (purely) forced responses may be questionable due to the limited sampling allowed by the relatively small size of the ensembles. Accordingly, each ensemble-average trajectory may not be truly representative of the expected simulated climate response to the imposed forcing and under the given background conditions. A full sampling of the phase space is anyway not viable due to the complexity of the simulated climate. An alternative approach to assess limited sampling drawbacks on our inferences is to systematically investigate ensemble unperturbed evolutions in the control run for different ensemble sizes. Figure 11 proposes examples of results of a superposed epoch analysis performed on the unperturbed AMOC, where three ensembles are created by randomly sampling the eruption dates along the control run (the analysis is otherwise conducted as for Figure 9a). The figure demonstrates that especially, but not only, small ensemble sizes can lead to (a) clustering of ensemble trajectories, (b) trajectories nonrepresentative of the expected evolution, or (c) ensemble trajectories temporarily exceeding the natural variability range. Increasing the ensemble size certainly helps constrain the uncertainty around the expectation value but obviously does not guarantee against clustering or against exceed of the natural variability range (see, e.g., the 100-member ensemble). We cannot therefore unambiguously state that all the diagnosed simulated post-Tambora decadal behaviors are purely (volcanically) forced signals modulated by the background conditions. They may at least partly reflect the sampling of different natural oceanic evolutions. Of course, internal variability more strongly imprints on the post-eruption decadal climate evolution in the VO2 ensemble (e.g., Figure 9d). On the other hand, the forced component dominates the post-eruption variability in the AF ensemble, where ensemble anomalies largely exceed the natural variability range (e.g., Figures 8c and 9c). Whereas the sampling-related intrinsic ambiguity prevents fully disentangling the forced and unperturbed components of simulated post-eruption decadal climate variability, it does not alter our general conclusion that background conditions contribute in generating multiple pathways of simulated post-eruption decadal climate evolution.
3.5 Clustered Responses in an Idealized System
[32] Background conditions, and especially the initial ocean state, thus largely constrain the simulated decadal responses to the Tambora eruption. The diversification of post-eruption evolutions especially concerns processes that indirectly respond to the radiative perturbation through atmosphere-ocean coupled interactions and feedbacks. As a compendium to the high complexity of the simulated climate, Figure 12 sketches for the purpose of illustration how background conditions can affect the response to a perturbation in a simplified idealized system (see Appendix B). In short, the system is described by two independent and naturally oscillating processes (v1 and v2), both responding linearly to external forcing but with different sign and inertia, and by a third process (v3) which is a nonlinear combination of the other two and is not directly perturbed (see Appendix B). The left panels in Figure 12 illustrate the individual variable evolutions encompassing one full oscillation of the system both under unperturbed conditions (gray thick lines) and when an arbitrary series of four identical bell-shaped perturbations is imposed (black lines, forcing illustrated in the top panel). In all variables, the responses to individual perturbations encompass both negligible and prominent oscillations/fluctuations of opposite sign. The response is clearly influenced by the initial state as well as by the phase of the natural background oscillation. Looking at the perturbed evolutions only, it would be difficult to envisage the forced signal generated by the imposed four perturbations even for the simple linear response characterizing v1 and v2.
[33] The right panels of Figure 12 illustrate the perturbed system evolutions from an ensemble of events obtained by iteratively perturbing the system with the same bell-shaped forcing at each step within a full oscillation of the system. Mapped are probabilities of occurrence of perturbed states (described as post-perturbation anomalies, following a superposed epoch analysis) at each lag for individual variables. The possible post-perturbation states of v1 and v2 are near uniformly distributed around the mean (expected) response. The ensemble average thus meaningfully converges toward the expected response as the number of perturbations increases. Conversely, the nonuniform probability of possible states of v3 delineates a few preferred trajectories (i.e., the darker regions in the map) that only partially follow the mean response, which arises as being potentially deceptive. Accordingly, an appropriate sampling of background conditions would allow for the disclosure of multiple “typical” (forced) responses that are typical only under certain specific circumstances. In this aspect, this simple nonlinear system seems to qualitatively behave as the climate simulated by a full-complexity Earth system model.
4 Discussion
[34] In summary, our results support the hypothesis that differences in the background conditions (i.e., initial climate state and additional external forcings) can yield substantial externally induced modifications in the simulated post-eruption near-surface regional climate variability and ocean dynamics (Figures 2, 5, and 7-10). Our comparative assessment of ensemble-simulated and reconstructed surface temperature estimates for the early 19th century (Figure 2) confirm that, beyond accurate forcing estimates, knowledge about background conditions is pivotal for reproducing observed/reconstructed climate responses to individual historical volcanic eruptions [Zanchettin et al., 2013]. This complicates our assessment of the accuracy of volcanic effects as they are represented in numerical climate models and particularly the magnitude of the simulated cooling [e.g., Brohan et al., 2012]. If our hypothesis is valid, then the climate not necessarily always responded to individual historical volcanic eruptions as one would expect under the imposed forcing. Stepping back from a truth-centered approach for simulation-ensemble interpretation, i.e., assuming that each ensemble member is sampled from a distribution centered around the truth, in favor of a probabilistic one, i.e., assuming that each of the members is considered to be “exchangeable” with the other members and with the real system [e.g., Knutti et al., 2010; Bothe et al., 2012], is even more encouraged for the investigated period when constrains on the ensemble-simulated climate response to the Tambora eruption are further weakened by uncertainties on the 1809 eruption (for instance, the season of the eruption is still a tentative estimate; see Cole-Dai et al. [2009]).
[35] The North Atlantic/Arctic Ocean, the focus region of this study, turns out as symptomatic for the diversification of simulated decadal climate responses to a volcanic eruption. In particular, we note the impact of background conditions on sea ice properties, upper ocean heat content, and transport (Figures 8-10). A colder background global climate is a primary source of inter-ensemble differentiation, as this allows the volcanic radiative perturbation to trigger stronger feedbacks. The generality of this statement is, however, constrained by the large spread which characterizes the dependency of simulated responses on background conditions in an ensemble of different eruptions during the last millennium [Zanchettin et al., 2012a].
[36] Background conditions exert their effects in addition to a general response mechanism. The latter in our case holds under all circumstances and includes temporary strengthening of global ocean heat loss and of the thermohaline circulation, as well as a temporary increase in the northern hemispheric sea-ice cover. This general mechanism corresponds to the signal amplification described by Zanchettin et al. [2012a], and it allows for a consistently large signal-to-noise ratio in the ensemble response of winter regional climates to volcanic forcing. Therefore, the North Atlantic/Arctic sector simultaneously appears as a source of simulated post-eruption signal coherency and diversification. This statement is only seemingly self-contradicting. It simply highlights the dependency of our inferences on the experimental design. The ensemble-average post-eruption AMOC evolution indicates robust strengthening of the thermohaline circulation under all circumstances, i.e., in all ensembles (Figure 9a). Individual eruptions in climate simulations of the last millennium performed with the same model show a considerable range of possible AMOC responses, including post-eruption AMOC weakening [Zanchettin et al., 2012a]. Furthermore, as shown here, background conditions may affect our inferences on ensemble-simulated post-eruption behaviors as representative of a purely forced signal (compare Figure 11). We therefore interpret the robust AMOC strengthening as a specific feature of these simulations. Hence, our results do not necessarily contrast with hypotheses of volcanically forced decadal climate variability entailing an AMOC weakening for specific events [Mignot et al., 2011; Miller et al., 2012] or an overall nondetectable AMOC response [Hofer et al., 2011]. The simulated ensemble-average oceanic evolution of upper ocean temperature (Figure 8c) and oceanic net heat loss to the atmosphere in the North Atlantic (Figure 8d), and oceanic net heat loss in the Arctic Ocean/Nordic Seas (Figure 9d) consistently describe a ~20 year AMOC-type fluctuation. The hypothesis that volcanic forcing constrains the phasing of multidecadal variability thus seems to concern more processes and phenomena than those represented by the AMOC alone [Otterå et al., 2010; Zanchettin et al., 2012a].
[37] The robustness of post-Tambora AMOC strengthening contrasts with the diversity of post-Tambora ocean heat transport under different background conditions (Figures 9b and 10), which may, in turn, affect post-eruption hemispheric and, especially, regional climate evolutions (compare Figures 3 and 7). This means that inferences from the general climate system response to volcanic radiative perturbations may be too simplistic for inferences about regional decadal post-eruption climate variability. The large-scale response in, e.g., the AMOC or the global average temperature may be clear cut, but regionally the response may be dominated by other processes or feedbacks which depend more strongly on diverse boundary or background conditions. This was evidenced here by the effects of gyre circulation–driven changes in the total meridional heat transport in the North Atlantic (Figure 10b). Background conditions are not merely a source of additive noise for volcanically forced variability. A simplistic ensemble average may smear out relevant information about post-eruption dynamics. This, of course, applies in principle also to the within-ensemble variability of our ensembles. It follows that uncertainty in the background conditions may hinder conclusive assessments about the response pathways compatible with a certain near-surface temporal and spatial volcanic signature. So even as volcanic forcing sets the phase of some ocean processes, as for the AMOC, this may be of only marginal help for the dynamical interpretation of reconstructed regional climate variability [see also Zanchettin et al., 2012a, 2013] and, similarly, for decadal regional climate predictability if the background and boundary conditions are unknown.
[38] The Tambora eruption produces ensemble-average short-term near-surface hemispheric signatures depicting a robust large-scale atmospheric circulation anomaly only in the absence of additional external disturbances. The signature is confined to the American continent (Figures 7b and 7e) and corresponds to a positive PNA phase. Under full-forcing conditions, the post-Tambora winter atmospheric circulation is dominated by a thermodynamically driven drop of geopotential heights which stretches hemispherically (Figure 3e), and the positive PNA-like response becomes hardly distinguishable. Seasonally resolved data confirm the lack, in all ensemble averages, of apparent atmospheric dynamical features over the North Atlantic/European sector, including both short-term post-eruption winter warming [e.g., Fischer et al., 2007; Hegerl et al., 2011] and delayed winter warming [Zanchettin et al., 2012a, 2013]. The limited representation of stratospheric processes and stratospheric-tropospheric interactions in the Earth system model used here arguably affects the dynamics related to the short-term post-eruption winter warming, especially the downward signal propagation of the volcanically induced polar vortex strengthening [e.g., Zanchettin et al., 2012a]. Accordingly, the radiative surface cooling becomes more relevant for the ensemble-average tropospheric response. Despite known deficiencies like this one, the Earth system model used here can, however, simulate both features [Zanchettin et al., 2012a, 2013]. As exemplified for the short-term post-eruption winter warming (Figure 3), the fact that they are not represented in the ensemble-mean simulated responses can be interpreted as eruption and background specific rather than model specific. As shown for simulated and observed global/reconstructed hemispheric temperatures in the early 19th century, the interplay between external forcing and internally generated variability can imprint uniquely also on decadal-scale climate evolutions. The decade around 1830 is exemplary in showing the possible amplitude of the range of simulated decadal anomalous states and their different agreement with the reconstructed state (Figure 2b). It remains therefore debatable whether simulated ensemble-average responses are appropriate for investigating dynamical climate features and/or for validating ensemble evolutions against a target.
5 Conclusions
[39] Our systematic investigation of post-Tambora climate evolution in a set of Earth system model simulation ensembles focusing on the early 19th century climate clearly indicates that uncertainty on background conditions (including initial climate state and presence and strength of additional external forcings) complicates the assessment of simulated decadal climate responses to short-lived radiative perturbations such as those induced by strong tropical volcanic eruptions. The effects of background conditions superpose on a general decadal response mechanism involving primarily temporary strengthening of ocean-to-atmosphere heat fluxes and strengthening of the thermohaline circulation. This superposition is, however, not merely adding noise to the forced response. Rather, it results in a clustering of simulated post-eruption decadal climate evolutions, particularly of ocean heat transport/budget and sea-ice cover. Acknowledging the existence of deterministic background influences on forced responses implies that background conditions must be appropriately accounted for in the dynamical interpretation of reconstructed and observed regional climate variability and, similarly, in the design of decadal climate prediction tools.
Acknowledgments
[43] The authors thank Nils Fischer who did the internal review at MPI-M, and Jochem Marotzke, Max Popp, and Daniela Matei who provided useful comments on earlier versions of the manuscript. We also thank three anonymous reviewers for their constructive comments. This work was carried out as part of the MPI-M integrated projects “Millennium” and “Super Volcano.” This work was supported by the Federal Ministry for Education and Research in Germany (BMBF) through the research program “MiKlip” (FKZ:01LP1158A(DZ):/01LP1130A(CT)) and by the European Union Seventh Framework Program under grant agreement n.212643 (EU-THOR). O.B. was supported through the Cluster of Excellence “CliSAP,” University of Hamburg, funded through the German Science Foundation (DFG). J.L. acknowledges support from the DFG Projects PRIME 2 (Precipitation In past Millennia in Europe-extension back to Roman times) within the Priority Program “INTERDYNAMIK.” J.L. also acknowledges support from the DFG Project “Historical Climatology of the Middle East based on Arabic sources back to AD 800” and from the EU/FP7 project ACQWA (NO212250).
Appendix A:: Ocean basins
[40] Figure A1 illustrates the ocean basins as they are defined in this study: the North Atlantic Ocean, the Nordic Seas, and the Arctic Ocean.
Appendix B:: The arbitrary idealized system
(1)
(2)
(3)
(4)
[42] In the absence of forcing F, v1 and v2 are sinusoidal functions of t with parameters a1 = 1, a2 = −1, ω1 = 0.02, ω2 = 0.03, and ϕ1 = ϕ2 = 0. Both variables respond linearly to F but with opposite sign. Inertia in the response of v1 to F is implemented in equation 1 by smoothing F with a backward moving average. F consists of a series of i events (four in Figure 12), each described by a Gaussian function as in equation 4, with parameters a3 = 1.5, β = 500 governing its amplitude and duration, and α determining its position in time (in Figure 12, α = 150, 300, 550, 700). Variable v3 is a nonlinear combination of v1 and v2, and does not respond directly to F.
