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

Abstract

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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ and${\mathit{\delta }}^{\mathrm{18}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{18}}\mathrm{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 ${}^{\circ }$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 ${}^{\circ }$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 ${}^{\circ }$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 (${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 (${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$) and oxygen (${\mathit{\delta }}^{\mathrm{18}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{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, ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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_{\text{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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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.

			Volcanic
Study period	Date of eruption	Volcano	explosivity
(CE)	(Day/month/year)	name	index (VEI)	Location, coordinates	References
520–560	NA/NA/535	Unknown	?	Unknown	Stothers (1984)
	NA/NA/540	Unknown	?	Unknown	Sigl et al. (2015);
					Toohey and Sigl (2017)
1242–1286	NA/05–10/1257	Samalas	7	Indonesia, 8.42${}^{\circ }$ N, 116.47${}^{\circ }$ E	Stothers (2000);
					Lavigne et al. (2013);
					Sigl et al. (2015)
1625–1660	26/12/1640	Parker	5	Philippines, 6${}^{\circ }$ N, 124${}^{\circ }$ E	Zielinski et al. (1994, 2000)
1790–1835	10/04/1815	Tambora	7	Indonesia, 8${}^{\circ }$ S, 118${}^{\circ }$ E	Zielinski et al. (1994, 2000)
1950–2000	15/06/1991	Pinatubo	6	Philippines, 15${}^{\circ }$ N, 120${}^{\circ }$ E	Zielinski et al. (1994),
					Sigl et al. (2015)

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, ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$, and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{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${}^{\circ }$ N, 148${}^{\circ }$ E) and eastern Taimyr (TAY: 70${}^{\circ }$ N, 103${}^{\circ }$ E) and mountain permafrost in Altai (ALT: 50${}^{\circ }$ N, 89${}^{\circ }$ 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>\mathrm{0.91}$; $p<\mathrm{0.05}$) with gridded data (http://climexp.knmi.nl, last access: January 2018). However, poor correlation is found with precipitation data ($r<\mathrm{0.45}$; $p<\mathrm{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${}^{-\mathrm{1}}$) at all sites. The growing season calculated with the tree growth threshold of $+$5 ${}^{\circ }$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${}^{-\mathrm{1}}$ in summer) compared to ALT (ca. 18 h day${}^{-\mathrm{1}}$ 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.

[Figure omitted. See PDF]

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$, ${}^{\circ }$C), precipitation ($P$, mm), sunshine duration ($S$, h month${}^{-\mathrm{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
Site	Tree	Location	Weather	$T$ (${}^{\circ }$C)	$P$ (mm)	$S$ (h	VPD	season	depth	temperature	precipitation
	species		station	(periods)		month${}^{-\mathrm{1}}$)	(kPa)	(day)	(max, cm)	(${}^{\circ }$C)	(mm)
YAK	Larix cajanderi Mayr.	69${}^{\circ }$ N, 148${}^{\circ }$ E	Chokurdach 62${}^{\circ }$ N, 147${}^{\circ }$ E, 61 m a.s.l.	1950–2000	1966–2000	1961–2000	1950–2000	50–70${}^{\text{a}}$	20–50${}^{\text{a}}$	$-\mathrm{14.7}$	205
TAY	Larix gmelinii Rupr.	70${}^{\circ }$ N, 103${}^{\circ }$ E	Khatanga 71${}^{\circ }$ N, 102${}^{\circ }$ E, 33 m a.s.l.	1950–2000	1966–2000	1961–2000	1950–2000	90${}^{\text{b}}$	40–60${}^{\text{b}}$	$-\mathrm{13.2}$	269
ALT	Larix sibirica Ledeb.	50${}^{\circ }$ N, 89${}^{\circ }$ E	Mugur Aksy 50${}^{\circ }$ N, 90${}^{\circ }$ E 1850 m a.s.l.	1963–2000	1966–2000			90–120${}^{\text{c}}$	80–100${}^{\text{c}}$	$-\mathrm{2.7}$	153
			Kosh-Agach 50${}^{\circ }$ N, 88${}^{\circ }$ E 1758 m a.s.l.			1961–2000	1950–2000

${}^{\text{a}}$ Abaimov et al. (1996); Hughes et al. (1999); Churakova (Sidorova) et al. (2016); ${}^{\text{b}}$ Naurzbaev et al. (2002); ${}^{\text{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, ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$, and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{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=\mathrm{4}$) and stable isotopes ($n=\mathrm{4}$) than for TRW ($n=\mathrm{12}$) or MXD ($n=\mathrm{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>\mathrm{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${}^{-\mathrm{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 (${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$) and oxygen (${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$) isotopes in tree-ring cellulose

During photosynthetic ${\mathrm{CO}}_{\mathrm{2}}$ assimilation ${}^{\mathrm{13}}{\mathrm{CO}}_{\mathrm{2}}$ is discriminated against ${}^{\mathrm{12}}{\mathrm{CO}}_{\mathrm{2}}$, leaving the newly produced assimilates depleted of ${}^{\mathrm{13}}\mathrm{C}$. The carbon isotope discrimination (${}^{\mathrm{13}}\mathrm{\Delta }$) is partitioned in the diffusional component with $a=\mathrm{4.4}$ ‰ and the biochemical fractionation with $b=\mathrm{27}$ ‰ for C${}_{\mathrm{3}}$ plants during carboxylation via Rubisco. The ${}^{\mathrm{13}}\mathrm{\Delta }$ is directly proportional to the $c_{\text{i}}/c_{\text{a}}$ ratio, where $c_{\text{i}}$ is the leaf intercellular, and $c_{\text{a}}$ the ambient ${\mathrm{CO}}_{\mathrm{2}}$ concentration. This ratio reflects the balance between stomatal conductance ($g_{\text{l}}$) and photosynthetic rate ($A_{\text{N}}$). A decrease in $g_{\text{l}}$ at a given $A_{\text{N}}$ results in a decrease in ${}^{\mathrm{13}}\mathrm{\Delta }$ as $c_{\text{i}}/c_{\text{a}}$ decreases and vice versa. The same is true when $A_{\text{N}}$ increases or decreases at a given $g_{\text{l}}$. Since ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ gas exchange is strongly interlinked with C-isotope fractionation, ${}^{\mathrm{13}}\mathrm{\Delta }$ is controlled by the same environmental variables, i.e., PAR, ${\mathrm{CO}}_{\mathrm{2}}$, VPD, and temperature (Farquhar et al., 1982, 1989; Cernusak et al., 2013). The oxygen isotopic compositions of tree-ring cellulose record the ${\mathit{\delta }}^{\mathrm{18}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{18}}\mathrm{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. $\mathit{\alpha }$-cellulose extraction was performed according to the method described by Boettger et al. (2007). For analyses of the ${}^{\mathrm{13}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ and ${}^{\mathrm{18}}\mathrm{O}{/}^{\mathrm{16}}\mathrm{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 ${}^{\mathrm{13}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ ratio was determined separately by combustion under oxygen excess at a reactor temperature of 1020 ${}^{\circ }$C. Samples for ${}^{\mathrm{18}}\mathrm{O}{/}^{\mathrm{16}}\mathrm{O}$ ratio measurements were pyrolyzed to CO at 1080 ${}^{\circ }$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): 1$$\mathit{\delta }\text{sample}=R_{\text{sample}}/R_{\text{standard}}-\mathrm{1},$$ where $R_{\text{sample}}$ is the molar fraction of the ${}^{\mathrm{13}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ or ${}^{\mathrm{18}}\mathrm{O}{/}^{\mathrm{16}}\mathrm{O}$ ratio of the sample and $R_{\text{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 $\mathit{\sigma }\pm \mathrm{0.1}$ ‰ for carbon and $\mathit{\sigma }\pm \mathrm{0.2}$ ‰ for oxygen. To remove the atmospheric ${\mathit{\delta }}^{\mathrm{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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{18}}\mathrm{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), ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ (red), and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{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.

[Figure omitted. See PDF]

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 ($-\mathrm{1.8}\mathit{\sigma }$) and 541 ($-\mathrm{2.4}\mathit{\sigma }$) 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 ($-\mathrm{3.1}\mathit{\sigma }$).

Furthermore, we found decreasing MXD values at ALT ($-\mathrm{4.4}\mathit{\sigma }$) in 537 CE and YAK ($-\mathrm{2.8}\mathit{\sigma }$) 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 ($-\mathrm{1.9}\mathit{\sigma }$) and 541 ($-\mathrm{2.4}\mathit{\sigma }$), whereas the response was smaller in TAY and ALT for the same years (Fig. 2). The ALT ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ chronology recorded a drastic decrease in 536 CE with $-\mathrm{4.8}\mathit{\sigma }$ (Figs. 2 and S1). A ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ decrease for YAK was found after the 1257 CE Samalas eruption in 1258 CE ($-\mathrm{1.5}\mathit{\sigma }$) and in 1259 ($-\mathrm{2.9}\mathit{\sigma }$), which is opposite to the increased ${\mathit{\delta }}^{\mathrm{18}}\mathrm{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 ($-\mathrm{2.3}\mathit{\sigma }$). 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, ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$, while it is less extreme in CWT and ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ value ($-\mathrm{3.3}\mathit{\sigma }$) in 1817 CE and YAK MXD ($-\mathrm{2.4}\mathit{\sigma }$) 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 ($-\mathrm{2.8}\mathit{\sigma }$) and CWT ($-\mathrm{2.2}\mathit{\sigma }$) 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 ${\mathit{\delta }}^{\mathrm{18}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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, ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$, ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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<\mathrm{0.05}$).

[Figure omitted. See PDF]

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<\mathrm{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=\mathrm{0.57}$) and YAK ($r=\mathrm{0.55}$), respectively (Fig. 4). July temperatures appear to be a key factor for determining tree growth as they significantly impact CWT, ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$, and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ (with the exception of TAY for the latter) chronologies ($r=\mathrm{0.28}$–0.60) at YAK and ALT.

Correlation analysis between July temperature and July sunshine duration indicates significant ($p<\mathrm{0.05}$) correlation for YAK ($r=\mathrm{0.56}$) and ALT ($r=\mathrm{0.34}$). July sunshine duration is strongly and positively correlated with ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ in larch tree-ring cellulose chronologies from YAK ($r=\mathrm{0.73}$) and ALT ($r=\mathrm{0.51}$) for the period 1961–2000 (available sunshine duration dataset).

Figure 4

Significant correlation coefficients among the tree-ring parameters TRW, MXD, CWT, ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$, and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{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.

[Figure omitted. See PDF]

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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ chronologies ($p<\mathrm{0.05}$). In addition, the ALT data show a significant relationship ($p<\mathrm{0.05}$) between March precipitation and TRW ($r=\mathrm{0.37}$) and MXD ($r=\mathrm{0.32}$), whereas April precipitation correlates positively with CWT ($r=\mathrm{0.34}$). At YAK, July precipitation showed a negative relationship with ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ in tree-ring cellulose ($r=-\mathrm{0.34}$; $p<\mathrm{0.05}$) only.

3.2.3 Vapor pressure deficit

June VPD is significantly and positively correlated with the ${\mathit{\delta }}^{\mathrm{18}}$O chronology from ALT ($r=$0.67 and $p<\mathrm{0.05}$, respectively) for the period 1950–2000. The ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ in tree-ring cellulose from YAK correlates with July VPD only ($r=\mathrm{0.69}$, $p<\mathrm{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.

[Figure omitted. See PDF]

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 ${\mathit{\delta }}^{\mathrm{18}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{18}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ values can be explained by a reduction in photosynthesis caused by volcanic dust veils. For the distinction whether ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ is predominantly determined by $A_{\text{N}}$ or $g_{\text{l}}$ the combined evaluation with ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ or TRW is needed. High ${\mathit{\delta }}^{\mathrm{18}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{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).

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 ${\mathit{\delta }}^{\mathrm{18}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{18}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ is an indirect proxy for needle temperature, low ${\mathit{\delta }}^{\mathrm{18}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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, ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$, and ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${}^{\circ }$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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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.

Data availability

Research data are available upon request from the corresponding author Olga V. Churakova (Sidorova) ([email protected]) and will be provided via PAGES ISO2K.

The supplement related to this article is available online at: https://doi.org/10.5194/cp-15-685-2019-supplement.

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.

Competing interests

The authors declare that they have no conflict of interest.

Acknowledgements

This work was supported by a Marie Curie International Incoming Fellowship (EU_ISOTREC 235122), a Reintegration Marie Curie Fellowship (909122), and grants to the following: a UFZ scholarship (2006), RFBR (09-05-98015_r_sibir_a), granted to Olga V. Churakova (Sidorova); SNSF to Matthias Saurer (200021_121838/1); an Era.Net RusPlus project granted to Markus Stoffel (SNF IZRPZ0_164735); RFBR (no. 16-55-76012 Era_a) granted to Eugene A. Vaganov; and a project grant to Vladimir S. Myglan RNF, Russian Scientific Fund (no. 15-14-30011). Alexander V. Kirdyanov was supported by the Ministry of Education and Science of the Russian Federation (no. 5.3508.2017/4.6) and RSF (no. 14-14-00295). We acknowledge a Scientific School (3297.2014.4) grant to Eugene A. Vaganov, US National Science Foundation (NSF) grants (no. 9413327, no. 970966, no. 0308525) to Malcolm K. Hughes, and US CRDF grant no. RC1-279 to Malcolm K. Hughes and Eugene A. Vaganov. We thank Tatjana Boettger for her support and access to the stable isotope facilities within the framework of the UFZ Haale/Saale scholarship 2006 and stable isotope facilities at the Paul Scherrer Institute (PSI), Switzerland; we thank Anne Verstege and Daniel Nievergelt for their help with sample preparation for the MXD and Paolo Cherubini for providing lab access at the Swiss Federal Institute for Forest, Snow and Landscape Research (WSL).

We thank two anonymous reviewers and the handling editor Juerg Luterbacher for their constructive comments on this paper.

Review statement

This paper was edited by Jürg Luterbacher and reviewed by two anonymous referees.