Siberian tree-ring and stable isotope proxies as indicators of temperature and moisture changes after major stratospheric volcanic eruptions

Stratospheric volcanic eruptions have far-reaching impacts on global climate and society. Tree rings can provide valuable climatic information on these impacts across different spatial and temporal scales. To detect temperature and hydroclimatic changes after strong stratospheric Common Era (CE) volcanic eruptions for the last 1500 years (535 CE unknown, 540 CE unknown, 1257 CE Samalas, 1640 CE Parker, 1815 CE Tambora, and 1991 CE Pinatubo), we measured and analyzed tree-ring width (TRW), maximum latewood density (MXD), cell wall thickness (CWT), and δ¹³C and δ¹⁸O in tree-ring cellulose chronologies of climate-sensitive larch trees from three different Siberian regions (northeastern Yakutia – YAK, eastern Taimyr – TAY, and Russian Altai – ALT).

All tree-ring proxies proved to encode a significant and specific climatic signal of the growing season. Our findings suggest that TRW, MXD, and CWT show strong negative summer air temperature anomalies in 536, 541–542, and 1258–1259 at all studied regions. Based on δ¹³C, 536 was extremely humid at YAK, as was 537–538 in TAY. No extreme hydroclimatic anomalies occurred in Siberia after the volcanic eruptions in 1640, 1815, and 1991, except for 1817 at ALT. The signal stored in δ¹⁸O indicated significantly lower summer sunshine duration in 542 and 1258–1259 at YAK and 536 at ALT. Our results show that trees growing at YAK and ALT mainly responded the first year after the eruptions, whereas at TAY, the growth response occurred after 2 years.

The fact that differences exist in climate responses to volcanic eruptions – both in space and time – underlines the added value of a multiple tree-ring proxy assessment. As such, the various indicators used clearly help to provide a more realistic picture of the impact of volcanic eruption on past climate dynamics, which is fundamental for an improved understanding of climate dynamics, but also for the validation of global climate models.

1 Introduction

Major stratospheric volcanic eruptions can modify the Earth's radiative balance and substantially cool the troposphere. This is due to the massive injection of sulfate aerosols, which reduce surface temperatures on timescales ranging from months to years (Robock, 2000). Volcanic aerosols significantly absorb terrestrial radiation and scatter incoming solar radiation, resulting in a cooling that has been estimated to about 0.5 ^∘C during the 2 years following the Mount Pinatubo eruption in June 1991 (Hansen et al., 1996).

Since trees – as living organisms – are impacted in their metabolism by environmental changes, their responses to these changes are recorded in the biomass, as is found in tree-ring parameters (Schweingruber, 1996). The decoding of tree-ring archives is used to reconstruct past climates. A summer cooling of the Northern Hemisphere ranging from 0.6 to 1.3 ^∘C has been reported after the strongest known volcanic eruptions of the past 1500 years (1257 CE Samalas, 1815 Tambora, and 1991 Pinatubo) based on temperature reconstructions using tree-ring width (TRW) and maximum latewood density (MXD) records (Briffa et al., 1998; Schneider et al., 2015; Stoffel et al., 2015; Wilson et al., 2016; Esper et al., 2017, 2018; Guillet et al., 2017; Barinov et al., 2018).

Climate simulations show significant changes in the precipitation regime after large volcanic eruptions. These include, among others, rainfall deficit in monsoon-prone regions and in southern Europe (Joseph and Zeng, 2011), as well as wetter than normal conditions in northern Europe (Robock and Liu 1994; Gillet et al., 2004; Peng et al., 2009; Meronen et al., 2012; Iles et al., 2013; Wegmann et al., 2014). However, despite recent advances in the field, the impacts of stratospheric volcanic eruptions on hydroclimatic variability at regional scales remain largely unknown. Therefore, further knowledge about moisture anomalies is critically needed, especially at high-latitude sites where tree growth is mainly limited by summer temperatures.

As dust and aerosol particles from large volcanic eruptions primarily affect the radiation regime, three major drivers of plant growth (i.e., photosynthetically active radiation – PAR, temperature, and vapor pressure deficit – VPD) will be affected by volcanic activity. This is reflected in low TRW as a result of reduced photosynthesis but even more so due to low temperature. As cell division is temperature dependent, its rate (tree-ring growth) will exponentially decrease with decreasing temperature below +3 ^∘C (Körner, 2015), outweighing the “low light–low photosynthesis” effect by far.

Furthermore, over the last years some studies using mainly carbon isotopic signals (δ¹³C) in tree rings showed ecophysiological responses of trees to volcanic eruptions at the middle (Battipaglia et al., 2007) or high (Gennaretti et al., 2017) latitudes. By contrast, a combination of both carbon (δ¹³C) and oxygen (δ¹⁸O) isotopes in tree rings has been employed only rarely to trace volcanic eruptions in high-latitude or high-altitude proxy records (Churakova (Sidorova) et al., 2014).

