Bayesian multiproxy temperature reconstruction with black spruce ring widths and stable isotopes from the northern Quebec taiga

TLDR: In this article, a tree-ring dataset for temperature reconstructions over the last millennium was developed in the northern Quebec taiga, which is composed of one δ18O and six ring width chronologies.

Abstract: Northeastern North America has very few millennium-long, high-resolution climate proxy records. However, very recently, a new tree-ring dataset suitable for temperature reconstructions over the last millennium was developed in the northern Quebec taiga. This dataset is composed of one δ18O and six ring width chronologies. Until now, these chronologies have only been used in independent temperature reconstructions (from δ18O or ring width) showing some differences. Here, we added to the dataset a δ13C chronology and developed a significantly improved millennium-long multiproxy reconstruction (997–2006 CE) accounting for uncertainties with a Bayesian approach that evaluates the likelihood of each proxy model. We also undertook a methodological sensitivity analysis to assess the different responses of each proxy to abrupt forcings such as strong volcanic eruptions. Ring width showed a larger response to single eruptions and a larger cumulative impact of multiple eruptions during active volcanic periods, δ18O showed intermediate responses, and δ13C was mostly insensitive to volcanic eruptions. We conclude that all reconstructions based on a single proxy can be misleading because of the possible reduced or amplified responses to specific forcing agents.

Full text

Vol.:(0123456789)
1 3
Clim Dyn (2017) 49:4107–4119
DOI 10.1007/s00382-017-3565-5
Bayesian multiproxy temperature reconstruction withblack
spruce ring widths andstable isotopes fromthenorthern Quebec
taiga
FabioGennaretti
1
· DavidHuard
2
· MaudNaulier
3
· MartineSavard
4
·
ChristianBégin
4
· DominiqueArseneault
5
· JoelGuiot
1

Received: 21 July 2016 / Accepted: 31 January 2017 / Published online: 1 March 2017
© The Author(s) 2017. This article is published with open access at Springerlink.com
Ring width showed a larger response to single eruptions
and a larger cumulative impact of multiple eruptions dur-
ing active volcanic periods, δ
18
O showed intermediate
responses, and δ
13
C was mostly insensitive to volcanic
eruptions. We conclude that all reconstructions based on
a single proxy can be misleading because of the possible
reduced or amplified responses to specific forcing agents.
Keywords Tree-ring· Oxygen isotopes· Carbon
isotopes· Last millennium· Summer temperature·
Volcanic impact· Proxy sensitivity
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. How-
ever, these chronologies are much more abundant in the
Abstract Northeastern North America has very few
millennium-long, high-resolution climate proxy records.
However, very recently, a new tree-ring dataset suitable
for temperature reconstructions over the last millennium
was developed in the northern Quebec taiga. This data-
set is composed of one δ
18
O and six ring width chronolo-
gies. Until now, these chronologies have only been used
in independent temperature reconstructions (from δ
18
O
or ring width) showing some differences. Here, we added
to the dataset a δ
13
C chronology and developed a signifi-
cantly improved millennium-long multiproxy reconstruc-
tion (997–2006 CE) accounting for uncertainties with a
Bayesian approach that evaluates the likelihood of each
proxy model. We also undertook a methodological sen-
sitivity analysis to assess the different responses of each
proxy to abrupt forcings such as strong volcanic eruptions.
Electronic supplementary material The online version of this
article (doi:
10.1007/s00382-017-3565-5) contains supplementary
material, which is available to authorized users.
* Fabio Gennaretti
gennaretti@cerege.fr
David Huard
huard.david@ouranos.ca
Maud Naulier
maud.naulier@irsn.fr
Martine Savard
martinem.savard@canada.ca
Christian Bégin
christian.begin@canada.ca
Dominique Arseneault
Dominique_Arseneault@uqar.ca
Joel Guiot
guiot@cerege.fr
1
Aix Marseille Univ, CNRS, IRD, Coll France, CEREGE,
13545Aix-en-Provence, France
2
Ouranos Consortium, 550 Rue Sherbrooke O,
H3A1B9Montréal, Canada
3
Institut de Radioprotection et de Sûreté Nucléaire, CEA
Cadarache, 13108Saint-Paul-Lez-Durance, France
4
Geological Survey ofCanada, Natural Resources Canada,
490 Rue de la Couronne, G1K9A9Québec, Canada
5
Département de biologie, chimie et géographie, Université
du Québec à Rimouski, 300 allée des Ursulines,
G5L3A1Rimouski, Canada

