1 Introduction

As recorded in polar ice cores, the largest volcanic sulphur release to the atmosphere in at least the past 2500 years occurred between 1257 and 1259 CE (hereinafter all calendar dates refer to the Common Era; CE) (Palais et al. 1992; Oppenheimer 2003a; Sigl et al. 2015). The source of this eruption was eventually identified as Samalas volcano, part of the Rinjani complex on Lombok Island, Indonesia (Lavigne et al. 2013). Though the eruption’s great magnitude and intensity (Lavigne et al. 2013; Vidal et al. 2015) and glaciochemical and petrological evidence (Sigl et al. 2015; Vidal et al. 2016, respectively) for a sulphur yield exceeding that of the Tambora eruption in 1815 (Oppenheimer 2003b) are undisputed, Samalas’ atmospheric and climatic responses have been subject to enduring debated (Stothers 2000; D’Arrigo et al. 2005; Timmreck et al. 2009; Timmreck 2012; Schmidt and Robock 2015; Guillet et al. 2017; Wade et al. 2020). Moreover, climate models, which are sensitive to parametrisation of the amplitude and timing of radiative forcing (Jungclaus et al. 2017), suggest a stronger and more prolonged Northern Hemisphere (NH) summer cooling than that found in tree ring-based temperature reconstructions (Anchukaitis et al. 2012; Mann et al. 2012).

In the absence of a firm historical date or cosmogenic tie-point for radiocarbon dating of relict wood samples (Büntgen et al. 2018; Reinig et al. 2021), fixing the year of an eruption typically combines evidence from ice cores (e.g., sulphur and/or ash), historical sources (e.g., records attesting to stratospheric aerosol veils), proxy-based climate reconstructions, and model-based climate simulations. Yet in the case of Samalas, there is no robust combination of evidence that can discriminate between 1257 or 1258 as the year of eruption. Stothers (2000) argued for early 1258, while Oppenheimer (2003a) and Lavigne et al. (2013) have suggested 1257. Although the earlier date is the more widely quoted, the arguments in favour are not watertight. For instance, they include documentary evidence for mild January conditions in 1258 over France (Lavigne et al. 2013) without mentioning that the association of large volcanic eruptions and NH winter warming (Zambri et al. 2017) is still disputed (Polvani and Camargo 2020).

Even less coherent is the available climate proxy and model evidence of the regional and global hydroclimatic response to sulphur-rich eruptions (Rao et al. 2017; Liu et al. 2020). While models simulate regional-scale precipitation deficits in South, East and Southeast Asia (Schneider et al. 2009; Iles et al. 2013), the sparse annually-resolved reconstructions in those regions show generally weak signals (Anchukaitis et al. 2010). Further complications in untangling the post-volcanic response patterns arise from the dynamical interactions amongst other modes of climate variability on radiative forcing, including the position of the Intertropical Convergence Zone (ITCZ) (Tejedor et al. 2021a), and both the phasing and intensity of the El-Niño Southern Oscillation (ENSO) (Dee et al. 2020; Predybaylo et al. 2020) and Atlantic Multidecadal Oscillation (AMO) (Mann et al. 2021).

The latest generation of annually-resolved and absolutely-dated summer temperature reconstructions for different parts of the NH, either using tree-ring width (TRW) (Büntgen et al. 2020, 2021), maximum latewood density (MXD) (Schneider et al. 2015), or a combination thereof (Stoffel et al. 2015; Wilson et al. 2016), all register distinct cold spells in 1258 and 1259 (Fig. 1a, b). The MXD-based cooling signal recovers rapidly after 1258 (Esper et al. 2015), the TRW-MXD hybrids reveal a 2-year post-eruption cold period, and the TRW-only reconstruction exhibits lowest temperatures in 1259. These differences arguably reflect the slightly ‘redder’ power spectrum in TRW than MXD chronologies (Franke et al. 2013; Büntgen et al. 2020), since the latter are generally less affected by biological memory (Büntgen et al. 2006; Esper et al. 2013). While too much persistence in specific TRW chronologies can be a possible disadvantage for reconstructing the amplitude and duration of climatic extremes (Ljungqvist et al. 2020), there are still too few MXD chronologies spanning the past millennium (Schneider et al. 2015; Esper et al. 2016).