The application of TRW, MXD, and cell wall thickness (CWT), as well as δ¹³C and δ¹⁸O, in tree cellulose chronologies is a promising tool to disentangle hydroclimatic variability as well as winter and early spring temperatures at high-latitude and high-altitude sites (Kirdyanov et al., 2008; Sidorova et al., 2008, 2010, 2011; Churakova (Sidorova) et al., 2014; Castagneri et al., 2017). In that sense, recent CWT measurements allowed for the generation of high-resolution, seasonal information on water and carbon limitations on growth during spring and summer (Panyushkina et al., 2003; Sidorova et al., 2011; Fonti et al., 2013; Bryukhanova et al., 2015). Depending on site conditions, δ¹³C variations reflect light (stand density) (Loader et al., 2013) and water availability (soil properties) and air humidity (proximity to open waters, i.e., rivers, lakes, swamps, and orography) as these parameters have been recognized to modulate the stomatal conductance (g_l) controlling carbon isotopic discrimination.

Depending on the study site, a decrease in the carbon isotope ratio can be expected after stratospheric volcanic eruptions due to limited photosynthetic activity and higher stomatal conductance, which in turn would be the result of decreased temperatures, VPD, and a reduction in light intensity. By contrast, volcanic eruptions have also been credited for an increase in photosynthesis as dust and aerosol particles cause increased light scattering, compensating for the light reduction (Gu et al., 2003). A significant increase in δ¹³C values in tree-ring cellulose should be interpreted as an indicator of drought (stomatal closure) or high photosynthesis (Farquhar et al., 1982). In the past, very little attention has been paid to the elemental and isotopic composition of tree rings for years during which they may have been subjected to the climatic influence of powerful, but remote and often tropical, volcanic eruptions.

Table 1 List of stratospheric volcanic eruptions used in the study.

NA – not available.

In this study, we aim to fill this gap by investigating the response of different components of the Siberian climate system (i.e., temperature, precipitations, VPD, and sunshine duration) to stratospheric volcanic events of the last 1500 years. By doing so, we seek to extend our understanding of the effects of volcanic eruptions on climate by combining multiple climate-sensitive variables measured in tree rings that were clustered around the time of the major volcanic eruptions (Table 1). We focus our investigation on remote tree-ring sites in Siberia, two at high latitudes (northeastern Yakutia – YAK and eastern Taimyr – TAY) and one at high altitude (Russian Altai – ALT), for which long tree-ring chronologies were previously developed with highly climate-sensitive trees. We assemble a dataset from five tree-ring proxies: TRW, MXD, CWT, δ¹³C, and δ¹⁸O in larch tree-ring cellulose chronologies. This is done in order to (1) determine the major climatic drivers of the tree-ring proxies, and to evaluate their individual and integrative response to climate change, and to (2) reconstruct the climatic impacts of volcanic eruptions over specific periods of the past (Table 1).

2 Material and methods

2.1 Study sites

The study sites are situated in Siberia (Russian Federation), far away from industrial centers (and 1500–3400 km apart from each other), in zones of continuous permafrost in northeastern Yakutia (YAK: 69^∘ N, 148^∘ E) and eastern Taimyr (TAY: 70^∘ N, 103^∘ E) and mountain permafrost in Altai (ALT: 50^∘ N, 89^∘ E) (Fig. 1a, Table 2). Tree-ring samples were collected during several field trips and included old relict wood and living larch trees: Larix cajanderi Mayr (up to 1216 years) in YAK, Larix gmelinii Rupr. (max. 640 years), in TAY and Larix sibirica Ldb. (max. 950 years) in ALT. TRW chronologies have been developed and published in the past (Fig. 1; Hughes et al., 1999; Sidorova and Naurzbaev, 2002; Sidorova, 2003, for YAK; Naurzbaev et al., 2002; Panyushkina et al., 2003, for TAY; Myglan et al., 2008, for ALT).

