1 Introduction

In North America, the network of proxy records used for high-resolution climate reconstructions over the last two millennia is dominated by tree-ring chronologies. However, these chronologies are much more abundant in the west than in the east and, importantly, there are only five chronologies spanning more than a millennium north of the 50th parallel (Pages 2k Consortium 2013; http://pastglobalchanges.org/ini/wg/2k-network/intro). In northeastern North America, developing tree-ring chronologies is highly challenging due to short tree longevity, the high frequency and severity of wildfires and the remoteness of many areas (Arseneault et al. 2013). To improve the quality of reconstructions in such regions with scarce proxy data, we can both increase series replication for a single proxy by creating new chronologies, and combine proxies with an independent response to climate forcing. In this paper, we adopted the solution of combining proxies using datasets from the northern Quebec taiga. Six highly replicated millennium-long ring width chronologies were recently developed in this region by Gennaretti et al. (2014) with black spruce [Picea mariana (Mill.) B.S.P] subfossil trees preserved in six lakes. Samples from one of these sites were also used to obtain a millennium-long record of oxygen isotopic ratios (δ18O) in tree-ring cellulose (Naulier et al. 2015). Until now, these ring width and δ18O data have only been used separately in independent summer temperature reconstructions. In principle, merging these data in a multiproxy approach should provide more robust estimates of the past regional climate variations. Indeed, within a multiproxy framework, we can exploit the distinct climatic responses embedded in each proxy to reduce the impact of individual proxy errors.

Multiproxy approaches at the regional scale have already been used with successful results (Sidorova et al. 2012, 2013). For example, McCarroll et al. (2013) improved their temperature reconstruction by combining chronologies of ring width, ring density and annual tree height growth from Scandinavian sites. Boucher et al. (2011) reconstructed the full spectra of a drought index in southern South America using only the reliable periodicities specific to each of their proxies. Tolwinski-Ward et al. (2015) improved their local temperature-soil moisture reconstructions using both ring width and isotopic data within a hierarchical Bayesian approach, allowing a better understanding of the temporal changes in the climatic controls on the proxies. Such Bayesian methods are especially useful for leveraging the multi-proxy information and for assessing uncertainties with a probabilistic perspective. Bayesian models have thus been used to (1) improve climate field reconstructions from multi-proxy networks (Tingley and Huybers 2010), (2) better infer climate variability from nonlinear proxies (Emile-Geay and Tingley 2016), (3) define spatially varying proxy-climate relationships (Tierney and Tingley 2014) or (4) investigate the mechanistic climate controls on the proxies (Tolwinski-Ward et al. 2013).

In this study, we present a new millennium-long δ13C chronology in tree-ring cellulose, which enhances the temperature-sensitive proxy dataset from the northern Quebec taiga. Thus, we use an ensemble of three proxies (ring width, δ18O and δ13C) to provide the first multiproxy regional high-resolution summer temperature reconstruction in northeastern North America over the last millennium (997–2006 CE). A linear Bayesian approach is used to generate sharp and reliable confidence intervals based on the likelihood and the convergence of the proxy models. Finally, the sensitivity of individual proxies to temperature perturbations is evaluated and discussed, focusing the analysis on the response to strong volcanic eruptions because they produce abrupt perturbations after key-dates.

2 Materials and methods

2.1 Proxy and climate data

The tree-ring data from the northern Quebec taiga are composed of six millennium-long ring width chronologies developed with series from 1782 subfossil stem segments and 150 living black spruces (Gennaretti et al. 2014). Subfossil logs were sampled from the water and sediments of the littoral zone of six boreal lakes (the coordinates of the central point are 54.23 N and 71.39 W), while living trees were selected in the lakeshore forest of the same sites. Their cross-sections are stored at the University of Quebec in Rimouski, and the series are already in the public domain (http://www.ncdc.noaa.gov/paleo). Here, we used the median of the 6 millennium-long site-specific chronologies keeping for each site only periods with sample depth greater than five. Individual chronologies were built with the regional curve standardization (RCS) pivot correction method to reduce the impact of varying sampling heights (Autin et al. 2015). Low and high frequencies of the median chronology were treated separately. Low frequencies (LFs) were obtained with a 9-year triangular filter to produce a chronology comparable to that of the stable isotopes (see below). High frequencies (HFs) were obtained by subtracting the LFs from the raw chronology. The bandwidth of the LF filter at the 50% threshold was 0.07 cycles/year. The LF chronology was also transformed to obtain a quasi-perfect Gaussian distribution with the inverse transform sampling technique. This technique is based on a quantile-based transformation and, in some cases, should improve the linearity of the relationship between proxies and normally distributed climate variables (Emile-Geay and Tingley 2016; van Albada and Robinson 2007):

$${y'_i} = \sqrt 2 er{f^{ - 1}}(2P({y_i}) - 1),$$
(1)

where erf represents the Gauss error function, \(P\left( y \right)\) the proxy cumulative distribution function and \({y_i}\) and \(y{'_i}\) the proxy untransformed and transformed values, respectively, for year i. Fig. S1 (see electronic supplementary material) shows the effect (quite low in this case) of this transformation on the LF chronology.

From one of the six sites from the northern Quebec taiga, 60 subfossil logs and 5 living trees were further analyzed at the Delta-lab of the Geological Survey of Canada to extract two millennium-long chronologies of stable isotope ratios in tree-ring cellulose. The δ18O chronology is already in the public domain (Naulier et al. 2015), whereas the δ13C chronology is detailed here and is accessible in the supplementary material (Dataset S1). These chronologies were built with the “offset-pool plus join-point” method (Gagen et al. 2012) to obtain an annual resolution from five tree replicates (this replication was proven to be adequate to obtain robust site chronologies; Naulier et al. 2014) and successive tree cohorts. Assuming that the five time series of each cohort are realizations of the same stochastic process, this method produces chronologies equivalent to time series smoothed by the aforementioned 9-year triangular filter. Indeed, within every cohort of five trees, the rings of each tree were divided into 5-year blocks for the isotopic measurements with an offset of 1 year among trees (Naulier et al. 2015). The δ13C values of the modern part of the chronology were also mathematically corrected for atmospheric δ13C CO2 changes due to fossil fuel combustion (Suess effect), and for plant response to increasing atmospheric CO2 concentrations (McCarroll et al. 2009; McCarroll and Loader 2004; Naulier et al. 2014). The final isotope ratio chronologies are compared with the ring width LF chronology in Fig. S2. As with the ring width chronology, the isotope chronologies were also transformed with the inverse transform sampling technique (Fig. S1).

The monthly climate data for our study area (1901–2010) were downloaded from the Climatic Research Unit (CRU) TS 3.23 climate dataset (Harris et al. 2014). A correlation analysis showed that the mean of the July and August temperature was the common best fit climate signal registered by our proxies (Fig. S3). Thus, this variable was retained for the climate reconstruction. This is consistent with previous findings showing that summer temperatures control the growth and isotopic values of these trees (Gennaretti et al. 2014; Naulier et al. 2014, 2015). Low and high temperature frequencies were treated separately and obtained with the same method as for the ring width chronology such that the LF chronology was comparable with that of the stable isotopes. The relationships between the temperature variable and the proxies over the last century are shown in Fig. S4.

2.2 Bayesian proxy analysis

Our objective was to infer the values of the mean July–August temperature for each year over a past period where the temperature is unknown. This inference was based on a set of known proxy values (D) and temperature observations (T) that overlap over a calibration period (cal; 1905–2006). Using Bayes’ theorem, the posterior distribution of temperature t at each year i can be obtained by:

$$p\left( {{t_i}|{D_i},~{D_{cal}},~{T_{cal}}} \right) = \frac{{p\left( {{D_i}|{t_i},~{D_{cal}},~{T_{cal}}} \right)p\left( {{t_i}|{T_{cal}},~{D_{cal}}} \right)}}{{p\left( {{D_i}|{T_{cal}},{D_{cal}}} \right)}}.$$
(2)

In Eq. (2), the first term in the numerator is the proxy likelihood, the second term of the numerator is the prior distribution for temperature, and the denominator is a normalization constant. If we assume that a model exists that estimates the proxy value given the July–August temperature and a vector of hyperparameters (\(\tau\)), then the proxy likelihood can be written as:

$$\begin{aligned} p\left( {{D_i}|{t_i},~{D_{cal}},~{T_{cal}}} \right) &= \mathop \smallint \nolimits p\left( {{D_i},~\tau |{t_i},{D_{cal}}.~{T_{cal}}} \right)d\tau \\ &= \mathop \smallint \nolimits p\left( {{D_i},~{D_{cal}}|\tau ,{t_i},{T_{cal}}} \right)p\left( {\tau |{t_i},~{T_{cal}}} \right)/p({D_{cal}}|~{t_i},{T_{cal}})d\tau \\ &\propto \mathop \smallint \nolimits p\left( {{D_i}|{t_i},\tau } \right)p\left( {{D_{cal}}|{T_{cal}},\tau } \right)p(\tau )d\tau .\end{aligned}$$
(3)

On the second line, the denominator term \(p({D_{cal}}|~{t_i},{T_{cal}})\) appearing through the second application of Bayes theorem is again assumed constant. The third line includes, from left to right, the likelihood of the proxy datum \({D_i}\), the likelihood of the calibration data \({D_{cal}}\), and the prior over-the-proxy model hyperparameters. In a sense, both the model calibration and the model prediction were merged into one formula.

The posterior can be solved by making assumptions about the likelihoods and the priors. Here, we assume that the prior over July–August temperature \({t_i}\) is a normal distribution for which parameters are given by the moments of the previous STREC reconstruction (Summer Temperature Reconstruction for Eastern Canada) based only on ring width series (Gennaretti et al. 2014):

$$p\left( {{t_i}|{T_{cal}},~{D_{cal}}} \right) = p\left( {{t_i}} \right) \equiv N({t_i};{\mu _{strec}},~{\sigma _{strec}}).$$
(4)

Next, considering that we did not find evidence of non-linearity in the proxy-climate relationships, the model for computing the proxy likelihood can be defined by a linear regression with normally distributed errors:

$$p\left( {{D_i}|{t_i},~\tau } \right) \equiv N\left( {{D_i};~\mu = \alpha {t_i} + \beta ,~\sigma } \right),~\quad \tau \equiv ~\left( {\mu ,~\alpha ,\sigma } \right).$$
(5)

A prior for the used hyperparameters (\(\tau\)) also needs to be defined, and here, we chose non-informative priors for simplicity (Fig. S5): a uniform prior for the intercept, a Jeffreys prior for the variance and a prior for the slope that respects the invariance over the choice of dependent and independent variables (i.e., uniform prior in sin(tan − 1 α)). We thus obtained the following minimally informative prior on the models:

$$p\left( \tau \right) \equiv p(\alpha ,\beta ,\sigma ) \propto \frac{{{{(1 + {\alpha ^2})}^{( - 3/2)}}}}{\sigma }.$$
(6)

The last step is to substitute the generic data \(D\)with our three proxy datasets (\(R\) for ring width, \(O\) for δ18O and \(C\) for δ13C) and to introduce a set of hyperparameters specific to each proxy, denoted by \({\tau _R},~{\tau _O},~{\tau _C}.~\)We can now write the posterior of Eq. (2) as:

$$\begin{aligned} & p\left( {{t_i}|{R_i},~{R_{cal}},~{O_i},~{O_{cal}},~{C_i},~{C_{cal}},~{T_{cal}}} \right) \propto \mathop \smallint \nolimits p\left( {{R_i}|{t_i},~{\tau _R}} \right)p\left( {{R_{cal}}|{T_{cal}},~{\tau _R}} \right)p\left( {{\tau _R}} \right)d{\tau _R}\\ & \quad \times \mathop \smallint \nolimits p\left( {{O_i}|{t_i},~{\tau _O}} \right)p\left( {{O_{cal}}|{T_{cal}},~{\tau _O}} \right)p\left( {{\tau _O}} \right)d{\tau _O}\\ & \quad \times \mathop \smallint \nolimits p\left( {{C_i}|{t_i},~{\tau _C}} \right)p\left( {{C_{cal}}|{T_{cal}},~{\tau _C}} \right)p\left( {{\tau _C}} \right)d{\tau _C}p({t_i}).\end{aligned}$$
(7)

In practice, Eq. (7) is solved using Markov Chain Monte Carlo (MCMC) sampling with Metropolis–Hastings steps. Instead of sampling all 10 dimensions (three hyperparameter vectors plus t i) at once, each hyperparameter vector is sampled independently according to \(p\left( {{D_{cal}}|{T_{cal}},~\tau } \right)p\left( \tau \right)\) (Fig. S5). The samples are then used independently to obtain a temperature posterior distribution from each proxy or all together to obtain a sharper distribution using the proxy ensemble (Figs. S6, S7). The spread of the distribution of each of the proxy models is an indication of the weight (confidence) of each proxy.

To account for the fact that isotopic series were smoothed in the measurement process and their HFs were lost, the temperature LFs (\({t^{low}}\)) were reconstructed with the three proxies (ring width, δ18O and δ13C), whereas the temperature HFs (\({t^{high}}\)) were reconstructed with ring widths only. The LF and HF components are independent, each described by a model with its set of hyperparameters (Fig. S5). The posterior distributions of \({t^{low}}\) and \({t^{high}}\) for each year i were then combined in the final reconstruction.

The posterior probability of models with different combinations of LF proxies was also evaluated with the following equation (Kruschke 2014):

$$p\left( {m|~{T_{cal}},{D_{cal}}} \right) = \frac{{p\left( {{T_{cal}}|~m,{D_{cal}}} \right)p(m)}}{{\mathop \sum \nolimits_{m = 1}^7 p\left( {{T_{cal}}|~m,{D_{cal}}} \right)p(m)}}.$$
(8)

In Eq. (8), \(m\) is an indexal parameter specific to each of the 7 possible combinations of proxies (R, O, C, R&O, R&C, O&C, R&O&C), \(p\left( {m|~{T_{cal}},{D_{cal}}} \right)\) is the posterior probability of the proxy model, \(p\left( {{T_{cal}}|~m,{D_{cal}}} \right)\) is the likelihood of the data given the model, and \(p(m)\) is the prior probability of the model. Equation (8) is solved considering uniform model prior probabilities \((p(m) = 1/7)\), and evaluating the likelihood of the temperature calibration data with \(p\left( {{T_{cal}}|m,~{D_{cal}}} \right) = \mathop \smallint \nolimits p\left( {{T_{cal}}|m,~{D_{cal}},~{\tau _m}} \right)p\left( {{\tau _m}} \right)d{\tau _m}\). To be clear, if \(m = 4\) and the proxy considered are R and O, then \(p\left( {{T_{cal}}|m = 4,~{D_{cal}}} \right) = \mathop \smallint \nolimits p\left( {{T_{cal}}|{R_{cal}},~{\tau _R}} \right)p\left( {{\tau _R}} \right)d{\tau _R}p\left( {{T_{cal}}|{O_{cal}},~{\tau _O}} \right)p\left( {{\tau _O}} \right)d{\tau _O}\).

2.3 Impact of uncertainties in the proxy chronologies

The proposed Bayesian model evaluates the capacity of the proxies to reconstruct temperature values based on the relationship over the calibration period and considers uncertainties in model parameter estimation. There is no source of uncertainty that depends on time in the model. However, we implicitly evaluated the impact of the time-varying uncertainties in the LF proxy chronologies, which depend on the spread of individual series. All the used proxy chronologies are built with an almost stable sample depth over the last millennium. For each year we computed the probability density of the proxy chronology values assuming that the available replicates (three to six site-specific chronologies for ring width and cohorts of five trees for isotopes) are normally distributed (Fig. S8). These probability densities were sampled 100 times to create 100 chronologies per proxy to be included in the Bayesian model. The spread of the resultant 100 temperature reconstructions represent the impact of the uncertainties in the proxy chronologies.

3 Results and discussion

3.1 Final reconstruction

The final reconstruction (Three Proxies Summer Temperature Reconstruction for Eastern Canada, hereafter 3P-STREC) and its confidence intervals were derived from the median, 5th and 95th percentiles of the temperature posterior densities of each year (Fig. 1). The R2 of the LF 3P-STREC versus the temperature of the last century was 0.81. This is a substantial improvement relative to previous reconstructions with the same ring width data (STREC; R2 = 0.64; Gennaretti et al. 2014) or with the same δ18O data (i-STREC; R2 = 0.58 or 0.64 if considering mean maximal temperature values such as in Naulier et al. 2015). Note that here and hereafter, the comparisons with previous data were performed under the same conditions (using the same climate variable, frequency components, smoothing algorithm and 1905–2006 or 997–2006 periods). The added value of the multiproxy approach is clearly shown in Fig. 2. The posterior probability of the model with the three proxies was higher than the probability obtained with any other combination of one or two proxies. The Bayesian framework also allows sharp and reliable confidence intervals to be produced by leveraging the multiproxy information. Indeed, 95.1% of the observed temperature values were inside the 90% 3P-STREC nominal confidence intervals despite the spread of the 3P-STREC temperature posterior distributions being only 39% of the overall spread of the distributions obtained using the three proxies independently (Fig. 3a). The Bayesian model performed satisfactorily also over independent validation periods when the full period with temperature observations (1905–2006) was spit for a cross-calibration validation exercise (Table S1). The confidence intervals of 3P-STREC did not increase back in time because there is no source of uncertainty that depends on time in the model. However, the confidence intervals did not significantly vary even when we evaluated the impact of the time-varying uncertainties in the LF proxy chronologies (Fig. 3b). This is because the sample depths of the used chronologies were quite stable.

Fig. 1
figure 1

The entire 3P-STREC reconstruction (a; 997–2006) and a zoom over the last two centuries (b; 1800–2006) with low frequencies only (black bold lines) and with high frequencies added (black thin lines). The 90% confidence intervals are also shown (dark blue for the low frequency reconstruction and light blue for the final reconstruction). Red lines are mean July–August temperature values (CRU TS3.23; low frequencies only with bold lines and low plus high frequencies with thin lines)

Full size image
Fig. 2
figure 2

Posterior probability of models with different combinations of low frequency proxies. The models are ordered according to their probability. \(R\) for ring width, \(O\) for δ18O and \(C\) for δ13C

Full size image
Fig. 3
figure 3

Evaluation of uncertainties. a Reduction of reconstruction uncertainties with the proposed Bayesian methodology, which exploits the convergence of the three proxies. The figure shows the 90% confidence intervals of the low frequency 3P-STREC reconstruction (blue) and those derived from the sum of three independent temperature posterior densities for each year obtained from single proxies (light gray). b 90% confidence intervals of the reconstruction when considering an independent time-varying source of uncertainty due to individual series spread in each low frequency proxy chronology (dark gray)

Full size image

When the HFs (reconstructed with ring width only) were added to the LF 3P-STREC (Fig. 1), the R2 of the final reconstruction versus the temperature of the last century was 0.41. Again, this is a notable improvement over STREC (R2 = 0.30). This improvement mostly reflects the more realistic decadal and longer term signals achieved with the multiproxy approach. Indeed, HFs in 3P-STREC did not improve because they were only based on ring width data (R2 = 0.11 with HF July–August temperature values; Fig. S4) and were very similar to those of STREC. A coherency plot (Fig. S9) also showed that the fidelity of 3P-STREC is especially strong at medium (frequency = 0.33 corresponding to 3 year periodicities) and low frequencies (frequency < 0.1 corresponding to decadal and longer time-scales). For the final 3P-STREC, 92.16% of the observed temperature values were within the 90% nominal confidence intervals. However, the inter-annual variability was underestimated because ring width HFs have a low predictive power. The development of wood density chronologies has the potential to improve the HF results.

The previous reconstructions with the same data (STREC from ring width and i-STREC from δ18O) have suggested that the Medieval Climate Anomaly (MCA) was particularly warm in the Quebec taiga (Fig. 4a, b). The LF 3P-STREC confirmed this result, but it also suggested that the 15-year period between 1992 and 2006 was the warmest of the last millennium (Table S2). However, when HFs were added, some years during the Middle Ages and even during the Little Ice Age (LIA, 1300–1850) became as warm as the last decade (Fig. 1; Table S2). In addition, considering that the inter-annual variability of the final 3P-STREC was underestimated, more extreme warm years have probably occurred in the past. These findings are consistent with the fact that the MCA and the LIA were not necessarily synchronous everywhere in the World, and their intensity also varied spatially (Luterbacher et al. 2016; Mann et al. 2009). These spatial variations could be explained by regional climate feedbacks, such as those related to the Greenland sea ice and Labrador current dynamics in the case of the Quebec–Labrador peninsula (Miller et al. 2012; Schleussner and Feulner 2013; Stenchikov et al. 2009; Zanchettin et al. 2012; Zhong et al. 2011). The 3P-STREC reconstruction also allowed a better definition of the cold periods during the LIA in the Quebec taiga. According to i-STREC (Naulier et al. 2015), the coldest period in this region was between 1660 and 1700 during the Maunder solar minimum, while according to STREC (Gennaretti et al. 2014), the coldest period was between 1810 and 1860 coinciding with a part of the Dalton minimum and some active volcanic decades. The 3P-STREC reconstruction supported this second alternative that the first half of the nineteenth century was likely the coldest period of the LIA and of the last millennium in the Quebec taiga (Fig. 1). Another improvement achieved with the new Bayesian multiproxy reconstruction, was to confirm a more realistic long-term cooling trend over the last millennium (−0.75 ± 0.07 °C per 1000 years; estimate ± SE; computed on the final 3P-STREC). This trend, probably due to a combination of orbital changes (Kaufman et al. 2009; Esper et al. 2012) and volcanic activity (Miller et al. 2012), is more similar to that previously reconstructed with the same δ18O data (−0.52 ± 0.04 °C per 1000 years; i-STREC) than with the same ring width data (−1.60 ± 0.11 °C per 1000 years; STREC) and is quite consistent with values already published for nearby regions (Miller et al. 2013; Viau et al. 2012).

Fig. 4
figure 4

Comparison with previous and other up-to-date reconstructions. a Low frequency 3P-STREC (black) compared with i-STREC (blue, Naulier et al. 2015, reconstruction using oxygen isotopes only). b Comparison with STREC (red, Gennaretti et al. 2014, reconstruction using ring width only) and a STREC variant (gray) where the chronologies are built with the pivot correction RCS method (Autin et al. 2015). c Comparison with N-TREND summer temperature reconstruction (dark green; Wilson et al. 2016) based on 54 Northern Hemisphere tree-ring records including our ring width record. The median of all records in N-TREND excluding ours is shown in light green. d Comparison with the Northern Hemisphere mid-latitude temperature reconstruction based on 15 maximum latewood density chronologies (violet; Schneider et al. 2015). All plotted series are transformed in z-scores over common periods

Full size image

The concordance of 3P-STREC with the recent N-TREND Northern Hemisphere summer temperature reconstruction based on tree-ring width and density records (Wilson et al. 2016) is impressive if we exclude the last century (Fig. 4c), despite the very different spatial domain of these two reconstructions (a similar conclusion can also be drawn if the comparison is done with the reconstruction of Stoffel et al. 2015, as shown in Fig. S10). It has been suggested that few robust tree-ring chronologies along the northern North-Hemisphere treeline are sufficient to produce good estimates of reconstructions obtained with much larger networks of tree-ring data (D’arrigo and Jacoby 1993). Our study further suggests that a single regional reconstruction (maximum distance between sampling sites is 150 km), extremely well replicated at the site and the tree levels (six homogeneous sites and measurements from 1932 trees), may be very similar to an hemispheric one. This concordance can be explained by the following reasons: (1) the summer temperature variations of the northern Quebec taiga are consistent with the average hemispheric variations suggesting a significant influence of large scale global forcing (solar and volcanism) on regional climate; (2) tree growth is homogeneous across the Northern Hemisphere due to some large-scale climate factor forcing an overall synchronization regardless of the species. These two hypotheses are not mutually exclusive but an argument in favor of the second one is that N-TREND is more similar to the reconstruction based on ring width only than to the other reconstructions based on δ18O or δ13C (Fig. S11). This concordance with N-TREND also highlights the importance of synchronized coolings following volcanic eruptions at both the regional and hemispheric scale (see next section). The temperature of the last century in the northern Quebec taiga appears instead anomalously low with respect to the N-TREND reconstruction. This could probably be a characteristic of our study region in response to regional climate feedback or the result of a methodological problem in the recent part of the reconstructions. Nevertheless, the situation is different when we compared 3P-STREC with another recent Northern Hemisphere summer temperature reconstruction based on only maximum latewood density chronologies (Fig. 4d; Schneider et al. 2015). Some decadal-scale fluctuations mainly due to volcanic eruptions remained similar, but important differences were visible on the century-scale fluctuations, especially before the fifteenth century and during the twentieth century. This highlights the importance of better understanding the proxy behavior (here, ring width and isotopes versus density) and producing multiproxy reconstructions to reduce possible misrepresentations.

3.2 Proxy interpretation

Long tree ring chronologies built from a mix of living and dead trees are often standardized using the RCS approach in order to preserve low frequency variance. However, the resulting RCS chronologies are sensitive to various sampling procedures and data treatment approaches (Melvin et al. 2013; Matskovsky and Helama 2014). In particular, Autin et al. (2015) showed that ring width chronologies combining subfossil and living trees are prone to biases if they are built with common RCS techniques. These biases are generated when the samples originate from varying heights on the trees, and they are especially strong at the recent chronology end if old living trees are all sampled at the same height in contrast with what happens with subfossils. In the previous STREC reconstruction, living trees were standardized apart from subfossil stems in order to attenuate this sampling height problem; an approach also used elsewhere to combine heterogeneous subfossil and living tree datasets (Büntgen et al. 2010, 2011). In the present study all ring width chronologies were standardized using the “pivot correction”, a new variant of the RCS approach specifically designed to remove the sampling height bias in our material (Autin et al. 2015). When this method was used, the reconstruction from ring width over the last century became more similar to that from oxygen isotopes, suggesting that the previous STREC underestimated temperature during the twentieth century relative to Medieval times (Fig. 4b). Thus the ring width component included in the new 3P-STREC provides a better temperature reconstruction than the previous STREC, even if STREC is considered as a good, highly replicated dataset (Esper et al. 2016). This STREC versus 3P-STREC comparison also suggests that considering the sampling height issue elsewhere has the potential to improve RCS-based reconstructions using material with unknown or variable sampling heights. In spite of the pivot correction applied to our ring width series, the δ18O remained the most predictive LF proxy (Fig. 2, S4), generating more narrow temperature posterior densities (see Fig. S6).

Naulier et al. (2014) have previously shown that for black spruce trees of northeastern Canada, the summer mean and maximal temperatures were the most important climate parameters influencing variations of isotopic series (r = 0.50 and 0.39 for δ13C and r = 0.46 and 0.54 for δ18O, respectively). However, different climate and environmental signals can be carried by isotopic chronologies. The δ13C values depend on the ratio of leaf intercellular CO2 (ci) relative to the atmospheric CO2 pressure (ca), the δ13CCO2 values of ambient air (Farquhar et al. 1982), and soil moisture (Francey and Farquhar 1982). The ci/ca ratio is controlled by the photosynthetic capacity and the stomatal conductance, which itself depends on the water vapor deficit of the air (correlated with air temperature). In our case, we have already demonstrated that δ13C series are a good proxy for temperature reconstructions (Naulier et al. 2014). Gagen et al. (2011) have also found high correlations of their δ13C series with cloud cover and sunshine and with summer temperature. Young et al. (2010) showed that temperature and cloud cover may be in phase (positively correlated) or in opposition (negatively correlated), which explained the divergence between the δ13C series and temperature over the reconstructed 500-years period of their study. In northeastern Canada, we do not exclude the possibility that δ13C series could, additionally to temperature, be correlated to sunshine or cloud cover over the last millennium. However, sunshine series are not available for our study site and the short available cloud cover series (CRU) did not strongly correlate with our δ13C data.

The δ18O signal of tree-ring cellulose depends on the δ18O value of source water (soil). However, one of the main controls on the final δ18O values in tree-rings is the temperature prevailing regionally during cloud mass distillation, as registered in the raindrop signal and transferred to the source water in soils, then, through the root system, to the tree (McCarroll and Loader 2004). Additionally, δ18O values vary with the climatic factors influencing stomatal opening, such as temperature and moisture (McCarroll and Loader 2004; Naulier et al. 2014). For these reasons, δ18O series can show stronger correlation with summer temperature (in our case r = 46 and 0.54 for mean and maximal temperature, respectively, 1945–2005 period; Naulier et al. 2014, 2015) than with hydrological parameters like precipitation amounts (r = −0.41) or vapor pressure deficit (r = 0.44). The δ18O chronology used in this study was previously demonstrated to be a good proxy for temperature reconstruction (Naulier et al. 2014, 2015), a coherent finding owing to the fact that the studied area has a boreal climate not subjected to drought periods and that the collected trees were riparian and consequently, never under water stress.

In the end, it appears that ring widths, δ13C and δ18O chronologies have strengths and weaknesses as proxy of past temperature. However, they can generate complementary information, and when we used the three proxies together, the performance of the proposed Bayesian model increased (Fig. 2).

3.3 Proxy sensitivity to volcanic eruptions

The Bayesian framework can also be used to analyze the proxy differential sensitivity to climate forcings. Here we focused on the impact of volcanic eruptions because they produce abrupt temperature perturbations after known dates. Conversely, the impacts of the other forcings (solar and anthropogenic), although they are significant, are much more smoothed and difficult to isolate from other influences. In general, volcanic perturbations on temperature should be short lasting: 1–3 years for some authors (Fischer et al. 2007; Stoffel et al. 2015) or up to 10 years for others (Sigl et al. 2015). However, some uncertainties on these durations remain because the last century, having a robust observation network, is a period with relatively weak volcanism. Furthermore, some data and climate model experiments suggest that the impact of strong volcanic eruptions on temperature is much longer lasting if sustained by sea ice/ocean feedback, especially when two or more strong eruptions occur in close succession (Miller et al. 2012; Schleussner and Feulner 2013; Stenchikov et al. 2009; Zanchettin et al. 2012; Zhong et al. 2011). It was also hypothesized that a series of strong volcanic eruptions during the twelfth and thirteenth centuries triggered the onset of the LIA in the eastern Canadian Arctic and in northeastern North America (Miller et al. 2012). These suggestions appear to be supported by the volcano-induced regime shifts that we found in our ring width chronologies (Gennaretti et al. 2014). However, using tree-ring proxy data and especially ring width, it is hard to discriminate the volcanic impact on temperature (i.e., the object of the reconstruction) from other collateral volcanic impacts, including tree damage, reductions in solar irradiance and changes in diffuse radiation (Robock 2005). In particular, it is considered that ring width series show a smeared delayed response to volcanic eruptions in comparison with instrumental temperature values and tree-ring density data due to their greater autocorrelation (Esper et al. 2013, 2015). Consequently, ring width series are considered inappropriate to examine the response to volcanism at interannual scale (D’Arrigo et al. 2013). Similarly, our smoothed isotopic data do not seem to be appropriate. However, as STREC and 3P-STREC contain a low-frequency signal linked to volcanism in response to very strong eruptions (an influence that extends up to 20 years in some cases), it is interesting to compare the volcanic signal between the ring width and isotopic components of 3P-STREC, even if this comparison would be based on smoothed data (chronologies equivalent to series smoothed with a 9-year triangular filter).

To verify whether the isotope ratios show the same sensitivity to volcanic eruptions as the ring width, we first compared the responses from their independent LF reconstructions (i.e., based on a single proxy) to each of the 46 strong volcanic eruptions listed in Table 1 (Fig. 5). These eruptions were derived from Gao et al. (2008), Crowley and Unterman (2013), Esper et al. (2013) or Sigl et al. (2013, 2015). In the figure, the y-axis shows the subtraction between the post-eruption anomalies reconstructed from isotope ratios (average anomalies from δ18O and δ13C) and from ring width. Positive (negative) values indicate that post-eruption anomalies are more negative in the reconstruction from ring width (isotope ratios). We can thus see that ring width indexes often reconstruct more negative temperature anomalies after volcanic events (Fig. 5b) than isotopic proxies, confirming the different proxy sensitivity to volcanic forcing. However, the variability of responses after different eruptions is large. Indeed different volcano locations, atmospheric states during the eruptions, event intensities and days of the year may cause specific impacts over the Quebec-Labrador peninsula, which are more or less registered by the proxies.

Table 1 List of 46 eruption years derived from Gao et al. (2008; global total stratospheric sulfate aerosol injection > 15Tg), Crowley and Unterman (2013; satellite aerosol optical depth > 0.1), Esper et al. (2013; eruptions classified as large events) or Sigl et al. (2013, 2015; eruptions classified as large events)
Full size table
Fig. 5
figure 5

Comparison of responses from low frequency reconstructions based on ring width and isotopic proxies to individual strong volcanic eruptions of the last millennium (listed in Table 1). Plot a shows if the anomalies of each 15 post-eruption years (anomalies relative to the year before the eruption) are more negative in the reconstruction from ring width (positive values) or from stable isotope ratios (negative values; average anomalies from δ18O and δ13C). Eruption years are vertical lines. Plot b shows a Superimposed Epoch Analysis (median and 60% confidence intervals) based on the data and the eruption years of plot a

Full size image

Subsequently, we focused on the three periods of the last millennium where the volcanic impact should have been the strongest and for which Gennaretti et al. (2014) found significant regime shifts toward lower growth values in the ring width chronologies. These periods are the thirteenth century with the series of eruptions centered around the 1258 Samalas event and that supposedly triggered the onset of the LIA; the second half of the fifteenth century with the cooling episode following the 1458–1459 Kuwae eruption; and the first half of the nineteenth century, probably the coldest period of the last millennium in the Quebec taiga, with a series of eruptions centered around the 1815 Tambora event. This analysis (Fig. 6) reveals that δ13C of spruce trees growing on lakeshores is almost insensitive to volcanic impacts. It is also clear that ring width shows much more important volcanic impacts than δ18O, reconstructing the strongest coolings after single events and the strongest cumulative impacts from multiple events. Interestingly, during the thirteenth century at the onset of the LIA, all three proxies agree in reconstructing an overall cooling trend despite the amplitude of the trend and short-term temperature variations depending on the proxies. In this case, the impact of individual eruptions appears not correlated to their intensity. Indeed, following the 1258 Samalas eruption, which was likely the strongest of the last millennium (Lavigne et al. 2013), the proxies react as after other strong eruptions of the same period. During the second half of the fifteenth century, δ18O and ring width strongly disagree in the volcanic responses. While δ18O shows only moderate coolings (approximately −0.5 °C) lasting for no more than 3 years after the two strong eruptions of 1452 and 1458–1459, ring width displays an important cumulative impact of these two events, culminating in anomalies of approximately −3.5 °C (note that this discussion on cooling intensities and durations is based on smoothed data because only low frequencies are available for isotope ratios due to the data processing method). A good proxy agreement is instead observed during the first half of the nineteenth century where the decadal trends in δ18O and ring width are very similar despite ring width again showing larger temperature variations. The LF 3P-STREC reconstruction with the ensemble of the proxies always has an intermediate behavior between the δ18O and the ring width reconstructions with a preference towards δ18O which is likely the better estimate of past temperatures.

Fig. 6
figure 6

Proxy response during periods of strong volcanic impact: the thirteenth century (a), the second half of the fifteenth century (b) and the first half of the nineteenth century (c). Temperature anomalies with respect to the first considered year are shown on the left of each panel (medians and 90% confidence intervals), and anomaly posterior densities of the year with lowest anomalies are shown on the right (the selected year is highlighted by a vertical gray bar on the left plots). Black colors for the proxy ensemble low frequency reconstruction 3P-STREC, red for the reconstruction using ring width only, blue for δ18O and green for δ13C. Vertical dotted lines and arrows indicate eruption years listed in Table 1

Full size image

As this comparison of the proxy differential sensitivity to volcanic eruptions showed that ring width LFs responded more strongly than isotopic ones, it is possible that the ring width component of 3P-STREC amplifies the volcanic forcing after strong eruptions (see Büntgen et al. 2015). This may be the case especially in cold periods, such as the Little Ice Age, and/or if the time between successive eruptions is short (Sigl et al. 2013). Indeed, in cold periods with recurrent volcano-induced temperature reductions, tree growth is likely to suffer from the depletion of carbohydrate reserves and from the fact that the tissues and canopy can be frost-damaged. Conversely, it is also possible that the isotopic proxies are less sensitive to the meteorological perturbations related to volcanic events. Latewood maximum density data, which nevertheless could also have their limitations (see Stine and Huybers 2014; Tingley et al. 2014), may be useful to discriminate between these hypotheses.

4 Conclusion

In this paper, we presented the first annually resolved and millennium-long multiproxy temperature reconstruction for northeastern North America based on tree-ring width and stable isotope ratios (δ13C and δ18O). The climate information embedded in the three proxies was exploited using a Bayesian framework that allows for a rigorous uncertainty assessment. The results showed that the final reconstruction 3P-STREC was a clear improvement over previous attempts based only on ring width or δ18O values. Using the Bayesian methodology, we can also provide reliable confidence intervals that are much sharper than simply merging the three independent reconstructions from single proxies. At the moment, 3P-STREC represents the best estimate of the past summer temperature of our region (the mean of July–August temperature), but in the future, the results could be further improved by (1) adding additional proxies (e.g., ring maximum density) to reduce the remaining sources of LF errors and improve the HF calibration, (2) considering more complex mechanistic modeling of the proxy-climate relationships, (3) providing a more robust intra- and inter-site-explicit uncertainty analysis, or (4) taking into account the year to year memory of each component in the model as a fractional Gaussian process (Lovejoy et al. 2015) in order to properly integrate the information from different sources (Li et al. 2010) and investigate the impact of higher long-term persistence in ring width data compared to instrumental data (Zhang et al. 2015). However, concerning this last point, we want to stress out that 3P-STREC is already shown to be well calibrated in the LF domain.

We also examined different proxy sensitivities to climate forcing, focusing the analysis on the impact of volcanic eruptions. The reconstruction from ring width had a larger and longer response to single or consecutive eruptions, while those from isotope ratios showed an intermediate (δ18O) or nearly absent (δ13C) volcanic impact. It is difficult to know the real influence of past volcanic events on temperature at the regional scale because many uncertainties exist due to (1) the limited knowledge on the mechanisms behind the proxy responses, (2) the absence of major volcanic events during the period when the proxies are compared with instrumental data, (3) the low agreement in the volcanic datasets concerning timing and respective forcing of eruptions (Crowley and Unterman 2013; Esper et al. 2013; Gao et al. 2008; Sigl et al. 2013, 2015), and (4) the complexity of comparing the results of proxy-based reconstructions with simulations of past climate with plausible initial conditions, sulfate aerosols and ocean feedback (circulation, sea ice, relaxation time of subsurface temperature and sea level) during volcanic coolings (Stenchikov et al. 2009). In our reconstruction, the uncertainties came from the specific reduced or amplified responses to volcanic eruptions of the three temperature-sensitive proxies. We should bear in mind that each proxy has its advantages and drawbacks. Only a multi-proxy approach can allow us to take advantage of the proxy convergence and reduce individual sources of error caused by specific mechanistic responses of the proxies.