Evaluation of the 1816 Year Without a Summer in the Last Millennium Ensemble

Model performance is evaluated (see “Methods” section) with satellite retrievals and reanalyses based on late 19th and 20th century volcanic eruptions. The LME is found to reproduce key features of the satellite-estimated net top-of-atmosphere (TOA) radiative imbalance arising from the recent 1991 eruption of Mt. Pinatubo25 and basic aspects of volcanic responses in June through August (JJA) in the year following eruptions estimated from reanalyses, Here, we also use updated reconstructions in Europe and Asia (Fig. 1, refs. 10,27), where anecdotal accounts suggest the eruption respectively drove strong cooling and drought, to evaluate model fidelity.

Fig. 1 Proxy reconstructions and LME-simulated responses compared. European 1816 summer temperature anomalies relative to 1990–1999 in a proxy reconstructions26 and b from the LME, and c southeast Asian Palmer Drought Severity Index (PDSI), an indicator of drought10, and d their ensemble means over contemporaneous intervals from the LME. Stippling indicates regions where reconstructions lie outside of the ensemble model range for both temperature and PDSI Full size image

Simulating the 1816 Year Without a Summer

Reconstructions spanning JJA of 1816, the YWAS, show regions of agreement and disagreement with the LME forced response, obtained by averaging across 18 volcanically forced LME ensemble members (13 fully forced and 5 volcanic-only forced, Fig. 1). Reconstructed surface temperature anomalies are characterized by cooling as large as −4 K in western Europe that exhibits strong zonal structure and becomes small east of about 20°E, with warming in eastern Europe. The LME forced response exhibits both similarities, such as cooling in western Europe, and differences, including only modest cooling in northern Europe and a general lack of warming in eastern Europe. These contrasts between reconstructions and the forced response may arise from error in the reconstructions, error in the model forcing18 and forced response, or the influence of internal variability16,18. Given the uncertainties arising from internal variability and errors in both reconstructions and models neither can be taken as definitive28 and so the goal here is to assess the general consistency between the two. This reconstruction is based on a data set that includes a large number of homogenized and quality-checked instrumental data series, a number of reconstructed sea-ice and temperature indices derived from documentary records for earlier centuries, and a few seasonally resolved proxy temperature reconstructions from Greenland ice cores and tree rings from Scandinavia and Siberia27. A separate spatial reconstruction of European 1816 summer temperatures, based solely on instrumental data, indicates anomaly patterns much like the reconstruction shown in Fig. 1a 26, suggesting that model forcing and error, and the influence of internal variability, are the most salient sources of contrast, whereas other reconstructions29 agree more closely with the model forced response, suggesting a role for internal variability. Regarding the latter, a role for the negative phase of the North Atlantic Oscillation (NAO30) and its observed summer teleconnections (Fig. 2) is suggested as it contributes to anomalous cooling in northern Europe and enhances the zonal temperature gradient, helping reconcile these model-reconstruction differences to a degree. The LME mean NAO anomaly in the summer of 1816 (i.e., its forced response) is small but statistically significant (−0.21 ± 0.18, 5% conf.), while the actual value in nature varies considerably across reconstructions, particularly in summer. Uncertainties surrounding the models, proxies, and NAO reconstructions preclude a precise determination of its contribution.

Fig. 2 Regional temperature structure of the NAO. Regression of the negative phase of the NAO during June through August against contemporaneous surface temperature for the region considered in Fig. 1 (see “Methods”). Contributions to cooling from the negative phase of the NAO are greatest in regions where the model forced response fails to reproduce the cooling estimated in reconstructions Full size image

In southeast Asia, hydroclimate reconstructions10 show regions of strong drought spanning India, Indonesia, and southeast Asia, bounded to the north by wetter than normal conditions in Pakistan and western China. While uncertainty surrounding these reconstructions also remains considerable28,31, the degree of agreement with the LME forced response is notable particularly given the large spread across individual realizations in the model (see Supplementary Discussion and Supplementary Figs. 1, 2).

