Volcano-induced regime shifts in millennial tree-ring chronologies from northeastern North America

Our summer temperature reconstruction for Eastern Canada (hereafter STREC) closely reproduces the warming trend of July−August temperatures over the study area during the last century (Fig. 1A), even when the region experienced one of the strongest increases in summer temperatures worldwide according to the gridded temperature dataset Climate Research Unit (CRU) TS3.20 (26). The correlation between STREC and the gridded July−August temperature data (1905−2011; CRU TS3.20) is highly significant (0.61, P < 0.001), even when the instrumental period is split into two subsets (1905−1957 and 1958–2011; Fig. 1A). The skill of the reconstruction in reproducing temperature data is indicated by the positive values of the reduction-of-error statistics computed with different calibration periods and on unsmoothed and smoothed datasets (Table S1). Cross-calibration verification confirms that STREC is robust and better predicts low frequencies than high frequencies (Table S1). In our study, the ring-width series of the subfossil and living tree samples were homogenized to remove the bias caused by varying sampling heights (see Methods) and avoid divergence between temperatures and standardized tree growth (Figs. S1 and S2). Spatial correlations show that STREC is mostly valid over the central Quebec-Labrador Peninsula within the NENA sector (Fig. 2).

Volcanism has been the primary factor forcing changes in summer temperatures on the decadal time scale in Eastern Canada during the last millennium. Numerous cooling episodes in STREC are synchronous with episodes already reported in hemispheric temperature reconstructions and simulations in response to strong volcanic eruptions (Fig. S3). The agreement between peaks in the global stratospheric volcanic sulfate aerosol injections (9) and cold anomalies inferred by STREC is also striking (Fig. 3A). A Superimposed Epoch Analysis (SI Methods) demonstrates that the 10 strongest volcanic eruptions of the last millennium produced highly significant cooling episodes in Eastern Canada lasting for about two decades, while less intense volcanic eruptions had a shorter influence (Fig. 4). The 20 postevent summers were significantly colder than the preceding ones for 8 out of the 10 largest eruptions (P < 0.1; Table S2). For example, temperature anomalies ranged from −1.3 °C to −3.0 °C in response to the three strongest tropical eruptions of the last millennium (A.D. 1257, 1452/1453, and 1815). With about two cooling episodes per century lasting for 10–20 y with a consequent reduction in tree growth, the volcanic signature in the NENA taiga is comparable to the epidemic signal of the eastern spruce budworm (Choristoneura fumiferana), which is the most destructive insect in the commercial forest of Eastern Canada south of our study area (27). Therefore, volcanic forcing during the last millennium has clearly impacted net primary production and the carbon balance in the NENA forests, at least over the Quebec-Labrador taiga.

Volcanism also strongly influenced century-scale temperature variability in NENA, as shown by our Bayesian analysis of regime shifts in STREC (SI Methods). Two of the strongest eruptions of the last 1,100 y, the A.D. 1257 Samalas and A.D. 1815 Tambora events, which were followed and preceded by other eruptions, coincide exactly with the two most persistent regime shifts detected in STREC. According to the Schwarz criterion (28), STREC is best modeled with a four-state normal Bayesian hidden Markov model, with the two warmer and the two colder regimes largely dominating the A.D. 910–1257 and A.D. 1816–2011 time periods, respectively (Fig. 3B and Fig. S4). Other studies based on tree-ring data have already reported temperature reductions post-A.D. 1257 and post-A.D. 1815 in other regions of the world (see, for example, refs. 29 and 30), but with smaller amplitude and temporal extent than STREC. Indeed, the strong overall cooling trend obtained by fitting a linear regression model to STREC (−1.60 ± 0.11 °C per 1,000 y; estimate ± SE; P < 0.001) can be mostly attributed to the A.D. 1257 and A.D. 1815 shifts. Some proxy-based evidence has also shown that a long-term cooling trend due to orbital changes has characterized the climate of the last 2,000 y in many regions of the world (4, 31, 32), including NENA (33). However, this cooling trend has been estimated at about −0.3 °C to −0.5 °C per 1,000 y (4, 32, 33) and is much weaker than the one of STREC. These facts, in addition to the absence of any negative trends in STREC over the A.D. 1257–1815 time period (Fig. S5), highlight the importance of volcanic-induced temperature regime shifts in NENA. Along with volcanism and orbital forcing, additional factors may have contributed to the strong negative trend of STREC and the associated transition from warm to cold regimes (Fig. 3B and Figs. S4 and S5), including solar forcing (19), and the specific regional domain of STREC. Because STREC is based on several sites and because of our data management approach (see Methods), the A.D. 1257 and A.D. 1815 shifts do not seem to be influenced by local nonclimatic disturbances affecting our chronologies (Fig. S6).

