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 δ13C and δ18O 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 δ13C, 536 was extremely humid at YAK, as was 537–538 in TAY. No extreme hydroclimatic [...] CHURAKOVA, Olga, et al. Siberian tree-ring and stable isotope proxies as indicators of temperature and moisture changes after major stratospheric volcanic eruptions. Climate of the Past, 2019, vol. 15, no. 2, p. 685-700 DOI : 10.5194/cp-15-685-2019


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., 2017Esper et al., , 2018Guillet 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 (δ 13 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 (δ 13 C) and oxygen (δ 18 O) isotopes in tree rings has been employed only rarely to trace volcanic eruptions in high-latitude or highaltitude proxy records (Churakova (Sidorova) et al., 2014).
The application of TRW, MXD, and cell wall thickness (CWT), as well as δ 13 C and δ 18 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 Sidorova et al., 2008Sidorova et al., , 2010Sidorova et al., , 2011Churakova (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, δ 13 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 δ 13 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.
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 NA -not available. -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, δ 13 C, and δ 18 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).

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  Hughes et al., 1999;Sidorova and Naurzbaev, 2002;Sidorova, 2003, for YAK;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).
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  Sidorova and Naurzbaev, 2002;Sidorova, 2003), eastern Taimyr (TAY -green, c) , and Russian Altai (ALTred, d) (Myglan et al., 2009). Photos show the larch stands at the YAK, TAY (Mukhtar M. Naurzbaev), and ALT (Vladimir S. Myglan) sites. 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 treering 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).

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.

Image analysis of cell wall thickness (CWT)
Analysis of wood anatomy was performed for all studied periods with an AxioVision scanner (Carl Zeiss, Germany). 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 −1 ), and vapor pressure deficit (VPD, kPa) data were downloaded from the meteorological database: http://aisori.meteo.ru/ClimateR (last access: January 2018).

Length of
Thawing Meteorological parameters growing permafrost Annual air Annual 1950-2000 1966-2000 1961-2000 1950-2000 50-70 a 20-50 a −14.7 205 1950-2000 1966-2000 1961-2000 1950-  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).

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 −3 (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 (δ 13 C) and oxygen (δ 18 O) isotopes in tree-ring cellulose During photosynthetic CO 2 assimilation 13 CO 2 is discriminated against 12 CO 2 , leaving the newly produced assimilates depleted of 13 C. The carbon isotope discrimination ( 13 ) is partitioned in the diffusional component with a = 4.4 ‰ and the biochemical fractionation with b = 27 ‰ for C 3 plants during carboxylation via Rubisco. The 13 is directly proportional to the c i /c a ratio, where c i is the leaf intercellular, and c a the ambient CO 2 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 13 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 2 and H 2 O gas exchange is strongly interlinked with Cisotope fractionation, 13 is controlled by the same environmental variables, i.e., PAR, CO 2 , VPD, and temperature (Farquhar et al., 1982(Farquhar et al., , 1989Cernusak et al., 2013). The oxygen isotopic compositions of tree-ring cellulose record the δ 18 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 690 O. V. Churakova (Sidorova) et al.: Siberian tree-ring and stable isotope proxies 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 δ 18 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 13 C/ 12 C and 18 O/ 16 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 13 C/ 12 C ratio was determined separately by combustion under oxygen excess at a reactor temperature of 1020 • C. Samples for 18 O/ 16 O ratio measurements were pyrolyzed to CO at 1080 • C . 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): where R sample is the molar fraction of the 13 C/ 12 C or 18 O/ 16 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 (VS-MOW) for oxygen. The precision is σ ± 0.1 ‰ for carbon and σ ± 0.2 ‰ for oxygen. To remove the atmospheric δ 13 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 δ 13 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 δ 18 O cellulose series were not detrended.

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).

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 nonvolcanic years, we computed probability density functions (PDFs; Stirzaker, 2003) for each study site and for each treering 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.

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 δ 18 O chronology recorded a drastic decrease in 536 CE with −4.8σ (Figs. 2 and S1). A δ 18 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 δ 18 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, δ 18 O, while it is less extreme in CWT and δ 13 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 δ 13 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 δ 18 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 δ 13 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, δ 13 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).

Monthly air temperatures and sunshine duration
Bootstrapped functions calculated for the instrumental period  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 692 O. V. Churakova (Sidorova) et al.: Siberian tree-ring and stable isotope proxies appear to be a key factor for determining tree growth as they significantly impact CWT, δ 13 C, and δ 18 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 δ 18 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).

Vapor pressure deficit
June VPD is significantly and positively correlated with the δ 18 O chronology from ALT (r =0.67 and p < 0.05, respectively) for the period 1950-2000. The δ 13 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.

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).

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).

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 δ 18 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.

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 δ 13 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 δ 13 C values point to humid summers at YAK in 536, 1258536, , 1259536, , 1642536, , and 1643536, , 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 δ 13 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.

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

Discussion
In this paper, we analyze climatic anomalies in years following selected large volcanic eruptions using long-term treering multi-proxy chronologies for δ 13 C and δ 18 O, TRW, MXD, and CWT for high-latitude (YAK, TAY) and highaltitude (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. . Significant correlation coefficients among the tree-ring parameters TRW, MXD, CWT, δ 13 C, and δ 18 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.

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 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 lim-itation 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 δ 13 C values can be explained by a reduction in photosynthesis caused by volcanic dust veils. For the distinction whether δ 13 C is predominantly determined by A N or g l the combined evaluation with δ 18 O or TRW is needed. High δ 18 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 δ 13 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).  (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 nonextreme (> 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.

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.

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 δ 18 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 δ 18 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)  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, temperaturesensitive 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 δ 13 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 696 O. V. Churakova (Sidorova) et al.: Siberian tree-ring and stable isotope proxies by the positive relationships among VPD, δ 18 O, and δ 13 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).

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;). 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 δ 13 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.

Synthetized interpretation from the multi-parameter tree-ring proxies
Our analysis demonstrates the added value of a tree-ringderived multi-proxy approach to better capture climatic variability after large volcanic eruptions. Besides the welldocumented effects of temperature derived from TRW, MXD, and CWT, stable carbon and oxygen isotopes in treering 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 lo-cal 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.

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 treering 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 multiproxy 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. Data availability. Research data are available upon request from the corresponding author Olga V. Churakova (Sidorova) (olga.churakova@hotmail.com) and will be provided via PAGES ISO2K.
Author contributions. TRW analysis was performed at the V.N. Sukachev Institute of Forest SB RAS by OVCS, DVO, VSM, and OVN. CWT analysis was carried out at the V.N. Sukachev Institute of Forest SB RAS, Krasnoyarsk, Russia, by MVF and at the University of Arizona by IPP. Stable isotope analysis was conducted at the Paul Scherrer Institute (PSI) by OVCS, MS, and RS. MXD measurements were realized with a DENDRO Walesch 2003 densitometer at WSL and at the V.N. Sukachev Institute of Forest SB RAS, Krasnoyarsk, Russia, by OVCS and AVK. Samples from YAK and TAY were collected by MMN. All authors contributed significantly to the data analysis and paper writing.