Together the reconstructions in Fig. 1 suggest the skillful reproduction of key aspects of the YWAS climate response in the LME, particularly for Asian drought and cooling in portions of southern Europe, and underscore the eruption’s global reach and the need to develop a suitable large-scale process-based understanding of the event in order to begin to understand its modulation in a warmer climate. The LME simulations provide considerable insight into the main energetic responses governing surface impacts. For example, a useful measure of the eruption’s radiative effect is provided by the global clear-sky shortwave flux through the atmosphere (SW-Catm), computed as the difference between clear-sky downwelling solar radiation at the surface and TOA. This is shown in conjunction with simulated global mean 2-m air temperature (T) and precipitation (P) anomalies (Fig. 3a), and associated TOA radiative flux anomalies (Fig. 3b). In the months following the eruption, a decrease in SW-Catm of about 16 W m−2 is simulated, initially in the tropics but spreading globally in the months thereafter (Supplementary Fig. 3). The decrease in TOA absorbed shortwave radiation (SW) exceeds that in outgoing longwave radiation (also Supplementary Fig. 4) by a factor of about 2, such that a net TOA deficit of about 6.5 W m−2 is realized in mid-1815. These simulated radiative responses are in general agreement with those derived in other recent model ensembles32,33. In response to this radiative deficit, the surface cools and the hydrologic cycle weakens. By late 1815, T drops by 0.74 K in the ensemble mean, an amount on par with the climate system’s net warming in the entire 20th century, and rainfall drops by 0.085 mm day−1 (2.8%), as less sunlight reaches the surface and reduces the energy available for evaporation.

Fig. 3 Global mean-simulated responses to the 1815 eruption. Simulated anomalies relative to 1814 are shown in a global mean clear-sky shortwave flux through the atmosphere (SW-Catm), 2-m air temperature (T), and precipitation (P), and b SW-Catm, global mean TOA incoming net (R T ), shortwave (SW), and outgoing longwave (LW) radiation. Shading denotes the range of two standard errors. Years indicated on the abscissas begin at the corresponding tick marks Full size image

Cooling is particularly strong in the Northern Hemisphere high latitudes (Fig. 4a, b), where cryospheric responses, including increases in snow cover over land and sea ice over ocean, extend the persistence of aerosol-induced albedo increases (Supplementary Fig. 5). The simulated anomalies are broadly consistent with other recent model ensembles in their depiction of pervasive northern hemisphere cooling over land32,33. Anecdotal accounts of the YWAS associate cool temperatures with the lingering radiative effects of volcanic aerosols6 but in the LME, SW-Catm anomalies have largely dissipated by the summer of 1816 (Fig. 3, Supplementary Fig. 3). Rather, the simulations show that it is mainly a cold upper ocean, associated land amplification (discussed further below), and other responses that set the stage for the anomalously cool YWAS conditions. In addition to cryospheric changes, a forced El Niño event (Fig. 4b) is realized (in 12 of the 18 ensemble members) and it is associated with positive rainfall anomalies over land (Fig. 4c), contributing to additional cooling in the western US and southern Europe.

Fig. 4 Regional-simulated responses to the 1815 eruption. a Zonal mean T for land and ocean (L + O), and land, and ocean, separately and YWAS July through September near-surface winds and b T and c P. Mapped colors (b, c) and hatching (a) indicate anomalies exceeding twice the ensemble standard error. Stippled regions in a indicate anomalies between one and two times the standard error. Years indicated on the abscissas begin at the corresponding tick marks Full size image

Despite the severity of the simulated weakening in the hydrologic cycle, energy budget perturbations suggest that impacts of the eruption may have been much larger had it not been for the mediating role of the ocean. For example, a peak net TOA flux anomaly of −6.5 W m−2 is simulated following the eruption (Fig. 3b), a deficit that is roughly equal to 8% of the global mean evaporation (~85 W m−2) that drives the water cycle, yet the cycle’s weakening is simulated to be less than 3%, an amount on par with the simulated weakening in another recent ensemble33. One potential explanation for the discrepancy is that responses in other terms of the surface energy balance, such as in sensible heat or longwave radiation, compensate for the imposed surface shortwave perturbation and mitigate the surface energetic imbalance so as to lessen reductions in evaporation. However, in these other terms, simulated responses (shown in Fig. 5) act to accentuate shortwave cooling anomalies rather than lessen them, though the magnitude of each is relatively small.

