The unidentified volcanic eruption of 1809: why it remains a climatic cold case

The “1809 eruption” is one of the most recent unidentified volcanic eruptions with a global climate impact. Even though the eruption ranks as the 3rd largest since 1500 with an eruption magnitude estimated to be two times that of the 1991 eruption of Pinatubo, not much is known of it from historic sources. Based on a compilation of instrumental and reconstructed temperature time series, we show here that tropical temperatures show a significant drop in response to the ~1809 eruption, similar to that produced by the Mt. Tambora eruption in 1815, while the response of Northern Hemisphere (NH) boreal sum5 mer temperature is spatially heterogeneous. We test the sensitivity of the climate response simulated by the MPI Earth system model to a range of volcanic forcing estimates constructed using estimated volcanic stratospheric sulfur injections (VSSI) and uncertainties from ice core records. Three of the forcing reconstructions represent a tropical eruption with approximately symmetric hemispheric aerosol spread but different forcing magnitudes, while a fourth reflects a hemispherically asymmetric scenario without volcanic forcing in the NH extratropics. Observed and reconstructed post-volcanic surface NH summer tem10 perature anomalies lie within the range of all the scenario simulations. Therefore, assuming the model climate sensitivity is correct, the VSSI estimate is accurate within the uncertainty bounds. Comparison of observed and simulated tropical temperature anomalies suggests that the most likely VSSI for the 1809 eruption would be somewhere between 12 -19 Tg of sulfur. Model results show that NH large-scale climate modes are sensitive to both volcanic forcing strength and its spatial structure. While spatial correlations between the N-TREND NH temperature reconstruction and the model simulations are weak in terms 15 of the ensemble mean model results, individual model simulations show good correlation over North America and Europe, suggesting the spatial heterogeneity of the 1810 cooling could be due to internal climate variability.

