Temperature reconstructions from tree-ring densities overestimate volcanic cooling

1 Introduction

Tree-ring proxies are widely available at high northern latitudes and underpin many annually resolved reconstructions of climate over the Common Era [National Research Council, 2006; Jones et al., 2009; Masson-Delmotte et al., 2013]. As they represent the response of dynamic and complex organisms to the ambient environment, some nuance in their interoperation is to be expected. The failure of some tree-ring proxies to record post-1960 warming in the high northern latitudes [Briffa et al., 1998; D'Arrigo et al., 2008] is perhaps the most well-known complexity in interpreting these records solely as temperature proxies.

Stine and Huybers [2014] (hereafter SH14) propose sensitivity of tree-ring densities to light availability as an explanation for the modern divergence of tree-ring series from other markers of temperature. Furthermore, SH14 find a significant correlation between light limitation and the magnitude of the tree-ring density response to volcanic events. SH14 do not, however, explicitly compare the tree-ring response to temperature estimates or control for the spatial pattern of cooling itself.

In contrast, Mann et al. [2012] suggest that in the particularly cold years that generally follow large eruptions, trees in marginal conditions may fail to put down rings. Errors are then introduced into the age chronology, and Mann et al. [2012] therefore argue that tree-ring records underestimate the magnitude of cooling associated with the largest volcanic events. The key counterargument to Mann et al. [2012] is a lack of evidence of missing rings in dendrochronological records [e.g., Anchukaitis et al., 2012; Esper et al., 2013a, 2013b].

Here we further investigate the fidelity of the maximum latewood density response to volcanic cooling in latitudes north of 45°N and the roles played by climatological spatial patterns of the relative importance of light availability and temperature in limiting plant growth [Nemani et al., 2003]. In the Arctic, the intensity of daytime surface illumination is lower than in any other region of the Northern Hemisphere, both because of the low incident Sun angle and because of the prevalence of clouds. Furthermore, the Arctic light environment is generally more diffuse than the midlatitudes and becomes more diffuse in regions with stronger relative light limitation [Stine and Huybers, 2014]. As Arctic trees tend to grow with an open canopy structure, they are expected to be less sensitive than dense canopy trees to differences in the fraction of diffuse versus direct light in driving photosynthesis [Gu et al., 2002]. Thus, we hypothesize that the Arctic and near-Arctic tree-ring density response to volcanism will be amplified as a result of the combined influences of cooling and reduction in light availability, and that the amplification will be largest in regions where baseline growth conditions are most limited by light availability.

Analysis is confined to the 1850–1960 interval during which instrumental records are widely available but before the advent of modern divergence [Briffa et al., 1998]. We consider eruptions that feature nonzero Northern Hemisphere sulfate deposits [Gao et al., 2009a, 2009b], cooling in the instrumental record, and observed decreases in growing season solar radiation in the northern latitudes. The 1912 Alaskan eruption of Katmai, more accurately called Novarupta [Ohvril et al., 2009], and the 1883 Indonesian eruption of Krakatau are the two events that meet these criteria.

2 Data and Methods

Instrumental data are from the CRUTEM3v 5° by 5° gridded compilation (hereafter CRU) [Brohan et al., 2006], which begins in 1850. Proxy data are from the gridded maximum latewood density (MXD) compilation of Briffa et al. [2002a, 2002b], which is based on the Schweingruber compilation of tree-ring chronologies [Schweingruber and Briffa, 1996]. The MXD and CRU data are on the same spatial grid, and we use the 96 MXD grid boxes with centers poleward of 45°N. Following Briffa et al. [2002a, 2002b], we consider April–September averages of the instrumental record and use all CRU time series with at least 10 years of complete monthly data that are poleward of 45°N. These data were used in Tingley and Huybers [2013], and data preprocessing is described therein.

We use a Hierarchical Bayesian Model [Gelman et al., 2003] called BARCAST [Tingley and Huybers, 2010a, 2010b] to infer past temperature anomalies from the MXD and CRU data. The statistical model assumes that temperature anomalies are first-order autoregressive in time, with exponentially decaying spatial covariance. Instrumental observations are modeled as noisy versions of the true latent temperature anomalies, while the proxies are treated as linear in the true temperature anomalies with additive noise. Further details about the statistical model can be found in Tingley and Huybers [2010a, 2010b], Tingley et al. [2012], and Tingley and Huybers [2013]. An initial calibration analysis employs the full data set of Tingley and Huybers [2013], consisting of ice core and varve thickness data, the CRU data, and, to avoid possible biases associated with the late twentieth century divergence, the MXD data up to 1959. We then run BARCAST twice in a reduced mode, drawing parameters from the calibration run and using only the MXD and only the CRU data set, respectively, to produce in each case an ensemble of space-time histories of temperature anomalies. The spatial covariance model included in BARCAST allows for estimation of temperatures and uncertainties in grid boxes without observations, permitting for comparisons between tree-ring and instrumental temperature estimates despite their disparate spatial distributions (Figure 1). Results presented below are based on 4000-member ensembles of predictions, with 90% point-wise uncertainties formed from the 5th and 95th percentiles.

To assess the influence of relative light and temperature limitation on plant growth, we regrid the unitless estimates of Nemani et al. [2003], available on a 0.5° grid, to match the 5° spatial gridding of the CRU and MXD data. These climatological patterns indicate the relative importance of light and temperatures in limiting plant growth, with plant growth expected to be more sensitive to changes in the light environment or temperatures in regions of high radiation or temperature limitation. The spatial domain is then partitioned into two regions, corresponding to high and low radiation limitation, based on the medians of the 5° gridded values, and likewise for temperature limitation (Figure 1). Greenland and the northerly portion of the Canadian Archipelago are excluded from the analysis because of the absence of tree-ring data.

Analysis is confined to volcanic events that occur before the advent of modern divergence in about 1960 [Briffa et al., 1998], in order to reduce the influence of any confounding factors that are unique to the last 60 years. Although the years 1883, 1912, 1925, and 1943 each feature nonzero Northern Hemisphere excess sulfate according to Gao et al. [2009a, 2009b], the instrumental data only indicate cooling in response to the 1912 Novarupta and 1883 Krakatau events (Figure 2 and Figure S1 in the supporting information). Krakatau and Novarupta both featured Volcanic Explosivity Index values of 6 [Simkin and Siebert, 1994] and correspond to the largest Northern Hemisphere sulfate signals of high- and low-latitude origin, respectively. According to the CRU data, both produced domain-wide coolings in excess of 0.3°C sustained over two growing seasons (Figure 2). In contrast, 1925, 1943, and 1944 feature domain-wide warming, while the cooling in 1926 is small compared with those observed following Krakatau and Novarupta (Figure S1). Although volcanic influence may be present in 1925 or 1943 but masked by natural variability, they are excluded as we expect no interpretable signals in tree-ring densities absent identification of a signal in the instrumental record (Figure S1). Furthermore, 1925 and 1943 feature optical depth perturbations that appear at least an order of magnitude smaller than those associated with Krakatau and Novarupta [Sato et al., 1993; Stothers, 1996].

Krakatau and Novarupta both produced strong Boreal growing season reductions in solar radiation intensity at the Earth's surface according to early pyrheliometer observations [Kimball, 1918]. The August 1883 eruption of Krakatau in tropical Indonesia (6°S, 105°E) is notable for producing a long and sustained reduction in solar intensity at the surface of the Earth, reaching the northern midlatitudes in December of 1883 [Stothers, 1996], and persisting though October 1886. The 6 June 1912 eruption of Novarupta in what is now Katmai National Park, Alaska (58°N, 155°W), resulted in a northern latitude surface solar radiation reduction comparable in magnitude but shorter in duration than Krakatau (compare Kimball [1918, Figure 1], Sato et al. [1993, Table 2], Stothers [1996, Table 6], and Ohvril et al. [2009, Figure 3]). These features are consistent with the expectation of a longer-lived signal from a tropical eruption as compared with a summer season high-latitude eruption [Kravitz and Robock, 2011]. Although pyrheliometer observations indicate reduced light intensity in response to the October 1902 eruption of Santa Maria in Guatemala, the radiation signal diminishes rapidly over the first few months of 1903, before the start of the boreal growing season [Kimball, 1918; Stothers, 1996]. The corresponding ice core sulfate signal is small and confined to the Southern Hemipshere [Gao et al., 2009a, 2009b], and the instrumental data we use do not show cooling.

The 1912 Novarupta event is distinct from the primarily tropical or unknown-location eruptions considered by SH14. The single near-Arctic eruption included in SH14 was the extended eruption of Laki in 1783. In contrast, the Novarupta event was explosive and occurred early in the Boreal growing season when resulting variations in temperature and light availability are expected to have a maximal effect on tree growth. The strongest temperature response is seen for the year of eruption (Figure 2a).

The larger, tropical eruption of Krakatau produced a delayed temperature response at high northern latitudes, presumably due to slow meridional transport of stratospheric scatterers to the Arctic [Dyer and Hicks, 1965; Robock, 2000], with the largest growing season temperature response in 1884 (Figure 2e). Although Krakatau was included in the analysis of SH14, the tree-ring response was not compared directly against instrumental temperature records.

To isolate the volcanic cooling signal, we consider anomalies with respect to a baseline mean calculated over 3 years prior to the initialization of volcanic cooling, corresponding to a 1909–1911 for Novarupta and 1881–1883 for Krakatau. For predictions using either the MXD or the instrumental data set, anomalies are calculated separately for each location and ensemble member prior to averaging over various regions and calculating percentiles.

3 Results

The cooling associated with Novarupta, as inferred from the CRU data and averaged over the entire spatial domain, lasted for 2 years (Figure 2a). For 1912, there is a clear bias between the MXD and CRU inferences, with the MXD indicating stronger cooling. The bias between the tree-ring density and instrumental response is largest in magnitude when averaging over the high radiation limitation region (hereafter, “High R”), and smallest for the low radiation limitation (“Low R”) region (Figures 2b and 2c). For all three averaging regions, the bias is highly significant [Probability (Prob.) MXD colder than CRU ≥0.99; Figures 3a–3c and Table S1]. More notably, there is a significantly greater bias in the High R, as compared with Low R, region in all 4000 ensemble members (Figures 2d and 3d and Table S1).

Results are qualitatively similar, but somewhat weaker, for the response to Krakatau, which features a more complex temporal evolution (Figures 2e–2h). For domain-wide averages, we detect no bias between the MXD and CRU for 1884 or 1885, and strong, significant bias for 1886. The domain-wide response averaged over 1884–1886 features significant bias [Prob. MXD colder than CRU =0.99]. For consistency with the treatment of the Novarupta eruption, we reference the Krakatau response to 1884, corresponding to the strongest cooling signal. Doing so also ensures that results are not influenced by the small optical depth perturbations produced by the June 1886 eruption of Tarawera in New Zealand [Stothers, 1996]. For 1884, there is no significant bias for the Low R region, significant bias for the High R region [Prob. MXD colder than CRU =0.97], and significantly greater magnitude bias in the High R as compared with the Low R region [Prob. greater bias in High R =0.99; Figure 3h and Table S2]; these regional differences are qualitatively unchanged when averaging the Krakatau cooling response over 2 or 3 years.

Results support the SH14 hypothesis that reduction in light availability accentuates the MXD response to volcanic events in regions that are most light limited. For the response to both Novarupta and Krakatau, tree-ring density records show a larger magnitude response than the instrumental records, with significantly greater bias in the High R as compared with the Low R region. These results indicate that the larger temperature response inferred from the MXD for the High R region as compared with the Low R region is the result of greater bias in regions where the MXD data are more sensitive to light availability, rather than a feature of the spatial pattern of the cooling itself. Indeed, the CRU data show that the High R region is warmer than the Low R region following Novarupta, and that the High and Low R regions featured similar cooling following Krakatau (Figures 3d and 3h). In contrast to the effects of radiation limitation, differences in MXD-CRU biases are not significant when comparing averages over regions of high and low temperature limitation (Tables S1 and S2).

Despite these biases, there is notable agreement between the spatial patterns of inferred cooling from the MXD and CRU data sets (Figure 4). Indeed, it is the general fidelity of the MXD inferences in reproducing patterns seen in the instrumental record that permits for our exploration of second-order features associated with differences. For the 1912 response, the correlation between the mean CRU and MXD response fields is r = 0.56, and the correlation is positive for each of the 4000 ensemble members. Results for the 1884 response are similar, with a correlation of r = 0.55 between the mean response fields [Prob.(r > 0) > 0.99].

The MXD-inferred response patterns are similar for both Novarupta and Krakatau (Figures 4a and 4d), with cooling throughout much of Eurasia and Eastern Canada, spotty warming near Scandinavia and Europe, and more cohesive warming from northwestern Canada to far eastern Siberia (Figures 4a and 4d). The temperature responses are generally largest in regions that are remote from the moderating influences of the ocean, and the observed features are consistent with previous reports of the summer temperature anomaly associated with large volcanic eruptions [Robock and Mao, 1995]. The correlation between the two MXD mean response patterns is r = 0.60 and is highly significant [r > 0 for all ensemble members]. In contrast, the correlation between the two CRU-inferred mean response patterns is positive (r = 0.21) but not significant [Prob.(r > 0) = 0.82; Figures 4b and 4e]. These results are consistent with similarities between the two MXD-inferred response patterns not exclusively owing to a common temperature response, and with regional patterns of light limitation remaining similar between eruptions.

Nontemperature influences on the two MXD-inferred spatial response patterns can be explored in greater detail using the bias fields (Figures 4c and 4f). Both feature areas of significant overestimation of the strength of cooling in Northern Eurasia—an area that SH14 identify as featuring the strongest Eurasian volcanic reduction in MXD and the most pronounced late twentieth century divergence. In addition, the bias field for Krakatau features significant overestimation of cooling in Northern Canada and a warm bias in far eastern Eurasia. The bias field for the response to Novarupta suggests similar features, though they are not independently significant. Despite the low correlations between the instrumental-inferred coolings, the correlation of r = 0.49 between the two mean bias fields is significant [Prob.(r > 0) = 0.99], further supporting a nontemperature control on the pattern of MXD response. Results are qualitatively unchanged when averaging the Krakatau response over 1884 and 1885.

The warm MXD bias for 1884 near the Sea of Okhotsk occurs in a region where both CRU and MXD data are sparse (Figures 4d and 4e). Although the uncertainty in mapping temperature estimates to these regions is accounted for in BARCAST, it is important to note that the instrumental record at the center of the warm bias reports significant cooling, while the closest, but more distantly located, MXD records show a warming that is not statistically significant (Figures 4d and 4e). These features combine to produce an estimated warm bias for 1884 in a region where, albeit not deemed significant, there is also a warm bias in 1912. Note also that there are numerous data-rich regions that display significant biases of both signs, including Northern Eurasia in both years (cold bias), Europe in 1912 (cold bias), and Northern Scandinavia in 1912 (warm bias).

Two separate regression models confirm that radiation limitation, but not temperature limitation, is a significant predictor of the pattern of bias between the volcanic responses inferred from MXD and CRU records. Following the division of space into regions of high and low radiation and temperature limitation (Figure 3), we first define indicator variables, such that I_R takes on the value zero for locations where radiation limitation is below the median value and one for locations where it is above the median, and likewise for I_T. The regression model then takes the form,

(1)

where Δ is the MXD–CRU bias (in °C), μ is the spatial average bias, ε is an error or discrepancy term, and replication is over space. β_R gives the expected increase in bias in the region of high, as compared with low, radiation limitation, and likewise for β_T for temperature limitation. In the second regression formulation, the indicator variables are replaced with the 5° averages of the numerical values from Nemani et al. [2003]. Fitting each model separately for each of the 4000 realizations of the cooling bias in response to either Krakatau or Novarupta, the coefficients β_R are in all cases significantly below zero [Prob.(β_R < 0) ≥ 0.95], while 90% uncertainty intervals for β_T cover zero (Tables S3 and S4). These regressions indicate that light availability, but not temperature, significantly controls the spatial pattern of bias.

4 Discussion and Conclusions

The two eruptions considered here are of different character, with the 1912 Novarupta injecting scatterers directly into the atmosphere at high northern latitudes in the midst of the growing season. There is an immediate impact on climate in regions poleward of 45°N, with significant cooling indicated by instrumental and MXD records in 1912 (Figure 2, Table S1). The instrumental records, however, indicate persistent cooling through the 1913 season, whereas the MXD records show a concentrated response in 1912. Insomuch as MXD also records variations in light, it is possible that the direct radiative effects decayed within a year but that thermal inertia carried the cooling signal through 1913. Indeed, pyrheliometer records indicate that the reductions in light availability following the Novarupta eruption were confined to 1912 [Kimball, 1918; Stothers, 1996]. The response to Novarupta shows substantial bias between MXD and instrumental estimates in terms of magnitude and spatial distribution (Figure 4). The domain-wide MXD-inferred response is significantly larger in magnitude than indicated by the instruments, and the bias is significantly larger in regions of high as compared to low light limitation (Prob. >0.99; Figures 2 and 3d; Table S1).

The high-latitude growing season temperature anomaly in response to the larger, tropical eruption of Krakatau in 1883 does not appear until 1884. The cooling signal then persists through 1885 in the instrumental records and through 1886 in the MXD (Figure 2e). The bias between MXD and instrumental estimates is small in 1884 but becomes more evident in 1885 and 1886; the domain-wide bias averaged over these 3 years is significant [Prob. MXD colder than CRU =0.99]. Similar to Novarupta, there is significantly greater bias between the MXD and instruments in high as compared with low radiation limitation regions (Prob. =0.99; Figures 2 and 3h; Table S2). Although atmospheric circulation perturbations in response to volcanic eruptions are complex, depending on both the season and location of an eruption [Kravitz and Robock, 2011], we speculate that the temporal evolution of the bias in response to Krakatau (Figure 2) may be a result of both dynamical effects resulting from initial cooling in the tropics and a reduced poleward transport of heat, and, later, the direct influence of scatterers once they arrive in the Arctic and sub-Arctic stratosphere. Relative contributions to Arctic cooling of the dynamic response to equatorial cooling and the slow stratospheric transport of reflective sulfate could usefully be explored in the context of other volcanic events and model simulations.

The patterns of bias between the MXD and CRU-inferred responses to the volcanic eruptions point to the importance of spatial patterns of light limitation in modulating the MXD response. Although the instrumentally inferred cooling fields in response to Novarupta and Krakatau show no significant correlation, the two spatial fields of bias are significantly correlated (Prob.(r > 0) = 0.99; Figure 4), and we find that, for both eruptions, the sensitivity of plant growth to changes in the light environment is a significant predictor of the spatial bias pattern [Prob.(β_R < 0) ≥ 0.95; Tables S3 and S4].

Debate exists as to the net effect of volcanic scatterers in the atmosphere on global net primary productivity [Roderick et al., 2001; Angert et al., 2004], as scattering not only decreases the total sunlight reaching the surface but also converts some direct insolation into more photosynthetically efficient diffuse radiation [Roderick et al., 2001; Robock, 2005]. The larger magnitude of volcanic response we find for the MXD as compared with CRU temperatures is indicative of an overall decrease in productivity, and is consistent with the relatively diffuse Arctic light environment and the generally open canopy structure of Arctic trees making an increase in the diffuse fraction insufficient to compensate for an overall reduction in light availability. In regions of low light limitation, the bias between the MXD and CRU response is small or absent (Figures 3c and 3g), consistent with scattering leading to a greater increase in the diffuse fraction when more direct light is initially present. See SH14 for a more complete discussion.

Results are consistent with and extend upon the findings of SH14 in several important ways. By transforming the MXD response into temperature anomalies and mapping the spatial pattern of response, we are able to directly compare the patterns of MXD and instrumental volcanic responses. Furthermore, the transformation and mapping algorithm provides comprehensive uncertainty estimates whose coverage intervals have been verified in synthetic tests [Tingley and Huybers, 2010a, 2010b], permitting for direct comparison of instrumental and tree-ring density data. No such direct comparison was possible in SH14 due to the focus on preinstrumental volcanic events and the MXD being neither transformed nor mapped. In contrast, our results provide little evidence for or against the missing ring hypothesis of Mann et al. [2012]. The generally cold bias we identify does not necessarily contradict Mann et al. [2012] because the cooling bias may overwhelm the influence of any missing rings.

Our results complicate the interpretation of high-latitude tree-ring density records in terms of a pure temperature response. In particular, strong coolings in the Polar Ural regions seen following many volcanic eruptions [Briffa et al., 2002a, 2002b] may not be indicative of the true spatial pattern of cooling. Mechanistic models of tree-ring widths include moisture, temperature, and light availability as climatic inputs [Evans et al., 2006; Tolwinski-Ward et al., 2011]. A similarly comprehensive approach to modeling tree-ring densities may improve climate inferences and may, in conjunction with other proxies that feature responses to multiple environmental variables [Tolwinski-Ward et al., 2014], allow for simultaneous inferences on both temperature and light.