Fig. 1
figure 1

a Five tree ring-based summer temperatures reconstructions that either use ring widths (Büntgen et al. 2021), maximum latewood densities (Schneider et al. 2015), or different combinations thereof (Stoffel et al. 2015; Wilson et al. 2016). Each timeseries has been scaled against the 1961–90 June–August (JJA) average instrumental temperature anomalies over the 30°–70° N extra-tropical landmass. b The five reconstructions exhibit a sharp depression in 1258 followed by parameter-specific recovering rates of up to one decade. c, d Ice core-based estimates of the total annual Stratospheric Aerosol Optical Depth (SAOD) averaged over 30°–90° N (Toohey and Sigl 2017). e A double-stained, wood anatomical thin section of a Mountain pine (Pinus uncinata) from the upper treeline ~ 2200 m asl in the central Spanish Pyrenees (Piermattei et al. 2020), showing a lack of cell wall lignification and frost damage in 1258

The estimated annual, NH extra-tropical, Stratospheric Sulphate Aerosol Optical Depth (SAOD) signature from Samalas occurs in 1258 (Toohey and Sigl 2017). This is consistent with the strongest cooling seen in the tree ring-based reconstructions (Fig. 1c, d). The temperature-SAOD link has been corroborated by independent wood anatomical evidence from the Spanish Pyrenees (Fig. 1e), matching the occurrence of distinct Blue Rings to severe ephemeral cooling during summer 1258 (Büntgen et al. 2017; Piermattei et al. 2020). Considering all the available climate proxies for the NH extra-tropics (Fig. 1), and our general understanding of volcanic forcing of the Earth’s climate system (Gao et al. 2008; Timmreck 2012; Stenchikov 2021), global patterns in the temperature and hydroclimate response to the ~ 1257 Samalas eruption remain unclear. This adds to the challenge of unravelling the direct and indirect interrelationships between climate variability and human history in the aftermath of the event (Guillet et al. 2017; Di Cosmo et al. 2021).

Here, we present a near-global network of 249 TRW chronologies from six continents to reconstruct the temperature and precipitation fingerprints of the Samalas eruption. More specifically, we aim to (1) narrow the time window of the eruption, (2) detect hydroclimatic responses, (3) identify and quantify any response in the Southern Hemisphere (SH), and (4) describe the overall spatiotemporal pattern of climate variability following the eruption. Furthermore, we compare our high-resolution proxy evidence from the thirteenth century against output from state-of-the-art paleoclimate model simulations, discuss agreement and disagreement in the context of natural climate dynamics, and make recommendations for research priorities concerning the volcano-climate interface.

2 Materials and methods

2.1 Tree-ring width chronologies

We compiled raw TRW measurements between 1200 and 1300 from 249 sites on all continents except Antarctica (Fig. 2). Most of the data come from the ITRDB and represent North America (122 sites) and Europe (78 sites). Six sites are in the northern boreal forest of Siberia, and another 18 represent high-elevation treeline ecotones in central Asia. Only one site in Vietnam is in the tropics, ~ 2000 km north of Lombok (Fig. 2 insert photograph). Compared to the 232 sites in the NH, there are just 17 sites in the SH covering the entire thirteenth century, two in New Zealand and Chile, three in Tasmania, and ten in Argentina. With only four sites in Morocco and one in Algeria, Africa is considerably underrepresented. The average site elevation is 1550 m asl, with a maximum range from 0 to 4350 m asl. The lowest elevation trees are oaks in the Netherlands, and highest are junipers on the Tibetan Plateau. Most of the sites are located either below 500 m asl (93 sites) or between 2000 and 4000 m asl (102 sites). The network is composed of 205 conifer and 44 broadleaf species (see Fig. 2 for further details on site and species distribution).