4108 F.Gennaretti et al.
1 3
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://past-
globalchanges.org/ini/wg/2k-network/intro). In northeast-
ern 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 qual-
ity 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 prox-
ies 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] subfos-
sil trees preserved in six lakes. Samples from one of these
sites were also used to obtain a millennium-long record of
oxygen isotopic ratios (δ
18
O) in tree-ring cellulose (Naulier
etal. 2015). Until now, these ring width and δ
18
O data have
only been used separately in independent summer tempera-
ture 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 Bayes-
ian approach, allowing a better understanding of the tem-
poral 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) bet-
ter 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 δ
13
C
chronology in tree-ring cellulose, which enhances the
temperature-sensitive proxy dataset from the northern Que-
bec taiga. Thus, we use an ensemble of three proxies (ring
width, δ
18
O and δ
13
C) to provide the first multiproxy regional
high-resolution summer temperature reconstruction in north-
eastern 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 evalu-
ated and discussed, focusing the analysis on the response to
strong volcanic eruptions because they produce abrupt per-
turbations after key-dates.
2 Materials andmethods
2.1 Proxy andclimate data
The tree-ring data from the northern Quebec taiga are com-
posed of six millennium-long ring width chronologies devel-
oped 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 chronolo-
gies were built with the regional curve standardization (RCS)
pivot correction method to reduce the impact of varying sam-
pling heights (Autin etal. 2015). Low and high frequencies
of the median chronology were treated separately. Low fre-
quencies (LFs) were obtained with a 9-year triangular filter
to produce a chronology comparable to that of the stable iso-
topes (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 trans-
formation 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):
where erf represents the Gauss error function,
P
(
y
)
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)
(1)
y
�
i
=
√
2erf
−1
(2P(y
i
)−1)
,

4109Bayesian multiproxy temperature reconstruction withblack spruce ring widths andstable…
1 3
shows the effect (quite low in this case) of this transforma-
tion on the LF chronology.
From one of the six sites from the northern Quebec
taiga, 60 subfossil logs and 5 living trees were further ana-
lyzed at the Delta-lab of the Geological Survey of Canada
to extract two millennium-long chronologies of stable iso-
tope ratios in tree-ring cellulose. The δ
18
O chronology is
already in the public domain (Naulier etal. 2015), whereas
the δ
13
C chronology is detailed here and is accessible in
the supplementary material (Dataset S1). These chro-
nologies 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 sto-
chastic process, this method produces chronologies equiva-
lent 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 δ
13
C values of the modern
part of the chronology were also mathematically corrected
for atmospheric δ
13
C CO
2
changes due to fossil fuel com-
bustion (Suess effect), and for plant response to increasing
isotopes. The relationships between the temperature vari-
able 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 distribu-
tion of temperature t at each year i can be obtained by:
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 nor-
malization constant. If we assume that a model exists that
estimates the proxy value given the July–August tempera-
ture and a vector of hyperparameters (
𝜏
), then the proxy
likelihood can be written as:
(2)
p
(
t
i
|D
i
, D
cal
, T
cal
)
=
p
(
D
i
|t
i
, D
cal
, T
cal
)
p
(
t
i
|T
cal
, D
cal
)
p
(
D
i
|
T
cal
, D
cal
)
.
atmospheric CO
2
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 chronol-
ogy, 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 tem-
perature 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 previ-
ous 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 tempera-
ture 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
(3)
p
(
D
i
|t
i
, D
cal
, T
cal
)
=
∫
p
(
D
i
, 𝜏 |t
i
, D
cal
. T
cal
)
d𝜏
=
∫
p
(
D
i
, D
cal
|𝜏 , t
i
, T
cal
)
p
(
𝜏 |t
i
, T
cal
)
∕p(D
cal
| t
i
, T
cal
)d
𝜏
∝
∫
p
(
D
i
|
t
i
, 𝜏
)
p
(
D
cal
|
T
cal
, 𝜏
)
p(𝜏 )d𝜏 .
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 calibra-
tion 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 dis-
tribution for which parameters are given by the moments
of the previous STREC reconstruction (Summer Tempera-
ture Reconstruction for Eastern Canada) based only on ring
width series (Gennaretti etal. 2014):
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:
(4)
p
(
t
i
|
T
cal
, D
cal
)
= p
(
t
i
)
≡ N(t
i
;𝜇
strec
, 𝜎
strec
)
.