A well-expressed Medieval Climate Anomaly (A.D. 910–1257) occurred in NENA before the A.D. 1257 Samalas event. The warmest decades reconstructed by STREC occurred between A.D. 1141 and 1170 (positive anomalies ranging from 0.89 °C to 1.80 °C with respect to the last decade) and between A.D. 1061 and 1095 (0.87–1.19 °C; Table S3). The confidence intervals of STREC for these two periods were almost all higher than the mean temperatures of the last decade (Fig. 1B). The amplitude and timing of the Medieval Climate Anomaly reconstructed by STREC also resemble the results of a pollen-based temperature reconstruction for the North American forest tundra (5), and closely correspond to a period of ice-cap melting in the Eastern Canadian Arctic (22) (Fig. 3A). Collectively, these complementary data sources demonstrate that a major and prolonged climatic shift occurred over the NENA sector after a series of 13th century volcanic eruptions centered around the A.D. 1257 Samalas event. This shift marked the end of the Medieval Climate Anomaly and the beginning of the Little Ice Age in this sector.

Similarly, the series of eruptions centered around the A.D. 1815 Tambora event shifted summer temperatures to the coldest 40-y period of the last 1,100 y in NENA. The A.D. 1815–1857 episode was extremely cold in Eastern Canada, with decadal anomalies reconstructed by STREC ranging from −2.76 °C to −3.58 °C relative to the last decade (Table S3). Low solar activity during this period, that started during the Dalton Minimum, could have concurred to cause the cooling episode (34). However, the regime shift observed in STREC coincides with the A.D. 1815 Tambora eruption (Fig. 3B and Fig. S4). The same climate shift has already been observed in a 263-y tree-ring-based summer temperature reconstruction whose sampling sites are 400 km eastward from our study sites (7), whereas a pollen-based temperature reconstruction suggests that this period was the coldest of the last 2,000 y in the North American boreal forest and forest−tundra (5). This climate shift was probably limited to the NENA sector as STREC diverges from the ensemble of Northern Hemisphere reconstructions and simulations after A.D. 1816 (Fig. S3 C and D).

Although the last century reconstructed by STREC was warmer than the early 19th century, it was colder than the Medieval Climate Anomaly (the difference between the average summer temperature of the 12th century and of the last 100 y is 1.66 °C, P < 0.001, according to the one-tailed Wilcoxon rank−sum test). In fact, STREC shows that the NENA sector experienced relatively cold conditions until late into the 20th century. This persistence of cold conditions over NENA is also suggested by the strong warming trend of the last 100 y denoting a colder starting point (Fig. 1 A and B), as well as by permafrost growth during the mid-20th century on the southern shore of Hudson Strait (35) and by the lack of postfire forest recovery over the last 900 y at the northern Quebec treeline in contrast to what occurred during the Medieval Climate Anomaly (36, 37). The warming trend in our study area has accelerated over the last 30 y (+0.7 °C per decade according to the dataset CRU TS3.20). If this trend continues, then summer temperatures will be similar to the maximum of the last 750 y during the next decade and to the maximum of the last 1,100 y during the following one (based on STREC and the data of Table S3).