Fig. 2
figure 2

Global distribution of 249 regional tree-ring width (TRW) chronologies (dots), classified into 79 and 170 sites where either precipitation (blue) or temperature (red) is considered the predominant growth control (some dots overlap). The elevational range of the 79 precipitation-sensitive sites is 0–3500 m asl, with a mean of 435 and a median of 180 m asl. The elevational range of the 170 temperature-sensitive sites is 50–4350 m asl, with a mean of 2069 and a median of 2205 m asl. The composition of the 79 precipitation-sensitive sites includes Quercus sp. (40), Cupressaceae (21), Pinus sp. (6), Pinaceae (6), Juniperus sp. (5), and Larix sp. (1). The composition of the 170 temperature-sensitive sites is Pinus sp. (77), Pinaceae (29), Juniperus sp. (26), Cupressaceae (19), Larix sp. (15), and others (4)

After considering each species’ physiological requirements and biogeographic setting at each site, we divided the TRW chronology network into two climate response archetypes, one containing 170 sites where tree growth is predominantly controlled by changes in warm-season temperature and another containing 79 mainly precipitation-limited sites (Fritts 1976), respectively. Though one may argue the growth-climate association of some TRW datasets in this collection would have been better described as a mixed signal (see “Discussion” for details), our binary classification was necessary for comparison with independent model simulations. We therefore followed the general assumption that conifer growth at the upper elevational and latitudinal, species-specific distribution limit is mainly temperature dependent, whereas tree growth at lower elevations under temperate and arid conditions typically depends on precipitation-driven soil moisture availability (see St. George 2014 for a NH perspective).

During the thirteenth century (1200–1300), the average number of raw TRW measurement series in the individual datasets is 15 series, with an average minimum and maximum replication of 11 and 18 series, respectively. To remove effects of the biological age-trend, non-stationary noise and heteroscedasticity from the raw TRW measurements (Cook and Peters 1981), we used individual cubic smoothing spline functions with 50% frequency cut-off at 20 years. This rather conservative standardisation technique emphasizes interannual to decadal growth variability but suppresses any lower frequency information (Cook et al. 1995; Helama et al. 2004; Battipaglia et al. 2010). The high-frequency TRW site chronologies were derived from bi-weight robust means, and their temporal variance was stabilized (Osborn et al. 1997). All detrending and chronology development steps were performed by the latest version of the ARSTAN software (Cook et al. 2017). To further explore chronology robustness, we applied 40-year spline functions and tested for the effects of power-transformation before detrending (Cook and Peters 1997). Since neither the different spline functions nor the power-transformation had significant effects on the resulting TRW chronologies. We decided to proceed with the 20-year spline versions. The resulting chronologies were then normalized over their common period 1200–1300 (i.e., transformed to have a mean of zero and a standard deviation of one). The normalized TRW indices were plotted and mapped in QGIS 3.10.13 (QGIS Development Team 2020). To account for the different phasing of extra-tropical growing seasons between hemispheres (Schulman 1956), all TRW chronologies from the SH were lagged by approximately 6 months (Büntgen and Oppenheimer 2020).

2.2 Climate model simulations