Due to the remote location of our study sites, we used meteorological data from monitored weather stations located at distances ranging from 50 to 200 km from the sampled sites. Temperature data from these weather stations are significantly correlated (r>0.91; p<0.05) with gridded data (http://climexp.knmi.nl, last access: January 2018). However, poor correlation is found with precipitation data (r<0.45; p<0.05), which is most likely the result of local topography (Churakova (Sidorova) et al., 2016).

Mean annual air temperature is lower at the high-latitude YAK and TAY sites than at the high-altitude ALT site (Table 2). Annual precipitation is low (153–269 mm yr⁻¹) at all sites. The growing season calculated with the tree growth threshold of +5 ^∘C (Fritts, 1976; Schweingruber, 1996) is very short (50–120 days) at all locations (Table 2). Sunshine duration is higher at YAK and TAY (ca. 18–20 h day⁻¹ in summer) compared to ALT (ca. 18 h day⁻¹ in summer) (Sidorova et al., 2005, 2011; Myglan et al., 2008; Churakova (Sidorova) et al., 2014).

Figure 1 Location of the study sites (stars) and known volcanos from the tropics (black dots) considered in this study (a). Annual tree-ring width index (light lines) and smoothed by 51-year Hamming window (bold lines) from northeastern Yakutia (YAK – blue, b) (Hughes et al., 1999; Sidorova and Naurzbaev, 2002; Sidorova, 2003), eastern Taimyr (TAY – green, c) (Naurzbaev et al., 2002), and Russian Altai (ALT – red, d) (Myglan et al., 2009). Photos show the larch stands at the YAK, TAY (Mukhtar M. Naurzbaev), and ALT (Vladimir S. Myglan) sites.

Table 2 Tree-ring sites in northeastern Yakutia (YAK), eastern Taimyr (TAY), and Altai (ALT) and weather stations used in the study. Monthly air temperature (T, ^∘C), precipitation (P, mm), sunshine duration (S, h month⁻¹), and vapor pressure deficit (VPD, kPa) data were downloaded from the meteorological database: http://aisori.meteo.ru/ClimateR (last access: January 2018).

^a Abaimov et al. (1996); Hughes et al. (1999); Churakova (Sidorova) et al. (2016); ^b Naurzbaev et al. (2002); ^c Sidorova et al. (2011).

2.2 Selection of volcanic events and larch subsamples

Identification of the events used in this study was based on volcanic aerosols deposited in ice core records (Zielinski, 1994; Robock, 2000) and, more precisely, on Toohey and Sigl (2017), in which the authors listed the top 20 eruptions over the past 2000 years based on volcanic stratospheric sulfur injection (VSSI). From that list, we selected reconstructed VSSI and events that are well recorded in tree-ring proxies and may thus have had a noticeable impact on the forest ecosystems in high-latitude and high-altitude regions (Briffa et al., 1998; D'Arrigo et al., 2001; Churakova (Sidorova) et al., 2014; Büntgen et al., 2016; Gennaretti et al., 2017; Helama et al., 2018). Therefore, based on our previously published TRW and developed MXD, CWT, δ¹³C, and δ¹⁸O in tree-ring cellulose chronologies, we selected the periods 520–560, 1242–1286, 1625–1660, 1790–1835, and 1950–2000 CE with strong volcanic eruptions in 535, 540, 1257, 1640, 1815, and 1991 CE, as they have had far-reaching climatic effects (Table 1). The recent period 1950–2000 is used to calibrate the tree-ring proxy against available climate data.

Tree-ring material was prepared from the 2000-year TRW chronologies available at each site from previous studies (Fig. 1b–d). According to the level of conservation of the material, the largest possible number of samples was prepared for each of the proxies. Unlike TRW, which could be measured on virtually all samples, some of the material was not available with sufficient quality to allow for tree-ring anatomy and stable isotope analysis. We therefore use a smaller sample size for CWT (n=4) and stable isotopes (n=4) than for TRW (n=12) or MXD (n=12). Nonetheless, replications are still comparable with those used in reference papers on stable isotopes and CWT (Loader et al., 1997; Panyushkina et al., 2003).

2.3 Tree-ring width analysis

The ring width of 12 trees was remeasured for each selected period. Cross-dating was checked by comparison with the existing full-length 2000-year TRW chronologies (Fig. 1). The TRW series were standardized using the ARSTAN program (Cook and Krusic, 2008) with a negative exponential curve (k>0) or a linear regression (any slope) prior to bi-weight robust averaging (Cook and Kairiukstis, 1990). Signal strength in the regional TRW chronologies was assessed with the expressed population signal (EPS) statistics as it measures how well the finite sample chronology compares with a theoretical population chronology with an infinite number of trees (Wigley et al., 1984). Mean inter-series correlation (RBAR) and EPS values of stable isotope chronologies were calculated for the period 1950–2000, for which individual trees were analyzed separately. All series have RBAR ranges between 0.59 and 0.87, and the common signal exceeds the EPS threshold of 0.85. Before 1950, we used pooled cellulose only. For all other tree-ring parameters and studied periods, the EPS exceeds the threshold of 0.85, and RBAR values range from 0.63 to 0.94.

2.4 Image analysis of cell wall thickness (CWT)

Analysis of wood anatomy was performed for all studied periods with an AxioVision scanner (Carl Zeiss, Germany). Micro-sections were prepared using a sliding microtome and stained with methyl blue (Furst, 1979). Tracheids in each tree ring were measured along five radial files of cells (Munro et al., 1996; Vaganov et al., 2006) selected for their larger tangential cell diameter (T). For each tracheid, CWT was computed separately. In a second step, tracheid anatomical parameters were averaged for each tree ring. Site chronologies are presented for the complete annual ring chronology without standardization due to the lack of a low-frequency trend. CWT data from ALT for the periods 1790–1835 and 1950–2000 were used from past studies (Sidorova et al., 2011; Fonti et al., 2013) and for YAK for the period from 1600 to 1980 from Panyushkina et al. (2003). Unfortunately, the remaining sample material for the 536 CE ring at TAY was insufficient to produce a clear signal. As a result, CWT is missing for 536 CE at TAY (Fig. 2).

2.5 Maximum latewood density (MXD)

Maximum latewood density chronologies from ALT were available continuously for the period 600–2007 CE from Schneider et al. (2015) and Aleksander  V. Kirdyanov (personal communication, 2012), as well as from YAK and TAY for the period 1790–2004 CE from Sidorova et al. (2010). For any other periods, at least six cross sections and for 520–560 CE four sections are used. The wood is subsampled with a double-bladed saw at 1.2 mm thickness with the angle to the fiber direction. The samples were exposed to X-rays for 35–60 min (Schweingruber, 1996). MXD measurements were obtained at 0.01 mm resolution and brightness variations calibrated to g cm⁻³ (Lenz et al., 1976; Eschbach et al., 1995) using a Walesch X-ray densitometer 2003. All MXD series were detrended in the ARSTAN program by calculating deviation from the straight-line function (Fritts, 1976). Site MXD chronologies were developed for each volcanic period using the bi-weight robust averaging.

2.6 Stable carbon (δ¹³C) and oxygen (δ¹⁸O) isotopes in tree-ring cellulose

During photosynthetic CO₂ assimilation ¹³CO₂ is discriminated against ¹²CO₂, leaving the newly produced assimilates depleted of ¹³C. The carbon isotope discrimination (¹³Δ) is partitioned in the diffusional component with a=4.4 ‰ and the biochemical fractionation with b=27 ‰ for C₃ plants during carboxylation via Rubisco. The ¹³Δ is directly proportional to the c_i∕c_a ratio, where c_i is the leaf intercellular, and c_a the ambient CO₂ concentration. This ratio reflects the balance between stomatal conductance (g_l) and photosynthetic rate (A_N). A decrease in g_l at a given A_N results in a decrease in ¹³Δ as c_i∕c_a decreases and vice versa. The same is true when A_N increases or decreases at a given g_l. Since CO₂ and H₂O gas exchange is strongly interlinked with C-isotope fractionation, ¹³Δ is controlled by the same environmental variables, i.e., PAR, CO₂, VPD, and temperature (Farquhar et al., 1982, 1989; Cernusak et al., 2013). The oxygen isotopic compositions of tree-ring cellulose record the δ¹⁸O of the source water derived from precipitation, which itself is related to temperature variations at middle and high latitudes (Craig, 1961; Dansgaard, 1964). It is modulated by evaporation at the soil surface and to a larger degree by evaporative and diffusion processes in leaves; the process is largely controlled by the vapor pressure deficit (Dongmann et al., 1972; Farquhar and Loyd, 1993; Cernusak et al., 2016). A further step of fractionation occurs as sugar molecules are transferred to the locations of growth (Roden et al., 2000). During the formation of organic compounds biosynthetic fractionation leads to a positive shift of the δ¹⁸O values by 27 ‰ relative to the leaf water (Sternberg, 2009). The oxygen isotope variation in tree-ring cellulose therefore reflects mixed climate information, often dominated by a temperature, source water, or sunshine duration modulated by the VPD influence.

The cross sections of relict wood and cores from living trees used for the TRW, MXD, and CWT measurements were then selected for the isotope analyses. We analyzed four subsamples for each studied period according to the standards and criteria described in Loader et al. (2013). The first 50 years of each sample were excluded to limit juvenile effects (McCarroll and Loader, 2004). After splitting annual rings with a scalpel, the whole wood samples were enclosed in filter bags. α-cellulose extraction was performed according to the method described by Boettger et al. (2007). For analyses of the ¹³C∕¹²C and ¹⁸O∕¹⁶O isotope ratios, 0.2–0.3 and 0.5–0.6 mg of cellulose were weighed for each annual ring into tin and silver capsules, respectively. Carbon and oxygen isotopic ratios in cellulose were determined with an isotope ratio mass spectrometer (Delta-S; Finnigan MAT, Bremen, Germany) linked to two elemental analyzers (EA-1108 and EA-1110; Carlo Erba, Italy) via a variable open split interface (CONFLO-II; Finnigan MAT, Bremen, Germany). The ¹³C∕¹²C ratio was determined separately by combustion under oxygen excess at a reactor temperature of 1020 ^∘C. Samples for ¹⁸O∕¹⁶O ratio measurements were pyrolyzed to CO at 1080 ^∘C (Saurer et al., 1998). The instrument was operated in the continuous-flow mode for both the C and O isotopes. The isotopic values were expressed in delta notation multiplied by 1000 relative to the international standards (Eq. 1):

$$\begin{array}{}\text{(1)}& \mathit{\delta }\text{sample}=R_{\text{sample}}/R_{\text{standard}}-\mathrm{1},\end{array}$$

where R_sample is the molar fraction of the ¹³C∕¹²C or ¹⁸O∕¹⁶O ratio of the sample and R_standard the molar fraction of the standards, Vienna Pee Dee Belemnite (VPDB) for carbon and Vienna Standard Mean Ocean Water (VSMOW) for oxygen. The precision is σ±0.1 ‰ for carbon and σ±0.2 ‰ for oxygen. To remove the atmospheric δ¹³C trend after 1800 CE from the carbon isotope values in tree rings (i.e., the Suess effect, due to fossil fuel combustion) we used atmospheric δ¹³C data from Francey et al. (1999; http://www.cmdl.noaa.gov./info/ftpdata.html, last access: January 2018). These corrected series were used for all statistical analyses. The δ¹⁸O cellulose series were not detrended.

2.7 Climatic data

Meteorological series were obtained from local weather stations close to the study sites and used for the computation of correlation functions between tree-ring proxies and monthly climatic parameters (Table 2; http://aisori.meteo.ru/ClimateR, last access: January 2018).

2.8 Statistical analysis

All chronologies for each period were normalized to z scores (Fig. 2). To assess post-volcanic climate variability, we used superposed epoch analysis (SEA; Panofsky and Brier, 1958) with the five proxy chronologies available at each of the three study sites. In this study, intervals of 15 years before and 20 years after a volcanic eruption were analyzed. SEA is applied to the six annually dated volcanic eruptions (Table 1).

To test the sensitivity of the studied tree-ring parameters to climate, bootstrap correlation functions were computed between proxy chronologies and monthly climate predictors using the “bootRes” package of R software (R Core Team, 2016) for the period 1950 (1966)–2000.

To estimate whether volcanic years can be considered as extreme and how anomalous they are compared to non-volcanic years, we computed probability density functions (PDFs; Stirzaker, 2003) for each study site and for each tree-ring parameter over a period of 219 years for which measurements are available (Fig. S1). A year is considered (very) extreme if the value of a given parameter is below the (5th) 10th percentile of the PDF.

Figure 2 Normalized (z-score) individual tree-ring index chronologies (TRW, black), maximum latewood density (MXD, purple), cell wall thickness (CWT, green), δ¹³C (red), and δ¹⁸O (blue) in tree-ring cellulose chronologies from northeastern Yakutia (YAK), eastern Taimyr (TAY), and Altai (ALT) for the specific periods 520–560, 1242–1286, 1625–1660, 1790–1835, and 1950–2000 CE before and after the eruptions in 535, 540, 1257, 1640, 1815, and 1991 CE are presented. Vertical lines show the year of the eruptions.

3 Results

3.1 Anomalies in tree-ring proxy chronologies after stratospheric volcanic eruptions

Normalized TRW chronologies show negative deviations the year following the eruptions at all studied sites (Fig. 2). Regarding CWT, a strong decrease is observed in 537 CE at all study sites. Only two layers of cells were formed in 537 CE (−1.8σ) and 541 (−2.4σ) for YAK compared to the 11–20 layers of cells formed on average during “normal” years. In addition, we also observe the formation of frost rings in ALT between 536 and 538 CE, as well as in 1259. An abrupt CWT decrease is recorded in TAY in 537 (−3.1σ).

Furthermore, we found decreasing MXD values at ALT (−4.4σ) in 537 CE and YAK (−2.8σ) in 536 CE. However, for TAY, the data show a less pronounced pattern of MXD variation (Fig. 2). In this regard, the sharpest decrease was observed in the CWT chronologies from YAK in 540 CE (−1.9σ) and 541 (−2.4σ), whereas the response was smaller in TAY and ALT for the same years (Fig. 2). The ALT δ¹⁸O chronology recorded a drastic decrease in 536 CE with −4.8σ (Figs. 2 and S1). A δ¹⁸O decrease for YAK was found after the 1257 CE Samalas eruption in 1258 CE (−1.5σ) and in 1259 (−2.9σ), which is opposite to the increased δ¹⁸O value found in 1259 CE at ALT (Figs. 2 and S1).

Regarding the carbon isotope ratio, negative anomalies are already observed in ALT in 1258 (−2.3σ). The 540 CE eruption was less clearly recorded in tree-ring proxies from TAY compared to YAK and ALT (Fig. 2). With respect to the 1257 CE Samalas eruption (Fig. 2), the year following the eruption was recorded as very extreme in the TRW, MXD, δ¹⁸O, while it is less extreme in CWT and δ¹³C from YAK. ALT chronologies show a synchronous decrease for all proxies in the 2 years after the eruption (Figs. 2 and S1).

The impacts of the more recent 1640 CE Parker, 1815 Tambora, and 1991 Pinatubo eruptions are, by contrast, far less obvious. In 1642 CE, decreasing values are observed in all tree-ring proxies from the high-latitude sites YAK and TAY, whereas tree-ring proxies are not clearly affected at ALT (Figs. 2 and S1).

Hardly any strong anomalies are observed in 1816 CE in Siberia regardless of the site and the tree-ring parameter analyzed. The ALT δ¹³C value (−3.3σ) in 1817 CE and YAK MXD (−2.4σ) in 1816 can be seen as an exception to the rule here as they evidenced extreme values (Fig. S1).

Finally, the Pinatubo eruption is mainly captured by the MXD (−2.8σ) and CWT (−2.2σ) chronologies from YAK in 1992 CE. Simultaneous decreases in all tree-ring proxies from ALT are observed in 1993 (Fig. 2), which, however, cannot be classified as extreme (Fig. S1). Overall, the SEA (Fig. 3) shows that volcanic eruptions centered around 535, 540, 1257, 1640, 1815, and 1991 CE have led to decreasing values for all tree-ring proxies in the 2 years afterwards. A short-term response by 2 years after the eruptions is observed in the TRW and CWT proxies for TAY, while for YAK and ALT, the CWT decrease lasts longer (up to 5–6 years in ALT and YAK, respectively) (Fig. 3). The δ¹⁸O isotope chronologies (z-score) show a distinct decrease the year after the eruptions. At ALT, however, the duration of negative anomalies was shorter (5 years) than at the high-latitude TAY (12 years) and YAK (9 years) sites. At the YAK site, two negative years followed the events, intermitted with one positive value, to remain negative during the following 7 years. The duration of negative anomalies recorded in δ¹³C values (z-score) is also longer at the high-latitude YAK site – 10 years after the eruptions and 13 years at TAY compared to 7 years at ALT (Fig. 3).

The largest decrease in MXD values (in terms of z-score) is found at the high-latitude YAK site. The SEA for TRW, MXD, δ¹³C, and CWT from YAK, as well as TRW and MXD from ALT, shows a more drastic decrease in values during the first year when compared to other proxies and study sites (Fig. 3).

Figure 3 Superposed epoch analysis (SEA) of δ¹⁸O, δ¹³C, CWT, TRW, and MXD chronologies for the Yakutia (YAK), Taimyr (TAY), and Altai (ALT) sites, summarizing negative anomalies 15 years before and 20 years after the volcanic eruptions in 535, 540, 1257, 1640, 1815, and 1991 CE. Statistically negative anomalies are marked with a red star (p<0.05).

3.2 Tree-ring proxies versus meteorological series

3.2.1 Monthly air temperatures and sunshine duration

Bootstrapped functions calculated for the instrumental period (1950–2000) show significant positive correlations (p<0.05) between TRW and MXD chronologies and mean summer (June–July) temperatures at all sites. Temperatures at the beginning (June) and the end of the growing season (mid-August) influenced the MXD chronology in ALT (r=0.57) and YAK (r=0.55), respectively (Fig. 4). July temperatures appear to be a key factor for determining tree growth as they significantly impact CWT, δ¹³C, and δ¹⁸O (with the exception of TAY for the latter) chronologies (r=0.28–0.60) at YAK and ALT.

Correlation analysis between July temperature and July sunshine duration indicates significant (p<0.05) correlation for YAK (r=0.56) and ALT (r=0.34). July sunshine duration is strongly and positively correlated with δ¹⁸O in larch tree-ring cellulose chronologies from YAK (r=0.73) and ALT (r=0.51) for the period 1961–2000 (available sunshine duration dataset).

Figure 4 Significant correlation coefficients among the tree-ring parameters TRW, MXD, CWT, δ¹³C, and δ¹⁸O versus weather station data: temperature (T, red), precipitation (P, blue), vapor pressure deficit (VPD, green), and sunshine duration (S, yellow) from September of the previous year to August of the current year calculated for three study sites. Table 2 lists stations and periods used in the analysis.

3.2.2 Monthly precipitation

The strongest July precipitation signal is observed at ALT ($r=-\mathrm{0.54}$) and TAY ($r=-\mathrm{0.51}$) with δ¹³C chronologies (p<0.05). In addition, the ALT data show a significant relationship (p<0.05) between March precipitation and TRW (r=0.37) and MXD (r=0.32), whereas April precipitation correlates positively with CWT (r=0.34). At YAK, July precipitation showed a negative relationship with δ¹⁸O in tree-ring cellulose ($r=-\mathrm{0.34}$; p<0.05) only.

3.2.3 Vapor pressure deficit

June VPD is significantly and positively correlated with the δ¹⁸O chronology from ALT (r=0.67 and p<0.05, respectively) for the period 1950–2000. The δ¹³C in tree-ring cellulose from YAK correlates with July VPD only (r=0.69, p<0.05). We did not find a significant influence of VPD in TAY tree-ring and stable isotope parameters.

3.2.4 Synthesis of the climate data analysis

In summary, during the instrumental period of weather station observations (Table 2) summer temperature impacts TRW, MXD, and CWT at the high-latitude sites (YAK, TAY), while summer precipitation affects stable carbon and oxygen isotopes (YAK, TAY, ALT), sunshine duration (YAK, ALT), and vapor pressure deficit (YAK, ALT).

3.3 Response of Siberian larch trees to climatic changes after the major volcanic eruptions

Based on the statistical analysis above for the calibration period, we assumed that these relationships would not change over time and will provide information about climatic changes during past volcanic periods (Fig. 5).

Figure 5 Responses of larch trees from Yakutia (YAK), Taimyr (TAY), and Altai (ALT) to volcanic eruptions (Table 1). Squares, rhombs, circles, and triangles indicate the years following each eruption that can be considered very extreme (negative values < 5th percentile of the PDFs, intensive color), extreme (negative values > 5th, < 10th percentile of the PDFs, light color), and non-extreme (> 10th percentile of the PDFs, white). July temperature changes are presented with squares. Summer vapor pressure deficit (VPD) variability is shown with circles. July precipitation is presented with rhombs, and July sunshine duration is shown as triangles.

3.3.1 Temperature proxies

We found strong negative summer air temperature anomalies at all sites after the 535 and 1257 CE volcanic eruptions. The temperature decrease was found in the TRW and CWT datasets at all sites and also in the MXD datasets at YAK and ALT (Fig. 5). For the volcanic eruptions in later centuries, the evidence for a decrease in temperature was not as pronounced.

Whereas no strong decline in summer temperature was found at ALT in 1642 CE, we observe a slight decrease in TRW, MXD, and CWT values in 1643. By contrast, a cold summer was recorded by most tree-ring parameters at YAK, except for δ¹⁸O. The absence of strong cooling is even more striking during the years that followed the 1815 CE Tambora eruption. In 1816 CE, only the MXD from YAK shows colder than normal conditions (Fig. 5). 1992 CE was recorded as a cold year in MXD and CWT from YAK, but again not at the other regions or by other proxies.

3.3.2 Moisture proxies: precipitation and VPD

Based on climatological analysis with the local weather station data (Table 2, Fig. 4) for all studied sites we considered δ¹³C in tree-ring cellulose as a proxy for precipitation and vapor pressure deficit changes. Yet, CWT from ALT could be considered a proxy with a mixed temperature and precipitation signal (Fig. 4). Accordingly, the δ¹³C values point to humid summers at YAK in 536, 1258, 1259, 1642, and 1643, at TAY in 536–538 and 1259, and at ALT in the years of 541, 542, 1258, 1259, and 1817. Compared to other proxies and sites, the years 536–538 were neither extremely humid nor dry at ALT (Fig. 5). No negative hydrological anomalies were recorded after the Tambora and Pinatubo eruptions at the high-latitude sites (YAK, TAY). However, positive anomalies were recorded in δ¹³C values, pointing to dry conditions at TAY in 1817 CE (Fig. 2). A rather wet summer was reconstructed for the high-altitude ALT site in 1817 CE compared to 1816 (Fig. 5). Overall, there were mostly humid anomalies after the eruptions at YAK.

3.3.3 Sunshine duration proxies

Instrumental measurements of sunshine duration (Table 2) at YAK and ALT during the recent period showed a significant link with δ¹⁸O cellulose. The sunshine duration decreased after various eruptions at YAK (538, 542, 1258, and 1259) and in 536 at the ALT site.

4 Discussion

In this paper, we analyze climatic anomalies in years following selected large volcanic eruptions using long-term tree-ring multi-proxy chronologies for δ¹³C and δ¹⁸O, TRW, MXD, and CWT for high-latitude (YAK, TAY) and high-altitude (ALT) sites. Since trees as living organisms respond to various climatic impacts, the carbon assimilation and growth patterns accordingly leave unique “fingerprints” in the photosynthates, which is recorded in the wood in the tree rings specifically and individually for each proxy.

4.1 Evaluation of the applied proxies in Siberian tree-ring data

This study clearly shows that each proxy has to be analyzed and interpreted specifically for its validity at each studied site and evaluated for its suitability for the reconstruction of abrupt climatic changes.

The TRW in temperature-limited environments is an indirect proxy for summer temperature reconstructions, as growth is a temperature-controlled process. Temperature clearly determines the duration of the growing season and the rate of cell division (Cuny et al., 2014). Accordingly, low temperature in the growing season is recorded by narrow tree rings. The upper limit of temperature is specific to tree species and biomes. In most cases, tree growth is limited by drought rather than by high temperatures, since water shortage and VPD increase with increasing temperature. Still, this does not make TRW a suitable proxy to determine the influence of water availability and air humidity, especially at temperature-limited sites.

MXD chronologies obtained for the Eurasian subarctic record mainly a July–August temperature signal (Vaganov et al., 1999; Sidorova et al., 2010; Büntgen et al., 2016) and add valuable information about climate conditions toward the end of the growth season. Similarly, CWT is an anatomical parameter, which contains information on the carbon sink limitation of the cambium due to extreme cold conditions (Panyushkina et al., 2003; Fonti et al., 2013; Bryukhanova et al., 2015). There is a strong signal of low cell number within a growing season; for example, strongly decreasing CWT in 537 CE at YAK and the formation of frost rings at ALT (in 536–538 and 1259 CE) have been shown in our study.

Low δ¹³C values can be explained by a reduction in photosynthesis caused by volcanic dust veils. For the distinction whether δ¹³C is predominantly determined by A_N or g_l the combined evaluation with δ¹⁸O or TRW is needed. High δ¹⁸O values indicate high VPD, which induces a reduction in stomatal conductance, reducing the back diffusion of depleted water molecules from the ambient air. This confirms a sunny 1993 CE at ALT with mild weather conditions according to observational data from the closest weather station (Table 2). Interestingly, we also find less negative values for δ¹³C in the same period. This shows that the two isotopes correlate with each other and indicates the need for a combined evaluation of C and O isotopes (Scheidegger et al., 2000) taking into account precautions as suggested by Roden and Siegwolf (2012).

4.2 Lag between volcanic events and response in tree rings

Most of the events discussed suggest a lag between the eruption and the tree-ring response for 1 year or more (Fig. 3). This lag is explained by the tree's use of stored carbohydrates, which are the substrate for needle and early wood production. These stored carbohydrates carry the isotopic signal of previous years and depend on their remobilization; as such the signals may be masked in freshly produced biomass. The delayed signal could also reflect the time needed for the dust veil to be transported to the study regions.

4.3 Temperature and sunshine duration changes after stratospheric volcanic eruptions

Correlation functions show that MXD and CWT (with the exception of TAY in the latter case), and to a lesser extent TRW chronologies, portray the strongest signals for summer (June–August) temperatures. In addition, significant information about sunshine duration can be derived from the YAK and ALT δ¹⁸O series. Thus, we hypothesize that extremely narrow TRW and very negative anomalies observed in the MXD and CWT chronologies of YAK, and to a lesser extent at ALT, along with low δ¹⁸O values, reflect cold and low sunshine duration conditions in summer. Presumably, the temperatures were below the threshold values for growth over much of the growing season (Körner, 2015). This hypothesis of a generalized regional cooling after both eruptions is further confirmed by the occurrence of frost rings at the ALT site in 538 and 1259 CE (Myglan et al., 2008; Guillet et al., 2017), as well as in neighboring Mongolia (D'Arrigo et al., 2001). The unusual cooling in 536–542 CE is also evidenced by a very small number of cells formed at YAK (Churakova (Sidorova) et al., 2014). Although δ¹⁸O is an indirect proxy for needle temperature, low δ¹⁸O values in 538, 542, 1258, and 1259 CE for YAK and in 536 CE for ALT are a result of low irradiation, leading to low temperature and low VPD (high stomatal conductance), both likely a result of volcanic dust veils.

Similarly, in the aftermath of the Samalas eruption, the persistence of summer cooling is limited to 1258 and 1259 CE at the three studied sites, which is in line with the findings of Guillet et al. (2017). Interestingly, a slight decrease in oxygen isotope chronologies, which can be related to low levels of summer sunshine duration (i.e., low leaf temperatures), allows for the hypothesis that cool conditions could have prevailed.

For all later high-magnitude CE eruptions, temperature-sensitive tree-ring proxies do not evidence a generalized decrease in summer temperatures. Paradoxically, the impacts of the Tambora eruption, known for its triggering of a widespread “year without summer” (Harrington, 1992), only induced abnormal MXD at YAK in 1816, but no anomalies are observed at TAY and ALT, except for the positive deviation of δ¹³C at TAY and the negative anomaly at ALT in 1817 CE (Figs. 2, 5, and S1). While these findings may seem surprising, they are in line with the TRW and MXD reconstructions of Briffa et al. (1998) and Guillet et al. (2017), who found limited impacts of the 1815 CE Tambora eruption in eastern Siberia and Alaska using TRW and MXD data only. The inclusion of CWT chronologies by Barinov et al. (2018) confirms the absence of a significant cooling signal after the second largest eruption of the last millennium (1815 CE) in larch trees of the Altai–Sayan mountain region.

Finally, in 1992 CE, our results evidence cold conditions at YAK, which is consistent with weather observations showing that the below-average anomalies in summer temperatures (after Pinatubo eruption) were indeed limited to northeastern Siberia (Robock, 2000). As both isotopes indicate a reduction in stomatal conductance, we found that warm (in agreement with MXD and CWT) and dry conditions were prevalent at ALT at this time. This isotopic constellation was confirmed by the positive relationships among VPD, δ¹⁸O, and δ¹³C at ALT.

However, temperature and sunshine duration are not always highly coherent over time due to the influence of other factors, like Arctic oscillations as suggested for Fennoscandia regions by Loader et al. (2013).

4.4 Moisture changes

Water availability is a key parameter for Siberian trees as they are growing under extremely continental conditions with hot summers and cold winters, with very low annual precipitation (Table 2). Permafrost plays a crucial role and can be considered a buffer for additional water sources during hot summers (Sugimoto et al., 2002; Boike et al., 2013; Saurer et al., 2016). Yet, thawed permafrost water is not always available to roots due to the surficial structure of the root plate or an extremely cold water temperature (close to 0 ^∘C), which can hardly be utilized by trees (Churakova (Sidorova) et al., 2016). Thus, Siberian trees are highly susceptible to drought induced by dry and warm air during July, and therefore the stable carbon isotopes can be sensitive indicators of such conditions. After volcanic eruptions, however, low light intensity due to dust veils induces low temperatures and reduced VPD, the driver for evapotranspiration. Under such conditions drought stress is unlikely to occur. However, the transition phases with changes from cool and moist to warm and dry conditions are more critical when drought is more likely to occur.

In our study, higher δ¹³C values in tree-ring cellulose indicate increasing drought conditions as a consequence of reduced precipitation for 2 years after the 1815 CE volcanic eruption at the TAY site. No further extreme hydroclimatic anomalies occurred at Siberian sites in the aftermath of the Pinatubo eruption.

4.5 Synthetized interpretation from the multi-parameter tree-ring proxies

Our analysis demonstrates the added value of a tree-ring-derived multi-proxy approach to better capture climatic variability after large volcanic eruptions. Besides the well-documented effects of temperature derived from TRW, MXD, and CWT, stable carbon and oxygen isotopes in tree-ring cellulose provide important and complementary information about moisture and sunshine duration changes (an indirect proxy for leaf temperature effective air-to-leaf VPD) after stratospheric volcanic eruptions.

Our results reveal the complex behavior of the Siberian climatic system to the stratospheric volcanic eruptions of the Common Era. The 535 and 1257 CE Samalas eruptions caused substantial cooling – very likely induced by dust veils (Churakova (Sidorova) et al., 2014; Guillet et al., 2017; Helama et al., 2018) – as well as humid conditions at both high-latitude and high-altitude sites. Conversely, only local and limited climate responses were observed after the 1641 CE Parker, 1815 Tambora, and 1991 Pinatubo eruptions. Similar site-dependent impacts referred to the coldest summers of the last millennium in the Northern Hemisphere based on TRW and MXD reconstructions (Schneider et al., 2015; Stoffel et al., 2015; Wilson et al., 2016; Guillet et al., 2017). This absence of widespread and intense cooling or missing drastic changes in hydrological regime over vast regions of Siberia may result from the location and strength of the volcanic eruption, atmospheric transmissivity, and/or the modulation of radiative forcing effects by regional climate variability. These results are consistent with other regional studies that interpreted the spatial–temporal heterogeneity of tree responses to past volcanic events (Wiles et al., 2014; Esper et al., 2017; Barinov et al., 2018) in terms of regional climates.

5 Conclusions

In this study, we demonstrate that the consequences of large volcanic eruptions on climate are rather complex between sites and among events. The different locations and magnitudes of eruptions, but also regional climate variability, may explain some of this heterogeneity. We show that each tree-ring and isotope proxy alone cannot provide the full information for the volcanic impact on climate, but that they, when combined, contribute to the formation of the full picture, which is critical for a comprehensive description of climate dynamics induced by volcanism and the inclusion of these phenomena in global climate models.

Analyses with a larger number of samples in investigations of Siberian and other Northern Hemisphere sites will indeed provide higher certainty in terms of data interpretation of climatic dynamics in these boreal regions. However, the multi-proxy approach as applied in our study also provides a strong set of complementary information to the research field, as it allows for the refinement of the interpretations and thus improves our understanding of the heterogeneity of climatic signals after CE stratospheric volcanic eruptions, as recorded in multiple tree-ring and stable isotope parameters.