Fig. 5 Responses in terms of the surface energy budget and OHC to the 1815 eruption. Changes in terms of the surface energy budget over ocean (left axis) following the 1815 Mt. Tambora eruption including net upwelling longwave flux (LW), sensible heat (SH), latent heat (LH), and net surface shortwave absorption (SW). Also shown is the change in full-depth OHC (right axis) Full size image

LME simulations demonstrate instead that, as suggested by land-sea contrasts in the temperature response (Fig. 4a), ocean cooling plays a key role in mediating the eruption’s surface response, with a peak OHC reduction of 7 × 1022 J in 1817 (Fig. 5), followed by a gradual recovery, particularly in the deep ocean, lasting over a decade (Figs. 5, 6, Supplementary Fig. 6). The associated heat flux out of the ocean in the years following the eruption thus played a key role in sustaining global evaporation near pre-eruption levels and limiting the reduction in global mean rainfall. The long timescale of OHC recovery suggests the potential for a continuing influence of the earlier 1809 eruption34 on OHC contemporaneous with the 1815 Tambora response. LME simulations show that the associated surface anomalies had largely dissipated by 1815, however.

The simulated ocean cooling response is particularly pronounced at the surface in the subtropics (Fig. 6, Supplementary Fig. 6) and extends to depth on the poleward fringes of the subtropical overturning cells, where subduction and ventilation in the upper 1000 m are strong35 (Supplementary Fig. 7). In the upper 100 m, cold anomalies peak in early 1816 (Fig. 6, Supplementary Fig. 6) and cooling extends to depth (Fig. 6), underpinning the cool surface, particularly equatorward of the main YWAS regions in North America and Europe (Fig. 4b, Supplementary Fig. 6). Surface ocean temperatures drop disproportionately in the Northern Hemisphere (Fig. 4a, b), where integrated OHC anomalies are comparable to those in the Southern Hemisphere but ocean extent is less.

Fig. 6 Responses in OHC with depth to the 1815 eruption. a Evolution of SW-Catm (left axis) and OHC (right axis) across various ocean layers. b Anomalies in zonal mean ocean temperature by depth and latitude in 1816, with stippled regions indicating a consistent sign of anomalies across at least 14 of the 18 (80%) ensemble members. Shading in a denotes the range of two standard errors. Years indicated on the abscissas begin at the corresponding tick marks Full size image

The forced response of LME simulation to the 1815 eruption reproduces many of the YWAS accounts, such as drought in Asia and cooling in northeastern North America and western Europe (Fig. 4), and suggests several others. By January through March (JFM) of 1816, land regions south of 60°N are simulated to have cooled generally, with particularly strong cooling in central North America and southern Eurasia. The Aleutian low had also strengthened, and the tropical Walker circulation had weakened (Supplementary Figs. 8, 9). The weakened Walker circulation is consistent with strong reductions in rainfall surrounding the Maritime continent (Fig. 4, Supplementary Fig. 9) and a transition in the atmosphere, such that during the YWAS, El Niño SST anomalies, and associated cloud anomalies emerge16 (Supplementary Figs. 8, 10). On average, the response of clouds to the eruption is characterized both by regional anomalies, associated for example with a transition into El Niño conditions in 1816–1817 (Supplementary Fig. 10), and a global reduction in cloud amount (Table 1) and, specifically tropical high cloud amount, implying a positive feedback in response to the eruption due to increased LW emission to space.

Table 1 Cloud feedbacks under climate change and in response to the Tambora eruption Full size table

Simulating an analogous 2085 eruption