For proxy-model comparison (Braconnot et al. 2012; Phipps et al. 2013; PAGES2K-PMIP3 Group 2015; PAGES Hydro2k Consortium 2017; Zhu et al. 2020; Tejedor et al. 2021b), we used the average simulated June–August (JJA) temperature and precipitation output between 1200 and 1300, from a suite of 12 externally forced Earth System Model simulations of the last millennium (Table 1), generated by the Paleoclimate Modelling Inter-Comparison Project 3 (PMIP3http://www.pmip3.lsce.ipsl.fr; Braconnot et al. 2011; Lohmann et al. 2013). Varying in spatial resolution (Table 1), the individual models used one of three volcanic forcing datasets based on ice core evidence from Greenland and Antarctica (Crowley et al. 2008; Gao et al. 2008; Crowley and Unterman 2013).

For comparison with the individual TRW chronologies, spatial overrepresentation of the coarse output of the individual Earth System models was avoided by cluster analyses. We applied K-means cluster analysis on the latitude, longitude and growth response of all TRW chronologies in the three consecutive years of 1257, 1258 and 1259 to calculate independent spatial clusters for both, temperature and precipitation sensitive TRW sites. We therefore used a pre-defined number of ten clusters to avoid too much spatial heterogeneity, which would be inappropriate to compare with the generally coarse resolution of the climate model simulations that is in the order of hundreds of kilometres. The determination of the lower limit of clusters comes from the analysis of the spatial de-correlation distance of similar climatic characteristics reflected by the TRW sites (Jones et al. 1997). The selection of ten clusters also guaranteed a minimum of three model grid points as a robust basis for the proxy-model comparison. Using K-means cluster analysis of temperature- and precipitation-sensitive TRW chronologies, we calculated for each reginal TRW cluster and each summer between 1257 and 1259 the percentage of models showing the same direction of climate response that was either positive or negative. Although comparing the direction of climate response in both the model and TRW clusters facilitates model assessment at regional scales, this technique ignores differences in the magnitude of physical climate variables (Phipps et al. 2013).

In addition, Spearman’s Rank correlation coefficients were used to assess the strength of the relationship between the reconstructed and simulated summer temperature and precipitation changes for the years 1257, 1258 and 1259 at each individual TRW site (e.g., 170 TRW chronologies for temperature and 79 TRW chronologies for precipitation). All data were normalized over the common period 1200–1300 (mean of zero and standard deviation of one). Due to the spatial resolution of climate model data, different TRW sites fall within the same model grid square and may therefore introduce a reduced degree of freedom owing to spatial autocorrelation (Mets et al. 2017). This has been considered by using the effective degrees of freedom for the proxy and model data when estimating t statistics and significance levels of the Spearman’s Rank correlation coefficients (Table 2). To mitigate methodological constraints, the proxy and model data were also compared at the grid level.

3 Results

3.1 Tree-ring width responses

Growth response patterns of the 170 predominantly temperature-sensitive TRW chronologies in 1257 indicate generally positive anomalies across Eurasia, mixed responses in central Asia, and overall neutral to positive anomalies in the SH (Fig. 3). The first indication of a distinct cold spell is expressed by 24 high-elevation sites in the Sierra Nevada and Rocky Mountains of the western US and two sites in western Canada (Fig. 4). A total of 133 temperature-sensitive sites exhibits a growth decrease in 1258. The evidence for summer cooling in 1258 is most distinct in the mid-latitudes of North America and Eurasia. By contrast, southern sites in central Asia and at high-latitudes in both hemispheres express somewhat warmer than average conditions. In 1259 a total of 117 temperature-sensitive sites exhibits a growth increase from the previous year, and the between-site variation is much reduced. While many of the high-elevation and high-latitude TRW chronologies in North America reveal positive indices, growth depressions persist or become evident in most European, Scandinavian and Siberian chronologies. Sites in central Asia, as well as South America and Australasia show mostly enhanced growth. Comparing histograms of the climate response of all 170 temperature-sensitive TRW chronologies in 1257, 1258 and 1259 highlights 1258 as the most distinctly cold summer (Fig. 3).

Fig. 3
figure 3

Growth response of 170 temperature-sensitive tree-ring width (TRW) chronologies in 1257, 1258 and 1259 (top to bottom). Positive anomalies of the normalized proxy timeseries (with a mean of zero and a standard deviation of one for the period 1200–1300), indicative of warmer summer conditions, are shown in light orange, orange and red; negative anomalies refer to colder summer conditions and range from light to dark blue

Fig. 4
figure 4

Growth response of 26 temperature-sensitive tree-ring width (TRW) chronologies in the western US (#24) and Canada (#2) that reveal negative anomalies already in 1257. All 26 sites are located between 50 and 3660 m asl (with a mean and median of 2247 and 2364 m asl, respectively), and within 34°–52° N and 104°–126° E

Growth response patterns in 1257 of the 79 predominantly precipitation-sensitive TRW chronologies indicate generally positive anomalies across the network, with a few exceptions found in the eastern US and Chile (Fig. 5), for example. The most heterogeneous behaviour is found across Europe and the Mediterranean region. In 1258, 37 precipitation-sensitive sites display a growth decrease while 42 sites show a growth increase from the previous year. The strongest growth depressions are found in the western US, whereas the strongest growth increase is found in Vietnam. Sites in the eastern Mediterranean region, central Asia and the SH all exhibit positive anomalies. Heterogeneous patterns prevail across Europe. In 1259, most of the mid-latitude sites in North America and Eurasia show the opposite response to the previous year: positive anomalies in the western US, but negative values along the east coast of the US and in Pakistan. Again, TRW patterns in the European chronologies remain ambiguous, whereas sites in the eastern Mediterranean continue to show above average growth. Reviewing the spatial response distribution of all 79 precipitation-sensitive TRW chronologies between 1257 and 1259 reveals a similar picture for each year.

Fig. 5
figure 5

Growth response of 79 precipitation-sensitive tree-ring width (TRW) chronologies in 1257, 1258 and 1259 (top to bottom). Negative anomalies of the normalized timeseries (with a mean of zero and a standard deviation of one for the period 1200–1300), indicative of drier summer conditions, are shown in light orange, orange and red, whereas positive anomalies refer to wetter summer conditions and range from light to dark blue

In summary, of the 170 chronologies that contain a predominant temperature signal, the first indication of a cold spell in 1257 originates from distinct growth depressions at 26 high-elevation sites in the western US and Canada (Figs. 3, 4). Widespread summer cooling across the mid-latitudes in North America and Eurasia is, however, most pronounced in 1258, and many boreal records in Scandinavia and Siberia show negative growth anomalies in 1259. In contrast to the marked post-Samalas temperature response in the NH, there is no dramatic hydroclimatic fingerprint exhibited by the 79 precipitation-sensitive chronologies.

A more detailed assessment reveals a negative relationship between the growth response of the 170 temperature-sensitive TRW chronologies and their site elevation (Fig. 6a). The correlation between lower growth rates at higher elevations is most significant in 1258 compared to the preceding and following years. An inverse, though less distinct, association is exhibited by the 79 precipitation-sensitive TRW chronologies (Fig. 6b), which indicate overall positive trends between increasing growth rates with increasing site elevation. The distribution of the growth response of the 170 temperature-sensitive TRW chronologies in 1257 and 1259 resembles the overall pattern of the thirteenth century (Fig. 7a), whereas the distribution in 1258 is skewed towards negative anomalies. The two strongest growth decreases in 1258 are associated with two high-elevation pine sites in the western US (~ 3535 m asl, ~ 39° N and ~ 104° W). The strongest growth increase in 1259 is exhibited by a pine site at slightly lower elevation in the western US (Mount Goliath; ~ 1900 m asl, ~ 40° N and ~ 108° W). The distribution of the growth response of the 79 precipitation-sensitive TRW chronologies in 1257, 1258 and 1259 resembles the overall pattern of the thirteenth century (Fig. 7b). The only outstanding value is the anomalous growth increase in 1258 found in a cypress chronology from Vietnam (Bidoup Nui Ba National Park; ~ 1700 m asl, ~ 12° N and ~ 108° E).

Fig. 6
figure 6

a Relationship between the growth response and site elevation of 170 temperature-sensitive tree-ring width (TRW) chronologies in 1257, 1258 and 1259 (left to right). b Relationship between the growth response and site elevation of 79 precipitation-sensitive TRW chronologies in 1257, 1258 and 1259 (left to right). Z-scores were calculated relative to the period 1200–1300

Fig. 7
figure 7

a Distribution of the growth response of 170 temperature-sensitive tree-ring width (TRW) chronologies in 1257, 1258 and 1259 (top to bottom), relative to the full thirteenth century distribution (black line). b Distribution of the growth response of 79 precipitation-sensitive TRW chronologies in 1257, 1258 and 1259 (top to bottom), relative to their general distribution in the thirteenth century (black line)

3.2 Proxy-model comparisons

Relative to the thirteenth century mean, eight models produce their lowest JJA temperature anomalies in 1258, and three models in 1259. The lowest JJA temperature anomaly in the FGOALS-gl simulation is in 1276. Spearman’s Rank correlation coefficients between the temperature-sensitive TRW chronologies and the spatially corresponding output from 12 climate model simulations of regional summer temperature anomalies in 1257, 1258 and 1259 reveal four significantly (p < 0.05) positive values in 1258 (Table 2). This finding suggests a high degree of proxy-model coherency on regional-to continental scales due to large volcanic forcing. In 1257, the overall positive TRW temperature response in eight spatial clusters resembles most of the model output (Fig. 8a). Similar response agreements are found in 83% of the comparisons in central Europe and central Asia, in 67% of the comparisons in Alaska and North America, in 58% of the comparisons in north-central Siberia, Australasia and South America, and in 8% of the comparisons in northern Scandinavia. In 1258, the strongest proxy-model agreement is found in North America and central Asia (both 92%), as well as in central Europe (83%). However, poor model agreement with positive TRW anomalies is characteristic for the high-northern latitudes and the SH. In 1259, only the TRW chronologies from central Europe and north-central Siberia exhibit marked growth depressions, which are in line with 100% and 92% of the corresponding model data, respectively. Where TRW is positive, 25% or less of the clustered model output agrees with the reconstructed direction of response.

Fig. 8
figure 8

a Directional agreement between temperature-sensitive tree-ring width (TRW) chronologies in eight spatial clusters and corresponding output from 12 model simulations of regional summer temperature anomalies in 1257, 1258 and 1259 (see Table 1 for model details). The eight temperature clusters represent parts of Alaska (Alaska), North America (N America), central Europe (C Europe), northern Scandinavia (N Scan), north-central Siberia (NC Siberia), central Asia (C Asia), Australasia (Australasia) and South America (S America), and contain 7, 90, 16, 14, 5, 20, 4 and 11 TRW chronologies, respectively. b Directional agreement between precipitation-sensitive tree-ring width (TRW) chronologies in five spatial clusters and their corresponding output from 12 model simulations of regional summer precipitation anomalies in 1257, 1258 and 1259. The five precipitation clusters represent parts of the western U (W US), the eastern US (E US), western Europe (W Europe), eastern Europe (E Europe) and the eastern Mediterranean (E Med), and contain 4, 17, 25, 10 and 14 TRW chronologies, respectively

In contrast to the proxy-model associations of the eight temperature-sensitive TRW clusters, there are only five precipitation-sensitive TRW clusters for which sample size is sufficient to compare the directions of response (Fig. 8b). In 1257, positive anomalies in the precipitation-sensitive TRW chronologies from the western US, western Europe and the eastern Mediterranean are reflected in 50, 42 and 50% of the models, respectively. Model agreement with negative growth responses in the eastern US and eastern Europe are 25 and 17%, respectively. In 1258, positive anomalies in the precipitation-sensitive TRW chronologies in the eastern US, western Europe and the eastern Mediterranean are seen in 8, 42 and 67% of the models, respectively. Model agreement with marked negative growth responses in the western US and eastern Europe is 75%. In 1259, positive anomalies in the precipitation-sensitive TRW chronologies from the western US, western Europe and the eastern Mediterranean are found in 58, 58 and 75% of the models, respectively. Model agreement with distinct negative growth responses in the western US and eastern Europe is 67%. In summary, the spatially explicit TRW response in the NH provides a reasonable benchmark for model evaluation after volcanic forcing.

4 Discussion

Though our network contains 249 TRW chronologies from six continents, it remains spatially biased towards the western US and central Europe. Except for Algeria and Morocco, no data come from Africa, and there are only 17 sites from South America and Australasia that represent the SH (Neukom and Gergis 2011). In addition to the spatial limitations, only 79 TRW chronologies are considered precipitation sensitive, and many of the lower-elevation, mid-latitude sites in North America and Eurasia likely reflect mixed climatic signals (Galván et al. 2014; Hellmann et al. 2016). Chronologies from such sites are particularly prone to nonlinear growth responses under climatic extremes (Vitali et al. 2017; Camarero et al. 2021).

High-elevation sites in the western US and Canada show a strong growth decline in 1257 that may reflect the impact of radiative forcing caused by Samalas aerosol (Fig. 4). If this deflection is an indication of post-volcanic boreal summer cooling, it would suggest an eruption at least weeks before the end of the growing season of these trees, i.e., before around August 1257, which would be consistent with the eruption window suggested by Oppenheimer (2003a) and Lavigne et al. (2013). Our network, however, reveals the strongest cooling in 1258 across much of the NH mid-latitudes, consistent with previously published regional evidence (Büntgen et al. 2017; Guillet et al. 2017; Piermattei et al. 2020). The strongest growth depression across the northern boreal forest of Scandinavia and Siberia in 1259 is consistent with previous reports of delayed northern Eurasia summer cooling 2 years after volcanic eruptions (Esper et al. 2015), as well as with documentary evidence from the Altai (Guillet et al. 2017).

The few TRW chronologies from South America and Australasia mainly indicate warmer conditions between 1257 and 1259, which challenges the argument of a globally consistent cooling response to the Samalas eruption (Neukom et al. 2019). Warmer SH summers in 1258 and 1259 might be related to El Niño conditions and a positive phase of the Pacific North America pattern, as evidenced by increased growth rates in Alaska (Rohli and Vega 2018). The observed response pattern at North American sites contributes to the debate around connections between volcanism and ENSO (Mann et al. 2005; Stevenson et al. 2018; Dee et al. 2020; Robock 2020), and more generally on the possible linkages between climate oscillations and volcanic eruptions (Mann et al. 2021). In addition, interactions between changes in external climate forcing and direct climate effects on temperature and precipitation over the SH are considerably influenced by the internal state of background variability that is affected by the large oceanic coverage of the SH. The amount of sea-ice around Antarctica potentially also impacts the state of the Antarctic Oscillation and associated synoptic-scale climate variability in the aftermath of volcanic eruptions. This somewhat deviating and counter-intuitive response of SH climate to volcanic eruptions is further reflected in the spatial extent and temporal structure of anomalous cold periods during the last millennium (PAGES 2k Consortium 2013), such as the Little Ice Age.

A more complex picture is drawn by the 79 predominantly precipitation-sensitive TRW chronologies, which show a high level of spatiotemporal heterogeneity between 1257 and 1259. The chronologies from the eastern US indicate summer drying in 1257, whereas the central US sites show a mixed response in this year. Positive precipitation anomalies across the US east coast in 1258 coincide with summer drying in the southwest, which is generally in line with the post-eruptive response patterns simulated by the HadCM3 and CMIP5 climate models (Iles et al. 2013; Iles and Hegerl 2014) and independent proxy evidence (Tejedor et al. 2021a). The inverse response pattern in 1259 suggests North American internal climate variability played an important role in modulating precipitation regimes in the aftermath of the Samalas eruption (Tejedor et al. 2021a; Stevenson et al. 2018). The European chronologies indicate dryer northern and wetter southern European summers in 1259, which is in line with previous findings (Rao et al. 2017; Di Cosmo et al. 2021). This dipole pattern possibly reflects a negative phase of the East Atlantic Pattern (EAP), or a negative phase of the summer NAO (Liu et al. 2020). However, other studies did not find any dramatic central European hydroclimatic response to the Samalas eruption (Raible et al. 2016; Büntgen et al. 2021; Tejedor et al. 2021a). Though there is some evidence for wetter summers following large tropical eruptions in southeast Asia (Anchukaitis et al. 2010; Buckley et al. 2010), the overall limited data availability from Asia hinders any further resolution of the relationship, if any. Moreover, the enhanced NH summer cooling after the Samalas eruption may have shifted the ITCZ southward, resulting in wetter conditions in South America and Australasia (Tejedor et al. 2021a; Liu et al. 2016), though we reiterate that the paucity of TRW sites in the SH precludes definitive conclusions. Nevertheless, the precipitation-sensitive TRW chronologies in South America and Australia indicate wetter conditions in 1258, which is in line with previous findings (Tejedor et al. 2021a).

Depending on the volcanic forcing used (Table 1), most PMIP3 models have their lowest NH surface temperatures for the thirteenth century in 1258 (Crowley and Unterman 2013), whereas CCSM4, IPSL-CM5A-LR and BCC-CSM1.1 indicate their coldest summers in 1259. The only model that uses the Gao et al. (2008) forcing and produces its lowest value in 1258 is GISS-E2-Rp122. Based on an earlier volcanic forcing dataset (Crowley et al. 2008), FGOALS-gl simulates the lowest summer temperature in 1257. In addition to uncertainty in the timing of the eruption, which hinders the comparison with independent proxy data (Büntgen and Oppenheimer 2020), different background conditions and forcing datasets also affect the simulated climate response to the Samalas eruption (Lohmann et al. 2013; Schmidt et al. 2011). The overall weak proxy-model association may be best explained by a combination of parameterisation uncertainty in the climate models (Gettleman and Rood 2016; Diaz and Vera 2018), aerosol lifetime differences and microphysical processes that are variably represented in the models (Timmreck et al. 2009; Timmreck 2012), and mixed signals in the TRW chronologies (Fritts 1976).

Table 1 Overview of the 12 externally forced climate model simulations of the last millennium, generated by the Paleoclimate Modelling Inter-Comparison Project 3 (PMIP3; Braconnot et al. 2011) and used in this study for comparison against the dendrochronological evidence
Table 2 Spearman's Rank correlation coefficients and their corrected two-tailed significance for all possible pairings between (A) temperature and (B) precipitation sensitive TRW chronologies and the 12 externally forced climate model simulations

5 Conclusions

We have presented a near global network of 249 annually-resolved and absolutely-dated regional TRW chronologies that reflect high-frequency summer temperature or precipitation variability, span the period of the Samalas eruption, and capture its effects on surface climate. Our analysis contributes in several ways to understanding the timing and climatic impacts of the eruption: (1) We provide tree-ring evidence for a possible eruption date of Samalas between May and July of 1257, followed by a distinct cold spell across the mid-latitudes of the NH in 1258, delayed cooling at the high-northern latitudes in 1259, and some indication for generally warmer growing seasons in the SH between 1257 and 1259. (2) We do not identify any dramatic hydroclimate tree-ring responses to the Samalas eruption. (3) We find no tree-ring evidence for globally coherent and exceptionally colder or wetter summers in the aftermath of the Samalas eruption.

Due to spatial data biases and mixed climate signals, future research would profit from focussing on the development of multi-centennial to millennial-long tree-ring chronologies in those regions of both hemispheres that are still underrepresented. Such endeavours should include wood anatomical, density and isotopic measurements, which contain less biological memory and may represent different seasonal windows compared to the classical TRW data.