4110 F.Gennaretti et al.
1 3
A prior for the used hyperparameters (
𝜏
) also needs to be
defined, and here, we chose non-informative priors for sim-
plicity (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 independ-
ent variables (i.e., uniform prior in sin(tan
− 1
α)). We thus
obtained the following minimally informative prior on the
models:
The last step is to substitute the generic data
D
with
our three proxy datasets (
R
for ring width,
O
for δ
18
O
and
C
for δ
13
C) and to introduce a set of hyperparame-
ters specific to each proxy, denoted by
𝜏
R
, 𝜏
O
, 𝜏
C
.
We can
now write the posterior of Eq.(2) as:
In practice, Eq. (7) is solved using Markov Chain
Monte Carlo (MCMC) sampling with Metropolis–Hast-
ings steps. Instead of sampling all 10 dimensions (three
hyperparameter vectors plus t
i
) at once, each hyperpa-
rameter vector is sampled independently according to
p
(
D
cal
|
T
cal
, 𝜏
)
p
(
𝜏
)
(Fig. S5). The samples are then used
independently to obtain a temperature posterior distribu-
tion 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 mod-
els 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, δ
18
O and δ
13
C), whereas
the temperature HFs (
t
high
) were reconstructed with ring
widths only. The LF and HF components are independ-
ent, each described by a model with its set of hyperpa-
rameters (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):
(5)
p
(
D
i
|
t
i
, 𝜏
)
≡ N
(
D
i
; 𝜇 = 𝛼t
i
+ 𝛽, 𝜎
)
, 𝜏 ≡
(
𝜇, 𝛼, 𝜎
)
.
(6)
p
(
𝜏
)
≡ p(𝛼, 𝛽, 𝜎)∝
(1 + 𝛼
2
)
(−
3
∕
2
)
𝜎
.
(7)
p
(
t
i
|R
i
, R
cal
, O
i
, O
cal
, C
i
, C
cal
, T
cal
)
∝
∫
p
(
R
i
|t
i
, 𝜏
R
)
p
(
R
cal
|T
cal
, 𝜏
R
)
p
(
𝜏
R
)
d𝜏
R
×
∫
p
(
O
i
|t
i
, 𝜏
O
)
p
(
O
cal
|T
cal
, 𝜏
O
)
p
(
𝜏
O
)
d𝜏
O
×
∫
p
(
C
i
|
t
i
, 𝜏
C
)
p
(
C
cal
|
T
cal
, 𝜏
C
)
p
(
𝜏
C
)
d𝜏
C
p(t
i
).
(8)
p

m T
cal
, D
cal

=
p

T
cal
 m, D
cal

p(m)

7
m=1
p

T
cal

m, D
cal

p(m)
.
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
(
m
|
T
cal
, D
cal
)
is the posterior probability of the proxy model,
p
(
T
cal
|
m, D
cal
)
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
(
T
cal
|
m, D
cal
)
=
∫
p
(
T
cal
|
m, D
cal
, 𝜏
m
)
p
(
𝜏
m
)
d𝜏
m
. To be
clear, if
m = 4
and the proxy considered are R and O, then
p
(
T
cal
|
m = 4, D
cal
)
=
∫
p
(
T
cal
|
R
cal
, 𝜏
R
)
p
(
𝜏
R
)
d𝜏
R
p
(
T
cal
|
O
cal
, 𝜏
O
)
p
(
𝜏
O
)
d𝜏
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 consid-
ers 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 chronol-
ogies, 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 Tem-
perature Reconstruction for Eastern Canada, here-
after 3P-STREC) and its confidence intervals were
derived from the median, 5th and 95th percentiles of

4111Bayesian multiproxy temperature reconstruction withblack spruce ring widths andstable…
1 3
the temperature posterior densities of each year (Fig.1).
The R
2
of the LF 3P-STREC versus the temperature of
the last century was 0.81. This is a substantial improve-
ment relative to previous reconstructions with the same
ring width data (STREC; R
2
= 0.64; Gennaretti et al.
2014) or with the same δ
18
O data (i-STREC; R
2
= 0.58
or 0.64 if considering mean maximal temperature values
such as in Naulier etal. 2015). Note that here and hereaf-
ter, the comparisons with previous data were performed
under the same conditions (using the same climate vari-
able, 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 frame-
work 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 inter-
vals despite the spread of the 3P-STREC temperature
posterior distributions being only 39% of the overall
spread of the distributions obtained using the three prox-
ies independently (Fig. 3a). The Bayesian model per-
formed satisfactorily also over independent validation
periods when the full period with temperature observa-
tions (1905–2006) was spit for a cross-calibration vali-
dation 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 signif-
icantly vary even when we evaluated the impact of the
Year
1800 1850 1900 1950 2000
810121416
Te mp. (°C)
b
1000 1200 1400 1600 1800 2000
81012141
6
Te mp. (°C)
a
Fig. 1 The entire 3P-STREC reconstruction (a; 997–2006) and a
zoom over the last two centuries (b; 1800–2006) with low frequen-
cies 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 recon-
struction). Red lines are mean July–August temperature values (CRU
TS3.23; low frequencies only with bold lines and low plus high fre-
quencies with thin lines)
Proxy combination
Model poste
rior probability
CRO R_C O_C R_O R_O_C
0.00 0.05 0.10 0.15
Fig. 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 δ
18
O and
C
for δ
13
C
1000 1200 1400 1600 1800 2000
810121416
Te mp. (°C)
Year
b
1000 1200 1400 1600 1800 2000
810121416
Te mp. (°C)
a
Fig. 3 Evaluation of uncertainties. a Reduction of reconstruc-
tion 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 recon-
struction (blue) and those derived from the sum of three independ-
ent 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 chronol-
ogy (dark gray)