Several hypotheses have been suggested to explain why posteruption temperature anomalies reconstructed from tree-ring data generally express higher values than expected, including regional variations in response to volcanic events, autocorrelation in ring-width series (17, 18), failure of growth rings to form during volcano-induced cold summers (20), and increased tree growth caused by volcano-induced diffuse radiation (38). In this study, we observed a strong response of ring-width data to volcanic activity, with the amplitudes and duration of negative anomalies (including persistent regime shifts) similar to model predictions. Diagnostic light rings are frequent in black spruce trees of our study area (39), thus allowing for a rigorous control of ring dating and for identifying the occurrence of missing rings. In addition, although autocorrelation in ring-width chronologies of black spruce is high (Table 1), the results of a simple model show that its effects on the amplitude and duration of reconstructed negative temperature anomalies after volcanic events are low (Fig. S7). Instead, our results suggest that the climatic impacts of eruptions vary among regions of the Northern Hemisphere and that the NENA sector is especially sensitive to these impacts compared with Eurasia, where the majority of temperature-sensitive tree-ring chronologies have been previously developed. This idea is also supported by recent simulation experiments (12–16, 22), which show that large and successive eruptions may trigger cold episodes whose duration may be sustained by complex and variable sea−ice−ocean feedbacks in the North Atlantic and that the resultant northward heat transport would tend to be more severely attenuated in the NENA than in the Eurasian sectors.

Our sampling area is situated in the Eastern Canadian taiga between latitudes 53.8°N and 54.6°N and longitudes 70.2°W and 72.5°W (Table 1). According to the gridded temperature dataset CRU TS3.20 (26), this region has experienced one of the fastest temperature increases on Earth during the last century. The mean July−August temperature, which is the object of our reconstruction, has increased by an average of 0.68 ± 0.15 °C and 0.19 ± 0.02 °C (estimate ± SE) each 10 y during the last 30 and the last 110 y, respectively.

To implement our network of tree-ring chronologies, 1782 black spruce [Picea mariana (Mill.) B.S.P.] subfossil trees were sampled from six lakes of the study area and cross-dated to the calendar year. Particular care was taken in selecting and replicating sites and trees to construct a dataset suitable for Regional Curve Standardization (RCS). The RCS method preserves long-term climate trends in tree-ring chronologies built with short-lived species, but requires high replication of trees belonging to a homogeneous population not disturbed by external factors (40). The selected sites were all characterized by an old-growth riparian forest on the side of the lake protected from dominant winds, an abrupt forest−lake transition, large stocks of subfossil trees in the littoral zone, and a well-documented low fire recurrence (23, 41). This allowed the development of an exceptional network of climate-sensitive tree-ring data comprising six local chronologies in the same region (one per lake), each including from 75 to 586 subfossil trees and spanning from 1,238 to 1,440 y (Table 1). The living trees extending the chronologies to the present day (25 trees per site) were sampled at the same sites as the subfossil samples, and issues related to tree selection and sampling height were also considered (see below).

To produce the six local chronologies (Fig. S3), two radii were measured at a precision of 1/100 mm and then averaged for each tree. A spline with a 50% frequency cutoff and a time-varying response [starting at 10 y and increasing by one each year (42)] was used as a smoothing algorithm to generate the local growth curves needed for RCS standardization. All specimens had well-preserved piths, and pith offset could easily be accounted for in the chronology development.

To guarantee homogeneity between subfossil and living samples, 25 living trees at each site were selected from among those growing near the lakeshore, as they would be the most likely to have generated subfossil stems (32). Furthermore, we considered that the growth trends of subfossil and living samples can diverge depending on the sampling height of the living trees (43), thus generating biases in the RCS chronologies. To reduce these biases, we sampled all living trees at four meters above the soil surface (i.e., the sampling height that minimized the differences between the growth trends of living and subfossil trees) and we corrected all ring-width series before standardization by subtracting the difference between the local growth curves of the living and subfossil trees (see Figs. S1 and S2).