The large-scale perturbations driven by Tambora’s 1815 eruption (as demonstrated above for example in its oceanic, surface, and atmospheric responses) raises the question as to what the modulating influence of climate change might be if such an event were to unfold in the future. Here, we therefore compare and contrast key large-scale aspects of a hypothetical 2085 eruption (see “Methods”) with the 1815 response (Fig. 7). Many aspects of the radiative response closely resemble those during the 1815 event and these include a dramatic decrease in SW-Catm, an increase in clear-sky albedo at low latitudes (Supplementary Fig. 11), and reductions in TOA net SW and LW (Supplementary Figs. 12, 13) fluxes. Differences in the overall magnitude and spatial structure of the SW radiative responses are generally small. At the surface, however, a future eruption’s peak response in both global mean temperature and rainfall is simulated to be approximately 40% greater than in 1815 (Fig. 7a), though the uncertainty associated with the ensemble spread is also considerable and limits a precise quantitative estimation of the increase. Changes in the background climate state are likely to be instrumental in driving this increase. Key among the background changes is the fact that the upper ocean has warmed considerably in the simulated future climate, particularly in the tropics, contributing to an increase in ocean stratification (Fig. 7b). Much of the increase in surface and upper ocean cooling following the 2085 eruption is simulated at low latitudes in winter and spring, and in northern polar regions, where the cryospheric response increases in magnitude and is displaced poleward due to the warmer base state, thereby moderating its radiative impact such that the net change in global surface albedo is small (Supplementary Fig. 14). As a consequence of increased ocean stratification, temperature anomalies in the upper ocean do not penetrate to depth as efficiently as in 1815 (Fig. 7d), leading to greater surface cooling and reduced cooling at depth at all latitudes except from 30° to 45° in both hemispheres. Increases in upper-ocean stability are evident on a global scale and persist throughout the annual cycle. Changes in wind speed, salinity, and mixed layer depth are also likely to influence the depth of cooling, though the spatial distribution of these fields is more complex and seasonally dependent than for stability (Supplementary Figs. 15, 16). Quantifying the relative roles of each in modulating the volcanic response remains a work in progress. At high latitudes, SW responses strengthen during summer as the polar oceans have moved toward a more ice-free state allowing the formation of sea ice induced by the eruption to have a disproportionate effect on albedo (Supplementary Fig. 11). Despite its enhanced magnitude relative to the 1815 eruption, the surface cooling induced by a 2085 eruption (~1.1 K, Fig. 4a) remains much less than the associated warming of the background state (4.2 K) and hence the magnitude of this amplified, historically exceptional volcanic event is nonetheless unable to offset even a third of the accumulated transient background warming since the 1815 Tambora eruption.

Fig. 7 Changes in eruption responses in surface and ocean temperatures between 1815 and 2085. Ensemble mean-simulated a anomalies in global mean surface temperature (K) and rainfall (mm day−1, dashed, right axis) for 1815 (blue) and 2085 (red) eruptions, smoothed with a 12-mo running mean and shown with a 2 standard error range for temperature (shading), b annual mean net ocean warming by latitude and depth estimated from the 2085–2094 change from 19th century ensemble zonal mean, c differences in ensemble zonal mean 2-m air temperature anomalies between the 2085 and 1815 eruptions (stippled where anomalies exceed twice the standard error), and d the 2086–1816 difference in ensemble zonal mean ocean temperature anomalies by latitude and depth. Years indicated on the abscissas begin at the corresponding tick marks Full size image

Lastly, it is worth noting that volcanic eruptions have been proposed previously as a means of understanding the climate system’s base sensitivities. One key context relates to climate sensitivity and cloud feedbacks in response to forcing, such as for example increases in greenhouse gases36,37. Here, responses in the cloud field to volcanic eruptions and a warming climate are found to differ significantly in key respects. By regressing cloud fields against global mean surface temperature (Table 1) during both eruptions and transient background climate change, large contrasts are identified. Generally, simulated cloud changes are much greater per degree cooling following volcanic eruptions than those that accompany anthropogenic climate change. For example, under greenhouse warming, estimated from 1815 to 2085 differences, total tropical cloud amount decreases by −0.41% K−1 yet following volcanic eruptions, it also decreases (with cooling) by about 1.26% K−1 and changes are thus not only substantially different in magnitude but opposite in sign. For tropical high cloud the contrast is particularly large with reductions with greenhouse gas warming of about −0.14% K−1, whereas following volcanic eruptions high clouds decrease in extent by 2.33% K−1. As the primary uncertainty in determining climate sensitivity lies in quantifying the response of such clouds to warming38,39, these discrepancies suggest that attempts to constrain cloud feedback and climate sensitivity uncertainties through observations of cloud and radiative responses following volcanic eruptions is unlikely to provide a useful constraint. Volcanic eruptions may, however, offer insight into the capabilities and limitations of models in evaluating geoengineering strategies and associated follow-on effects. Continuing to confront climate models with a range of climate phenomena across timescales, including major volcanic eruptions of the last millennium, remains a useful exercise for improving our understanding of climate as it changes beyond the range of observed variability in Earth’s recent history.