. Volcanic radiative forcing. Global stratospheric aerosol optical depth (SAOD) at 0.55 µm based on eVolv2k VSSI estimates (Toohey and Sigl, 2017) and calculation with the volcanic forcing generator EVA (Toohey et al., 2016) for the four different forcing scenarios: "Best", "Low", "High" and "nNHP" for the 1809 eruption and the "Best" forcing scenario for the Mt. Tambora eruption. Bottom: Spatial and temporal distribution of a zonal mean stratospheric SAOD for the four experiments. reconstruction of the spatial forcing structure of volcanic aerosol is important to get a reliable climate response. In fact, given the current knowledge, an interhemispheric symmetric forcing is only one of several realistic possibilities for the 1809 eruption.
In recognition of these additional sources of uncertainty (and motivated in part by the apparent muted temperature response   (Brohan et al., 2012) and Indo-Pacific warm pool data (D'Arrigo et al., 2006). b) NH summer land temperatures from four tree-ring based reconstructions (Wilson et al. 2016 (N-TREND(N)), Anchukaitis et al., 2017 (N-TREND (S)), Guillet et al., 2017(NVOLC), Schneider et al., 2015). c-d) Monthly mean NH winter (c) and summer (d) temperature anomalies ( • C) from 53 station data averaged over different European regions   interior Alaska/Yukon (Briffa et al., 1994;Davi et al., 2003;Wilson et al., 2019), 1810 is one of the coldest summers identified over recent centuries. In earlier reconstructions of summer SAT in different regions of the western United States (Schweingruber et al., 1991;Briffa et al., 1992), 1810 was shown to be the third coldest summer in the British Columbia/Pacific Northwest region. Likewise, European tree-ring records show cooling after 1809 (e.g., Briffa et al., 1992;Wilson et al., 2016). In contrast, 260 tree-ring networks in certain regions such as Eastern Canada show minimal response after 1809 (Gennaretti et al., 2018).
Based on compilations of regional records, TR-based reconstructions of NH-mean land summer SAT show a large spread in hemispheric cooling after the 1809 eruption (Fig. 2b), with anomalies of -0.87 • C, -0.77 • C, -0.21 • C and -0.15 • C in 1810 for the N-TREND (S), NVOLC, N-TREND (N) and SCH15 reconstructions, respectively). Although using the same data set, the spatial N-TREND (S) and the nested N-TREND (N) reconstructions show a quite different behaviour. In N-TREND (S), 265 the nature of the spatial multiple regression modelling biases the input records to those that correlate most strongly with local temperatures which, when available, were are likely MXD data. In all four reconstructions, NH temperature does not return to the climatological mean after an initial drop in 1810, but remains low, or even exhibits a continued cooling trend until the Mt.
Tambora eruption in 1815 . The spatial variability of the reconstructed NH extratropical temperature response to the 1809 eruption is illustrated in Fig. 2d,e, based on the spatially resolved N-TREND (S) reconstruction (An-270 chukaitis et al., 2017), displaying zonal oscillations consistent with a "wave-2" structure, that are especially evident in 1810 but already appreciable in 1809. This hemispheric structure is in contrast with the relatively uniform cooling seen in TR records for Tambora (Fig. 2f) and indeed for many of the largest eruptions of the past millennium . Information about regional and seasonal mean NH temperature anomalies in the early 19th century can be obtained from different station data across Europe and from New England (Fig. 2c). In NH winter the measurements reflect the high variability 275 of local scale weather (Fig. 2c). Warm anomalies, an indication for post-eruption "NH winter warming", are clearly visible in 1816/1817 in the 2nd winter after the Mt. Tambora eruption in 1815. Northern Europe shows the largest warm anomaly for all regions (about 3 • C). Warm NH winter anomalies between 1.5 • C and 2 • C are seen in the winter 1809/1810 over Northern and Eastern Europe and over New England. Strong cooling is however found for the 1808/1809 winter in Northern and Central Europe. NH summer temperature anomalies are less variable than in winter (Fig. 2d). A local distinct cooling is found in Interestingly, warm anomalies in the order of 2 • C are found in summer 1811 over Central and Eastern Europe.

Simulations
Firstly, we compare the simulated evolutions of monthly mean near-surface (2m) air temperature anomalies between the four experiments globally, in the tropics and NH extratropics (Fig. 3). Ensemble mean global mean temperature anomalies grow through 1809, and reach peak values through 1810 in all experiments, before decaying towards climatological values (Fig. 3a).
Peak cooling reaches around 1.0 • C in the High experiment, compared to 0.5 • C in the Low and nNHP experiments. Peak 290 temperature anomalies across the experiments correlates with the magnitude of prescribed AOD (Fig. 1a), and the responses are qualitatively consistent with expectations, with the AOD for the Low and nNHP experiments, which is similar in magnitude to that from the observed 1991 Pinatubo eruption, leading to global mean temperature anomalies also similar to those observed after Pinatubo. Global ensemble mean near-surface temperature anomalies are close together in Low and nNHP over boreal summer but differ for boreal winter when the intrinsic variability is higher. Low is the only experiment for which large-scale 295 temperatures return within the 5th-95th percentile range percentile interval of the control run before the Mt. Tambora eruption extra-tropical northern hemisphere and c) tropical averages of near-surface air temperature with respect to the pre-eruption (1800-1808) climatology. All data are deseasoned using the respective annual average cycle from the control run. Thick (thin) black dashed lines are the 5th-95th percentile intervals for signal occurrence in the control run for the ensemble mean (ensemble spread). Bottom bars indicate periods when an ensemble member's monthly mean temperature (color code as for the time series plots) is significantly different (p = 0.05) from the control run according to the Mann-Whitney U test. Right: ensemble distributions (median, 25th-75th and 5th-95th percentile ranges) of seasonal mean anomalies for the first post-eruption winter (1809-1810, DJF) and summer (1810, JJA) following the 1809 eruption as well as for the pre-eruption period (1800-1808).
in 1815. Global mean temperature anomalies of the other three experiments return only within the 5th-95th percentile range of unperturbed variability by 1815. As expected, almost no significant near-surface temperature anomalies are found for the nNHP simulation in the NH extratropics except a few months in spring and autumn 1813 (Fig. 3b). The nNHP ensemble-mean values stay within the interquartile range of the control run but show a slight negative trend between 1809 and 1815. nNHP is 300 also the only experiment where the NH extratropical summer of 1814 is colder than the summer of 1809. Internal variability is relatively high in the NH extratropics in particular in NH winter spanning more than 1.5 • C. So, even the ensemble mean near surface temperature anomalies for Best and High almost reach the 5th-95th percentile range of the control run in the 1st post volcanic winters. Peak cooling appears for all experiments except nNHP in the summer 1810. In the tropics, Best, High and nNHP are outside the 5th-95th percentile range in the first four -post volcanic years while Low exceeds the 5th-95th percentile 305 range only for two years (Fig. 3c).
The ensemble distributions for the seasonal mean of winter 1809/1810 and summer 1810 illustrate the differences between the four experiments not only in the mean anomaly but also for the ensemble spread (  found in the tropics (30 • S-30 • N). In the inner tropics, the cooling disappears after one and a half years in Low and two to three years later in Best, High and nNHP. In the subtropics, a significant surface cooling signal is found over the ocean until 1815 in Best and High, while over land no significant cooling appears in 1814 (Fig. S2). A strong cooling signal is found in the NH extra tropics in Best, Low and High in summer 1810 and in High and to a small extent also in Best in summer 1811. In nNHP no surface cooling is detectable over the NH extratropics in the 1st four years after the eruption, consistent with the prescribed 320 volcanic forcing (see Fig. 1). However, a cooling anomaly is apparent around 60 • N in summer 1813, which is seen in the zonal mean over the ocean (Fig. S2) and likely due to decreased poleward ocean heat transport. Significant cooling south of 30 • S appears only in austral spring 1809.   latter with a narrow distribution as shown by the tight whisker plot. Although the PNA is mostly active in NH winter and is therefore negligible for the NH summer circulation, the weakening of the summer Aleutian Low is expressed by the positive 360 NPI index in nNHP which is compatible with the Western North American cooling in the reconstructed temperature pattern (Fig. 3). The SOI shows significant post-eruption anomalies in the High ensemble in winter (negative) and in the High and Best ensembles in summer, the latter with opposite signs (Fig. 6d). Overall, the results illustrate the substantial differences in the post-eruption evolution of continental and subcontinental climates that can be produced by internal climate variability and forcing structure through changes in the large-scale atmospheric circulation, as seen by the spread of responses within 365 individual ensembles (often as large as the range of pre-eruption variability) and by the possibility of non-overlapping response distributions generated by different forcings (seen for, e.g., the NPI and SOI indices in summer 1810). experiments.In 1810-1812, the cooling in the Best, High and nNHP experiments is stronger than that observed, and therefore 375 the results from the Low experiment are generally most consistent with the ship-borne measurements (Fig. 7a). A comparison of TROP with our four experiments reveals that all experiments lie within the 5th-95th percentile interval of the TROP reconstruction although the reconstructed SST response appears to be dampened in comparison to the model experiments (Fig.   7b). Although the long term trends of TROP and the model experiments are in general agreement, the dampened post volcanic cooling could reflect autocorrelative biases in the proxies (Lücke et al., 2019). Detailed scrutiny of high resolution tropical 380 SST proxies and their potential biases to robustly reflect volcanically forced cooling has not been made in the same way that has been performed for tree-ring archives over the last decade (Anchukaitis et al., 2012;D'Arrigo et al., 2013;Esper et al., 2015;Franke et al., 2013;Lücke et al., 2019). A similar behavior is found for the Indonesian warm pool (Fig. 7c). However, in contrast to the whole tropics, the differences between the different forcing experiments are much smaller for the warm pool region compared to the wider tropical regions and the volcanic signal is more pronounced in the reconstructed SST at least for 385 the unidentified 1809 eruption. Best show a distinct cooling in 1809 and nNHP in 1811. In contrast to the East Pacific, variability in the West Pacific is rather small (Fig. 8b). In all four experiments a clear volcanic signal is visible in the simulated ensemble mean SST anomaly after the 1809 eruption and the Mt. Tambora eruption, whereas only a weak signal appears for both eruptions in the reconstruction.

395
Interestingly, in the West Atlantic, two distinct positive SST anomalies appear in the reconstructions in the aftermath of the unidentified 1809 and the Mt. Tambora eruption, while the MPI-ESM simulations shows cooling (Fig. 8c). Reasons for the anticorrelated behaviour are not obvious per se and may be related to changes in either ocean circulation or other climate factors than SST that influence the coral record, such as salinity and precipitation. In the reconstruction, the Indian Ocean is the only region that displays a peak cold anomaly after the Mt. Tambora eruption, but the magnitude of this cooling is comparable to an 400 apparent cooling in 1807. A clear reference to the Mt. Tambora eruption is therefore difficult to establish. No large cooling is found in the coral data after 1809 over the Indian Ocean (Fig. 8d). show that aerosol quickly spreads uniformly across the tropics, it is unlikely that aerosol forcing from the 1809 eruption would be significantly different between the tropical Indian and Atlantic ocean basins. Therefore, differences in temperature response 420 in the model between the two regions seems more likely to be related to model sensitivity, which might particularly be linked to differences in ocean circulation and/or mixed layer depth. the Best experiment (Fig. S3) suggesting the influence of internal variability. For NH winter, both model and station data show higher variability than in NH summer (Fig. 14). Simulated and observed NH winter temperature anomalies agree quite well in the first three winters after the 1809 eruption. The only exception is New England (Fig. 14f) where, similar to NH summer 510 (Fig. 13f), less agreement is found between the station data and the four experiments. The strong cooling signal of more than -2 • C which is found at Northern, Western and Central European stations in winter 1813 /1814 is not reproduced by the model  et al., 2002b), and is reproduced by some forced simulations and years from the control run, suggesting it is a plausible result of natural variability. In contrast, the reconstructed temperature pattern over Asia is not produced by any model simulation, forced or control. Possible explanations for this model-data disagreement include deficiencies and uncertainties regarding both tools. In particular, the dendrochronological network remains weakly constrained over the central Asian region, especially with 625 respect to available MXD or density-like parameters, and instrumental observations are sparse, especially prior to the 1950s making calibration and validation difficult (Cook et al., 2013). Anchukaitis et al. (2017) clearly detailed poor reconstruction validation for many regions across North America as well as Northeast and central Asia. These regions must be targeted in ongoing proxy network development to update and develop new MXD or density-like parameter data-sets which are proven to be superior proxies of volcanically forced summer cooling (e.g., Anchukaitis et al., 2012;D'Arrigo et al., 2013;Esper et al., 630 2015; Lücke et al., 2019). Concerning the simulations, beyond model deficiencies that may bias the response, our study demonstrates that choices regarding both volcanic forcing strength and spatial structure can similarly affect reconstruction-simulation comparative assessments. Specifically, it was shown that for the muted response to the 1809 eruption in the NH extratropics, agreement between reconstructions and simulations improves by weakening the magnitude of the eruption and, alternatively, by preventing the volcanic aerosol to spread over the Northern Hemisphere. The forcing estimate of 1809, which is currently 635 based on ice core data, can only be improved by modelling experiments to narrow down the uncertainty range. This will be facilitated if further information on location and eruption season is identified. At a broad-scale, latitude can be identified using S isotopes (Burke et al., 2019) which could at least help constrain whether nNHP experiments are potentially valid to improve the characterization of the climatic imprint of the 1809 eruption. Finally, it is clear that an increase in reconstruction accuracy by improving spatial coverage, including also the southern hemisphere, is needed. Otherwise this eruption will remain most 640 likely a mystery and a climatic cold case.

Northern hemisphere extratropics
respective climatological value, with the suffix [x] for the region x. The NAO index is calculated based on the latitude-longitude two-box method from Stephenson et al. (2006) applied on Z500 data, i.e., as the pressure difference between spatial averages over 20-55 • N; 90 • W-60 • E and 55-90 • N; 90 • W-60 • E. The NPI is calculated using the definition from Trenberth and Hurrell (1994) applied to sea level pressure (SLP) data. The index is computed as spatial SLP averaged over 30-65 • N; 160-220 • E, so that positive phases of the index indicate a weaker-than-normal Aleutian low and the opposite holds for the 650 negative phases. The SOI is defined as the difference between the average SLP over the domains 20-15 • S; 147-152 • W and 15-10 • S; 128.5-133.5 • E.

Source
Station Start year End year lat lon N