All of the nonrobust time intervals of the local RCS chronologies were discarded, with the resulting gaps being reconstructed using the analog method (44). Most of these intervals were asynchronous among the lakes and were caused by wildfires that temporarily reduced the input of tree remains into the littoral zone (23, 41). We identified the nonrobust time intervals using the Rbar statistic (mean correlation between pairs of individual series) over 31-y moving windows and discarded all intervals with negative or noncomputable Rbar values due to low replication. At lake L16, the 1978–2011 period was also discarded, as the construction of a road raised the lake water level and disturbed the riparian trees. Over the 1,102 y retained for STREC (see below), we discarded 114, 70, 52, 19, 0, and 31 y for lakes L1, L12, L16, L18, L20, and L22, respectively. The analog method used to reconstruct these gaps is commonly used in paleoclimatology to estimate missing values in proxy-climate matrices. For years where missing values are present, the procedure identifies the most similar years in the calibration period, i.e., all other years without missing values, which are then weighted according to their similarity and used to reconstruct the gaps (44). The measure of similarity is based on the Euclidian distance.

We used the CRU TS3.20 climate dataset (26) to calibrate our reconstruction. The first four years (1901−1904) were excluded because of a poor fit with the tree-ring indexes due to the lack of operating weather stations at the beginning of the last century near the study area. The climate reconstruction method was based on a linear scaling procedure. Each local RCS chronology was rescaled so that its mean and SD matched those of the July−August mean temperature in the calibration period (1905−2011). The final STREC was then obtained by computing the median of the six local reconstructions. This approach attenuated the influence of local nonclimatic disturbances, as the median of the six reconstructions is not sensitive to outliers. It also obtained better cross-calibration verification results (1905−1957 vs. 1958–2011) than a reconstruction based on a partial least squares regression, especially when considering low frequencies (i.e., smoothed datasets; Table S1). Only reconstructed values subsequent to A.D. 910 were retained to limit the analyzed period to the statistically reliable interval (overall replication >53 individual series and Expressed Population Signal >0.85; Fig. 1 C and E). Due to the fast temperature increase in the study area during the 20th century, the reconstructed values in the calibration encompassed a range of 4.3 °C (from 9.6 °C to 14.0 °C), and 86% of the total STREC reconstructed values are inside this range.

We generated realistic time-varying confidence intervals considering two kinds of errors (Fig. 1 A and B). The first one depends on the nonperfect fit between observed and reconstructed values over the verification periods and was computed as the mean root-mean-square error between the two datasets. The second one is due to the variation over time of the strength of the common climate signal for each of the six local chronologies and was computed each year as the interval among the four central reconstructed values from the six chronologies. The consideration of this error produces a time-varying 60% confidence interval.

The volcanic signature in STREC was analyzed using Superimposed Epoch Analysis (SEA) and Bayesian hidden Markov models. We used SEA to test the agreement between the strongest volcanic eruptions of the last millennium and corresponding cooling episodes in STREC. Our SEA was performed in the R environment as described in SI Methods. We used Bayesian hidden Markov models to identify sudden changes in the STREC time series. Such models provide an explicit mechanism to represent transitions between different states and allowed the data to be classified into distinct regimes. The Schwarz criterion is used to identify the number of states and the probability distribution that best fits the data. This approach is briefly described in SI Methods.

All tree-ring and temperature data are included in the Dataset S1 and will also be available from the World Data Center for Paleoclimatology (http://www.ncdc.noaa.gov/paleo/paleo.html).

The authors wish to thank Julia Autin, Yves Bouthillier, Pierre-Paul Dion, Sébastien Dupuis, Benjamin Dy, and Joëlle Marion for field and laboratory assistance; Aurore Catalan for technical assistance; and Joel Guiot for its comments on the paper. This research is a contribution from the ARCHIVES project (http://archives.ete.inrs.ca/) and was financially supported by the Natural Sciences and Engineering Research Council of Canada, Hydro-Quebec, Ouranos, ArcticNet, the EnviroNorth training program, and the Centre for Northern Studies.