Paleoclimatic evidence is necessary to place the current warming and drying of the western Mediterranean basin in a long-term perspective of natural climate variability. Annually resolved and absolutely dated temperature proxies south of the European Alps that extend back into medieval times are, however, mainly limited to measurements of maximum latewood density (MXD) from high-elevation conifers. Here, the authors present the world’s best replicated MXD site chronology of 414 living and relict Pinus uncinata trees found >2200 m above mean sea level (MSL) in the Spanish central Pyrenees. This composite record correlates significantly (p ≤ 0.01) with May–June and August–September mean temperatures over most of the Iberian Peninsula and northern Africa (r = 0.72; 1950–2014). Spanning the period 1186–2014 of the Common Era (CE), the new reconstruction reveals overall warmer conditions around 1200 and 1400, and again after around 1850. The coldest reconstructed summer in 1258 (−4.4°C compared to 1961–90) followed the largest known volcanic eruption of the CE. The twentieth century is characterized by pronounced summer cooling in the 1970s, subsequently rising temperatures until 2003, and a slowdown of warming afterward. Little agreement is found with climate model simulations that consistently overestimate recent summer warming and underestimate preindustrial temperature changes. Interannual–multidecadal covariability with regional hydroclimate includes summer pluvials after large volcanic eruptions. This study demonstrates the relevance of updating MXD-based temperature reconstructions, not only back in time but also toward the present, and emphasizes the importance of comparing temperature and hydroclimatic proxies, as well as model simulations for understanding regional climate dynamics.
The Mediterranean basin, with its well-described west–east climate dipole (Xoplaki et al. 2004; Roberts et al. 2012; Seim et al. 2015), has been termed a “climate change hotspot” (Giorgi 2006), where climate warming has amplified drying since the 1970s (Kelley et al. 2012; Philandras et al. 2011; Hoerling et al. 2012; Cook et al. 2015, 2016). Drought-prone forest ecosystems already show signs of suffering from changes in the hydrological cycle (Büntgen et al. 2012a, 2013, 2015a,b; Galván et al. 2014, 2015). Further, climate projections almost uniformly point to a continuation of western Mediterranean drying (Seager et al. 2014; Fischer and Knutti 2015; Kelley et al. 2015).
To contextualize the ongoing and predicted changes in drought frequency and severity within a multicentennial-to-millenial perspective, paleoclimatic evidence is needed (Ljungqvist et al. 2016). A wide range of proxy archives with subdecadal-to-centennial resolution from the western Mediterranean basin allow comparison of the recent hydroclimatic conditions with those that may have occurred during the relatively warm medieval period [see Luterbacher et al. (2012) for a review]. Although several lake and marine sediment cores, and speleothem records, indicate relative drying during periods of medieval times (Martínez-Cortizas et al. 1999; Martín-Puertas et al. 2010; Morellón et al. 2011; Nieto-Moreno et al. 2011; Pérez-Sanz et al. 2014; Wassenburg et al. 2013; Servera Vives et al. 2014; Sánchez-López et al. 2016), some proxies from the Iberian Peninsula suggest wetter conditions approximately 1000–1200 yr ago (Benito et al. 2003; Corella et al. 2016; Barreiro-Lostres et al. 2014). The most detailed depiction of Mediterranean drought variability over the last 900 yr is presented in the Old World Drought Atlas (OWDA; Cook et al. 2015, 2016), a tree-ring-based field reconstruction of warm-season self-calibrated Palmer drought severity indices (scPDSI; van der Schrier et al. 2013). In contrast to latitudes ≳50°N, there are substantially fewer well resolved, high-quality temperature reconstructions from lower latitudes (Esper et al. 2016). Indeed, there are no temperature sensitive tree-ring width chronologies for the region south of the Alpine arc (Büntgen et al. 2012b), although inverse relationships have been reported for Iberia’s southeast (Dorado-Liñán et al. 2015).
The majority of annually resolved and absolutely dated, near-millennium-long temperature reconstructions from the Mediterranean region are based on maximum latewood density (MXD) measurements from high-elevation pine trees in the Spanish central Pyrenees (Büntgen et al. 2008, 2010; Dorado-Liñán et al. 2012; Galván et al. 2015). However, their limited sample size before ~1600 CE and modern terminations in 2005 hamper straightforward comparison with recent temperature trends and extremes to those that may have occurred, for instance, during medieval times. Despite the pioneering attempt by Dorado-Liñán et al. (2012), a combined assessment for southwest Europe of refined climate reconstructions, state-of-the-art climate model simulations, and precise climate forcing records is missing. A thorough analysis of the effects of volcanic eruptions on western Mediterranean temperature and hydroclimatic variability during industrial and preindustrial times is also needed (Anchukaitis et al. 2012; Esper et al. 2013a,b).
Here, we address the aforementioned issues by introducing the most recent and world’s best replicated MXD site chronology. This record is an update of 414 MXD measurement series from living and dead high-elevation mountain pines in the Spanish central Pyrenees (Büntgen et al. 2008). The new composite chronology robustly covers the period 1186–2014 CE, and reflects variation in May–June and August–September temperatures over the Iberian Peninsula and northern Africa.
2. Data and methods
Core and disc samples were collected from living, dry-dead, and subfossil mountain pines (Pinus uncinata Ram.) of diverse ages, within the treeline ecotone surrounding Lake Gerber in the Spanish central Pyrenees (Büntgen et al. 2008). This remote, undisturbed site is located at the northeastern limit of the Aigüestortes i Estany de Sant Maurici National Park, near Port de la Bonaigua (42°374′5.15″N, 00°59′37.18″E). The study area is difficult to access, with its lowest margin defined by Lake Gerber at 2166 m above mean sea level (MSL) and the highest trees growing ~2400 m MSL.
Mountain pine is a shade-intolerant conifer species that dominates the subalpine zone of the Pyrenees >1600 m MSL (Vigo and Ninot 1987). Radial growth of this species is enhanced by warm conditions during November in the year prior to tree-ring formation, and during May of the growth year (Tardif et al. 2003). This second response agrees with phenological studies of radial growth that found the highest growth rates to occur in May, and latewood formation takes place from July to September (Camarero et al. 1998).
All MXD measurements were made on a Walesch 2003 X-ray densitometer at the Swiss Federal Research Institute for Forest, Snow, and Landscape [Eidgenössische Forschungsanstalt für Wald, Schnee und Landschaft (WSL)] in Birmensdorf (Eschbach et al. 1995). To avoid biases from stem-internal changes in fiber direction, each sample was cut in 2–4-cm-long sections that were 1.23 mm thick, perpendicular to the tracheid’s longitudinal axis. Wood density—herein measured as the ratio between tracheid wall thickness and its transverse lumen area—was continuously recorded along the entire sample at a resolution of 0.01 mm. Using a five-step calibration wedge, the X-ray image brightness variations were converted into grams per centimeter cubed [correction factor (CF) = 0.878; Lenz et al. 1976]. Of the 414 samples, 377 were measured using an earlier version of the Dendro 2003 software (Büntgen et al. 2008), and 37 of the most recent update samples were measured using the latest version (WalDendro V1.10).
Different detrending methods and chronology development techniques were applied to evaluate the removal of nonclimatic, age-related trends inherent to the raw MXD measurements. Individual series detrending via cubic smoothing splines of 50% frequency-response cutoff at 40 and 200 yr (Cook and Peters 1981), and negative exponential functions were utilized for a first comparison of the high-frequency variability in the original (261 series) and updated (153 series) MXD chronologies (not shown). After consideration of pith-offset estimates for all core and disc samples that lacked pith, 14 slightly different regional curve standardization (RCS) chronologies (Esper et al. 2003) were produced so as to preserve lower-frequency information in the complete dataset of 414 series (Cook et al. 1995). Dimensionless indices were calculated as residuals after power transformation (PT; Cook and Peters 1997). The resulting index chronologies were truncated at a minimum replication of eight series, and their variance stabilized with a 31-yr moving window approach (Osborn et al. 1997).
Eight different temperature records were compiled to evaluate their suitability as calibration targets. Monthly-resolved temperature means from the closest grid points in the Berkeley (Rohde et al. 2013), HadCRUT4.4 SST 2-m temperature (T2m), HadCRUT4.4 SST filled-in T2m, CRU Time Series v3.23 (CRUTS3.23; Harris et al. 2014), and E-OBS v12 gridded instrumental temperature products (Haylock et al. 2008), as well as from adjusted Global Historical Climatology Network, version 3 (GHCNv3), station measurements (Lawrimore et al. 2011) at Zaragoza Airport (41.67°N, 1.02°E), Perpignan (42.73°N, 2.87°E), and Toulouse-Blag (43.63°N, 1.37°E), were transformed into anomalies with respect to the 1961–1990 mean climatology. A simple correlation matrix over the 1950–2015 common period, and moving 31-yr standard deviations were used to assess spatial measurement coherency and temporal variance changes. Because of a decreasing number of instrumental station records, together with an increasing amount of uncertainty associated with these measurements before the 1940s (Frank et al. 2007; Jones 2016), and particularly across the greater Pyrenees region (Büntgen et al. 2008, and references therein), the growth-climate response analyses and proxy-target calibration trials were restricted to the post-1950 period for which European-wide, high-resolution reanalysis data are available (E-OBS v12; Haylock et al. 2008; van den Besselaar et al. 2011).
The climate sensitivity of 14 near-similar RCS chronologies was evaluated against monthly and seasonally resolved minimum, mean, and maximum temperatures, as well as precipitation totals and sea level pressure indices from a spatial subset of the E-OBS network (42.5°–43.0°N, 0.5°–1.0°E). A simple scaling approach was used for calibration (Esper et al. 2005). Both the mean and variance of the RCS MXD chronology that produced the highest monthly correlations (hereafter PMR chronology) were adjusted to mean May–June and August–September (MJ&AS) temperatures over the 1950–2014 period. This procedure, the simplest among calibration techniques, is less susceptible to variance underestimation than direct regression (Bürger et al. 2006; Christiansen 2011), and it is not as prone to variance overestimation by indirect regression (Moberg 2012; Christiansen and Ljungqvist 2017). Uncertainty ranges associated with the new temperature reconstruction not only consider the amount of unexplained variance in the 1950–2014 calibration period [±1 root-mean-square error (RMSE)] but also reflect possible biases emerging from the MXD data (bootstrap confidence intervals), as well as the detrending (RCS) and index calculation (ratios or residuals after PT) techniques applied [see Esper et al. (2007) for error estimates].
All three previously published warm-season temperature reconstructions from the greater Pyrenees region (Büntgen et al. 2008; Dorado-Liñán et al. 2012; Esper et al. 2015), a recent MXD-based reconstruction of Northern Hemisphere summer temperature variability (Schneider et al. 2015), and a regional subset from the OWDA (Cook et al. 2015), centered over Lake Gerber (42.25°N, 0.75°E), were used for comparison with the new PMR chronology. A total of four simulated MJ&AS temperature means, CCSM4 (Gent et al. 2011; Landrum et al. 2013), GISS-E2-R (with global primary production of 121 PgC yr−1 at preindustrial CO2 concentrations; Schmidt et al. 2006), IPSL-CM5A-LR (Dufresne et al. 2013), and MPI-ESM-P (Bothe et al. 2013; Giorgetta et al. 2013) from phase 5 of the Coupled Model Intercomparison Project (CMIP5; Taylor et al. 2012) and extracted for the Pyrenees (centered over Lake Gerber), were also used for comparison. The approximate horizontal resolutions of the four CMIP5 models are 1.00° latitude (lat) × 1.15° longitude (lon) for CCSM4, 2.00° lat × 2.50° lon for GISS-E2-R, 1.90° lat × 3.75° lon for IPSL-CM5A-LR, and 1.90° lat × 1.90° lon for MPI-ESM-P. Additional simulations outside CMIP5 but with similar models in terms of the MPI-ESM-P “past1000” ensemble and the CCSM4 Last Millennium Ensemble (Otto-Bliesner et al. 2016) were considered for the assessment of specific temperature responses to the large tropical volcanic eruptions of Samalas (part of Rinjani volcano, Indonesia) in 1257 CE (Lavigne et al. 2013; Vidal et al. 2015) and 1809/1815 (unknown/Tambora, Indonesia; Stothers 1984; Oppenheimer 2003) (see the discussion section for details). All simulations include orbital, solar, volcanic, and greenhouse gas forcing (Schmidt et al. 2011, 2012). Most simulations also implement effects of transient land-use/land-cover changes.
Different estimates of total solar irradiance changes (Steinhilber et al. 2009, 2012), as well as a revised bipolar record of volcanic eruptions (Sigl et al. 2015), were utilized as a guideline for temporal dynamics in the external drivers of natural climate variability over the last millennium. Superposed Epoch Analysis (SEA; Chree 1913) was performed to quantify the behavior of the new MJ&AS temperature reconstruction and the regional OWDA subset after volcanic eruptions. A total of 20 known (11) or suspected (9) eruption dates for identified (7), putative (5), or unknown (8) source volcanoes have been selected for comparison with the new PMR chronology (see the results and discussion section for further details). This set of large tropical or Northern Hemisphere extratropical eruptions with estimated volcanic explosivity index (VEI) values thought to likely exceed 5 is based on the approximate coincidence of peaks in volcanic forcing time series (Sigl et al. 2015) and evident depressions in the new PMR chronology. Temporal mismatch of up to three years during the time of interest can be accommodated by small dating uncertainties in the ice core sulfate record (Esper et al. 2013b; M. Sigl 2016, personal communication). Another reason for slight temporal adjustments of the SEA may originate from an inconsistent treatment of eruption years and climate forcing years (e.g., the Tambora eruption in April 1815 is associated with the “year without a summer” in 1816; Stothers 1984; Luterbacher and Pfister 2015; Raible et al. 2016) (see also discussion). Split-period SEA was applied for 17 preindustrial and 3 industrial era eruptions, before and after 1850 CE, respectively. Because of low replication, caution is advised when interpreting the post-1850 response pattern. Significance levels of all correlation coefficients (p ≤ 0.01) were corrected for a reduction in the degrees of freedom due to positive first-order autocorrelation structures.
Substantial agreement exists between the original (261 series) and newly updated (153 series MXD chronologies from the Spanish central Pyrenees (Figs. 1a–c), which are entirely independent and have no data in common. The mean correlation coefficient between the two records after power transformation, RCS detrending, and variance stabilization is 0.86 for the period 1258–2005 CE, when both datasets contain at least eight individual measurement series. Moving 31-yr correlations between these records reveal a high level of coherency over most of the last seven centuries (r > 0.9). However, there is slightly lower agreement before ~1350, when sample size decreases, and in the middle of the eighteenth century, when material from living and relict trees overlap. The temporal distribution of the newly sampled 153 MXD series considerably improves replication of the original dataset before ~1600 and extends its modern termination to 2014 (Fig. 1d).
The 414 samples from living and relict Pinus uncinata trees >2200 m MSL from the Spanish central Pyrenees represent the world’s best replicated MXD site chronology (Fig. 2a). The average (relative) MXD per sample is 0.64 g cm−3 with minimum and maximum values of 0.44 and 0.84 g cm−3 respectively. MXD values (as herein defined) generally decrease from juvenile to mature growth (Fig. 2a, inset). Start and end dates of all 414 individual series are distributed between 924 and 2014 CE, and the mean segment length (MSL) ranges from 11 to 732 yr, with an average of 223 yr. The age structure of the new compilation is extremely homogeneous between the late thirteenth and the late twentieth century, whereas overall older and younger trees dominate the periods before and afterward (Fig. 2b). The combination of water logs, dry-dead snags, and living trees permits usage of the RCS methodology, producing a robust chronology between 1186 and 2014 (Fig. 2c). Despite possessing a seemingly large degree of interannual–multicentennial variability, the PMR chronology is characterized by an exceptionally low first-order autocorrelation (r = 0.1).
For northern Spain and southern France, homogenized and adjusted instrumental station measurements, as well as high-resolution gridded products, are available from 1950 onward (Fig. 3a). Gridded data are essentially consistent over the past six decades. The average correlations >0.84 between MJ&AS temperature variations are significant (p ≤ 0.01) (Fig. 3a, inset). The spatially interpolated data, reaching back to 1750, show an initial summer warming in the 1950s followed by two decades of relatively cool conditions, and increasing temperatures from the early 1980s to ~2000, with no warming afterward. The coldest and warmest summers since 1750 occurred in 1972 and 2003, respectively. Moving 31-yr standard deviations of the longer records reveal an obvious collapse in variance beginning in the first half of the twentieth century (Fig. 3b) and extending to their earliest values.
The strongest relationships between slightly different versions of the PMR chronology and monthly resolved E-OBS temperatures are seen for May (r ~ 0.61), June (r ~ 0.36), August (r ~ 0.54), and September (r ~ 0.35) (Fig. 4a). Nonsignificant correlations are found with temperatures in the year prior to ring formation, which is in agreement with the overall very low first-order autocorrelation structure of the PMR-chronology types. January–April, July, and October–December temperatures in the year of tree growth also have no relationships with MXD, which is confirmed by nonsignificant correlations between the monthly temperature means (from −0.08 to 0.39), as well as the nonsignificant lag-1 autocorrelation structures of the individual time series (from −0.04 to 0.23) (Table 1, top). Most interesting is the unimportance of July (see the discussion section for details), when monthly temperature means (14.4°C) are highest and precipitation totals (54.8 mm) are lowest (Fig. 4b) (based on the 1950–2015 measurement period). The annual mean temperature and total precipitation is 5.86°C and 1158 mm, respectively. Slightly different seasonal temperature means between May and September (excluding July) produce correlations ~0.60 with MXD (1950–2014).
Interannual–multidecadal variability of the scaled PMR chronology is in agreement with the measured MJ&AS temperature means of the nearest E-OBS grid box (Fig. 5a). Correlations between the PMR chronology and other gridded temperature products not only demonstrate the high quality of the E-OBS data (Table 1, bottom) but also indicate increasing disharmony between the PMR chronology and the different meteorological measurements back in time. Correlation between the PMR chronology and E-OBS is 0.72 (1950–2014), with higher and lower coefficients during the first and second half of proxy–target overlap (r = 0.88 and 0.53). Both time series reveal positive temperature anomalies in the 1950s, followed by almost two decades of relatively cool summers with negative extremes in 1972 and 1984 (Fig. 5a). Although temperature fluctuations since ~1990 occurred at a relatively high level, the overall warming trend until 2003 was only moderate. Most interesting is the fact that both the measured and reconstructed temperatures exhibit a plateau since the 1990s, with a constant decline after 2003 (see the discussion section for details). The first-order autocorrelation structure of the proxy and target time series is low (r = 0.10–0.12), and agreement in the timing of annual extremes is high. Split and full period calibration–verification statistics confirm a sufficient level of reconstruction skill (Table 2).
The PMR chronology explains a significant (p ≤ 0.01) portion of MJ&AS temperature variability over most of the western Mediterranean basin (Fig. 5b). The highest spatial field correlations cover northeastern Spain, southern France, and northern Algeria (r > 0.6), whereas no relationships are found north of 50°N and east of 15°E. This proxy-based pattern of explained summer temperature variability coincides fairly well with the results from the gridded data themselves (Fig. 5b, inset), that is, when correlating the MJ&AS temperature means from the single E-OBS grid box over the study site against all E-OBS grid boxes over Europe. Although the new proxy–target fit is an improvement over the initially obtained relationships from the 2008 work (Büntgen et al. 2008), there is a slightly decreasing linear trend in the annual residuals between the reconstructed and measured temperatures (Fig. 5c). This feature, consistent with a slightly increasing uncertainty range since around the 1990s, indicates some sort of recent disharmony in the proxy–target association. The density distribution histograms of the proxy and target data, however, show comparable shapes for the reconstructed and measured values (Fig. 5d).
The anomalous cold conditions of 1972 are well preserved in both the target and proxy data, which confirms the exceptional skill of the MXD parameter to capture realistically negative temperature extremes. However, the positive outlier of 2003, well reflected in the instrumental measurements, is less pronounced in the PMR chronology. Global-scale field correlations of the new MJ&AS temperature reconstruction reveal rapidly declining associations with increasing distance from the proxy site (Fig. 6a). While a large amount of summer temperature variability is explained over the western Mediterranean basin, no relevant teleconnection patterns are found. This regional scope is corroborated by spatial correlations of the MJ&AS temperature means from a single E-OBS grid box over the study site (Fig. 6b). The spatial correlation field suggests some positive associations with parts of central Asia and the North Atlantic, resembling the resulting circumpolar wave train teleconnection pattern of surface air temperature that modulates European and Eurasian summer climate.
The new MXD-based MJ&AS temperature reconstruction from the Spanish central Pyrenees covers the period 1186–2014 CE (Fig. 7a). Although the record reflects year-to-year but also longer-term fluctuations, its first-order autocorrelation is only 0.17. Clusters of particularly warm summers occurred at ~1200, ~1400, and ~1900 (Table 3), whereas negative extremes are more equally distributed throughout time. The coldest MJ&AS anomaly occurred in 1258 (−4.4° ± 1.5°C) after one of the largest volcanic eruptions of the last two millennia (Table 4). The overall coolest conditions are reconstructed for parts of the thirteenth century, and slightly below average temperatures continued from ~1600 until the mid-nineteenth century. The twentieth century, characterized by a distinct summer cooling in the 1970s, falls within the envelope of preindustrial variability, and the reconstructed MJ&AS temperatures, in line with instrumental measurements (Fig. 5a), started to decline again in 2003. The new reconstruction, while resembling the high- to low-frequency information of previous MXD-based records (Büntgen et al. 2008; Dorado-Liñán et al. 2012) (Fig. 7a, inset), suggests a century of relatively warmer summers centered on 1400 CE. This warm phase coincides with a decadally resolved, June–August temperature reconstruction based on detrended carbon isotopic (δ13C) values from the same site (Esper et al. 2015), which shows the largest departure from the MXD-based approaches in the late-nineteenth century (Fig. 7a).
The new MJ&AS temperature reconstruction from the Pyrenees challenges regional output from state-of-the-art climate model simulations (Fig. 7b). In comparison with the reconstructed preindustrial low-frequency summer temperature variability, all four model simulations of MJ&AS temperature reflect a somewhat dampened amplitude. A limited agreement between proxy-based reconstructions and model-based simulations on interannual time scales was expected, unless volcanic forcing dominates. This dissimilarity is reflected in overall nonsignificant correlations between the 20-yr low-pass filtered records (Fig. 7b, inset). In contrast to the rather weak proxy–model agreement, there is a remarkable synchrony between the Pyrenees data and a new MXD-based Northern Hemisphere summer temperature reconstruction (Schneider et al. 2015) (Fig. 7c). Correlations between the various reconstructed Pyrenees (Büntgen et al. 2008; Dorado-Liñán et al. 2012; this study) and NH (Schneider et al. 2015) temperature means range from 0.45 to 0.48 and from 0.59 to 0.70 when using the original and 20-yr low-pass filtered time series, respectively (Fig. 7c, inset). It should be noted in this regard, that the Schneider et al. (2015) NH record includes MXD measurements from Büntgen et al. (2008).
In contrast to the reconstructed low-frequency variation in western Mediterranean summer temperatures, reconstructions of hydroclimate are restricted to much shorter time scales (Figs. 8a,b). Most of the decadal-to-multidecadal pluvials in the regional OWDA subset coincide with relatively cold episodes. The observed recent drying since the 1970s falls well within the envelope of natural hydroclimatic variability since medieval times. The positive and negative scPDSI extremes are rather uniformly distributed throughout time (Table 3), whereas the five wettest decades correspond to periods of increased volcanic activity. The strong relationship between higher summer temperatures and increased drought is reflected in negative correlations between the original and smoothed time series over the 1186–2012 period of overlap (Fig. 8b, inset).
The overall correlation, over the full length of both the new temperature reconstructions and the regional OWDA subset, is r = −0.27 at annual resolution. This negative relationship increases to −0.44 after 20-yr low-pass filtering, remains at this level at lower-frequency domains, and also holds throughout time. An interesting, but not distracting, feature is the relatively high first-order autocorrelation of the more “reddish” ring-width-based OWDA data (r = 0.36) compared to the “whiter” MXD-based temperature reconstruction (see the discussion section for details).
Based on a set of 20 large volcanic eruptions that occurred between 1209 and 1992 (Table 4), the SEA shows negative and positive responses for the temperature and hydroclimatic reconstructions (Figs. 8c,d), respectively. Exceptional summer cooling together with above average precipitation totals are typical features after strong volcanic eruptions (see the discussion section for details). The reconstructed average MJ&AS temperature depression (with respect to the 1961–90 reference period) following 20 large volcanic eruptions is −1.96°C (±0.87°C). The abrupt cooling is most pronounced in the first summer after the assumed eruption date (see the methods section). The response pattern is even more pronounced when using the 17 preindustrial eruptions [−2.3°C (±0.9°C)], although it is absent for the three eruptions that occurred after 1850 [0.0°C (±0.7°C)] (Table 4). The reconstructed average summer scPDSI increase following 20 large volcanic eruptions is 2.21 scPDSI. Although the reconstructed pluvials are strongest in the first year after the eruption, there is a tendency toward prolonged precipitation surplus. A similar response pattern is found when using 17 eruptions before 1850 (2.43 scPDSI), but the response is clear for the three eruptions that occurred during Europe’s industrial period since 1850 (0.92 scPDSI). Spatial characteristics of the reconstructed summer hydroclimate emphasize a meridional structure in the postvolcanic response, with wet conditions in the Mediterranean region and northeastern Scandinavia in contrast to relatively dryer conditions in central Europe (Fig. 9). Distinct pluvials on the Iberian Peninsula, as well as over large parts of the Mediterranean basin, are evident in the first and second year after an eruption. Above average scPDSI values after large eruptions are also found over northern Fennoscandia (see the discussion section for details), whereas most of central Europe and the Baltic experience drought after eruptions.
4. Discussion and conclusions
After nearly one decade of collaboration between institutions in Spain, Switzerland, and Germany, the world’s best replicated MXD chronology of living and relict pines (Pinus uncinata) from a unique high-elevation site in the Spanish central Pyrenees–Lake Gerber (Büntgen et al. 2008) has been developed. The updated composite record now covers the period 1186–2014 CE with a minimum of eight and a maximum of 175 MXD measurement series per year. After standardization and trend preservation, the well-replicated PMR chronology correlates significantly positive (p ≤ 0.01) with MJ&AS temperatures over the Iberian Peninsula and parts of northern Africa. The new summer temperature reconstruction for the western Mediterranean basin and the last eight centuries shows a recent slowdown of warming after 2003 (Fyfe et al. 2016), as well as exceptional cooling following the 1257 eruption of Samalas volcano in Indonesia (Lavigne et al. 2013; Vidal et al. 2015). The new Pyrenees record describes a larger amplitude of low-frequency temperature variability than four climate model simulations for the same region (CMIP5; Taylor et al. 2012), which tend to consistently underestimate natural temperature variability and overestimate the recent warming trend in the western Mediterranean region. Importantly, the new MJ&AS temperature reconstruction reveals significant (p ≤ 0.01) covariability with a regional subset of the OWDA hydroclimate reconstruction at interannual–multidecadal time scales (Cook et al. 2015, 2016). Higher summer temperatures throughout the last eight centuries have in general resulted in an increasing risk for summer drought over the Iberian Peninsula.
The subsequent discussion, split into five sections, places our results in a broader paleoclimatic context and addresses selected data-related and methodological issues.
a. Proxy–target association
Although the new PMR chronology exhibits variability from interannual to multicentennial time scales, its nonsignificant first-order autocorrelation structure persists throughout the last eight centuries (Fig. 2c). Thus, the new MJ&AS temperature reconstruction is most certainly not biased toward an artificial amplification of lower-frequency variability due to tree physiological carry-over effects from previous growing seasons (Frank et al. 2007). That the new PMR chronology has a third less autocorrelation than the local OWDA reconstruction is thoroughly acceptable given the nature of scPDSI’s calculation (Cook et al. 2015), one that includes a number of month-to-month, even seasonal, carry-over affects (Gilman et al. 1963; Wells et al. 2004). The so-called biological memory that can potentially shift tree-ring-based climate reconstruction from “white” to “red” frequency domains (Frank et al. 2007; Bunde et al. 2013) plays a minor role in MXD formation in Pinus uncinata from the upper treeline >2200 m MSL in the Spanish central Pyrenees (Camarero et al. 1998; Büntgen et al. 2008, 2012b; de Andrés et al. 2015). It is important to note that the first-order autocorrelation structure of the new PMR chronology and instrumental measurements are nearly identical over their common period (0.12 vs 0.10, respectively) (Fig. 5a), which confirms the skill of the PMR chronology to capture year-to-year and longer-term changes in regional summer temperature.
Although the first-order autocorrelation structures of the MXD proxy and instrumental target time series are practically the same, a decline in reliable meteorological measurements across the Pyrenees is evident as early as the first half of the twentieth century (Büntgen et al. 2008, and references therein). Moreover, the available gridded products that effectively interpolate information from much of Europe to provide early data for the greater Pyrenees region are characterized by a substantial reduction of temperature variability back in time (Fig. 3). Since it is very unlikely that the quality of the PMR chronology is changing during the last 100–200 yr when sample size and growth coherency are stable (Fig. 2), the observed weakening in temperature sensitivity before ~1950 is probably related to an overall lower quality of the gridded temperature records for this region (Jones 2016). Simply scaling over the shorter 1950–2014 period thus appears most appropriate to preserve the full amplitude of MJ&AS temperature variability (Esper et al. 2005). July temperatures were excluded from the seasonal calibration window, because they do not play a role on MXD formation for Pinus uncinata at Lake Gerber (Fig. 4).
b. Growth response
Effects of summer drought, together with the seasonal course of tree metabolism, are probably relevant for the observed seasonal response pattern of the PMR chronology to early and late summer temperatures (Fig. 4). While May and June coincide with the reduction of stored nonstructural carbohydrates (NSCs) to form early wood (EW) and the beginning of latewood (LW) formation, the LW lignification is mainly restricted to August and September (Camarero et al. 1998). The absolute lignin content of the cell walls, however, can be neglected when interpreting the Walesch-based relative MXD values (Schweingruber et al. 1978), which are a simple reflection of the highest density value within the LW zone of an annual tree ring. The PMR chronology represents the ratio of cell thickness and cell lumen area, for which July temperatures are found to be irrelevant (Fig. 4). A closer look at ring formation, including cell division, cell expansion, cell wall thickening, lignification, and postmitotic senescence (Cuny et al. 2014), is necessary to fully explain the complex response pattern of the PMR chronology. Since such split behavior has typically been characterized in lowland Mediterranean sites where summer drought constrains xylogenesis (Camarero et al. 2010) we are not disturbed by the observed response pattern.
The end of winter dormancy and subsequent onset of the growing season is mainly controlled by species-specific levels of temperature and photoperiodic sensitivity (Basler and Körner 2014). Highest cell formation in alpine and boreal conifers between June and July (Deslauriers et al. 2003; Mäkinen et al. 2003) has been also described for Pinus uncinata at the upper treeline in the Spanish central Pyrenees (Camarero et al. 1998). Secondary growth at the high-elevation Lake Gerber commonly starts around early to mid-June (de Andrés et al. 2015) whereas the maximum daily growth rate, expressed as the number of cells per increment of radial stem enlargement, often continues until the end of July. Until this point of xylogenesis (Cuny et al. 2015), about 80% of the cells, which account for almost the entire EW, have already reached maturation (i.e., the final lignification of the cell walls). The LW development, however, continues during September (Camarero et al. 1998). Cell development depends not only on assimilates of the actual photosynthetic activity, but also on remobilized NSCs, which start being used at the beginning of the EW formation (Deslauriers et al. 2016).
The exceptionally high correlations of the PMR chronology to May temperatures—before LW is developed—might be explained by the ability of the evergreen Pinus uncinata to be photosynthetically active before the onset of secondary growth to produce and store NSCs (Fig. 4), which are ultimately needed to produce EW and then develop LW from mid-June to late July (Kirdyanov et al. 2007). The positive relationship between November temperatures of the year prior to tree growth and the ring width measurements of Pinus uncinata suggests the importance of carry-over effects (Frank et al. 2007). Such products and processes most likely include the climate-induced production of photosynthates during periods of cambial dormancy on EW production at the beginning of the following year (Tardif et al. 2003). A slightly different strategy has been reported for deciduous trees flushing their leaves before the onset of secondary growth, which remobilize NSCs from the previous growing season (Gessler and Treydte 2016). In the case of Pinus uncinata growing at high-elevation sites in the Spanish central Pyrenees, climatic influences of the year prior to the MXD formation, however, can be ignored.
Reduced biological memory in the case of MXD is further reflected by the low first-order autocorrelation of the PMR chronology. This rather white behavior supports our hypothesis that most of the carbohydrates are produced at the onset of the current vegetation period to the EW formation (Camarero et al. 1998). To realize maximum growth rates and start developing LW at the end of June and July, Pinus uncinata at Lake Gerber most likely utilizes the NSCs synthesized in May and June. This strategy would further explain the lack of correlations with the high summer temperature means and precipitation totals of July. The quasi early-late split warm-season MXD response to MJ&AS temperatures also reinforces the hypothesis that NSCs are not only highly mobile but also influence cell development throughout the entire vegetation period. The NSCs that are needed to produce EW and LW cell walls are generated from a mixture of seasonally “recent” and “old” NSCs stored in parenchyma cells (Kagawa et al. 2006). This process most likely starts as early as the beginning of the growing season in May and June, because photosynthetic activity is not sufficient to produce tissues without remobilized NSCs. Considering the effects of temperature and hydroclimate in all seasons, future research should focus the improvement of process-based, mechanistic models of continuous tree growth.
c. Proxy update
The fundamental value of all prehistoric climate records from proxy evidence is the long-term perspective they provide. Although there are hundreds, maybe thousands, of annually resolved, tree-ring reconstructions, there are only a handful of millennial length or longer that also capture a reasonable temperature signal (Esper et al. 2016). Toward understanding the “color” of climate change only the longest records are meaningful (Osborn and Briffa 2004). In the not so distant past a prevailing opinion was “Climate changes irregularly, for unknown reasons, on all timescales” (National Research Council 1991, p. 21). The recent temperature decline could be considered as one such example. It may be that some inertia in the climate system, from some yet unknown forcing or combinations of forcings, is responsible for what appears to be a random shift in mean temperatures. What is certain is that 1) the common period within the current global network of temperature measurements (a proxy themselves) is too short to investigate any lags greater than a couple decades, and that will not change for decades to come; 2) when it comes to disentangling natural variability from anthropogenically affected variability the vast majority of the instrumental record may be biased; and 3) of the few suitably long preindustrial period records that do exist, it is likely that perceived inconsistencies in their earliest measurements have been treated or removed. But even they are too short to study lags greater than three to four decades. Therefore, it is imperative we continue the pursuit of long tree-ring chronologies and continually improve extant collections.
Although the causes of the recently measured slowdown in global and regional warming during the last decade are still debated (Karl et al. 2015; Fyfe et al. 2016), our study provides the first long-term proxy evidence for this temperature decline over the western Mediterranean basin (Fig. 7). This finding is in line with local, regional, and subcontinental meteorological observations (Fig. 5) and consistent with the observations by Gleisner et al. (2015) that the post-2003 pause in rising mean surface temperatures is most strongly expressed at midlatitudes. Potential influences of the North Atlantic Ocean on warm-season temperature variability, including the recent slowdown of warming in the Pyrenees, can be excluded (Fig. 6). The only obvious large-scale teleconnection pattern is the circumpolar wave train (Jungclaus et al. 2014), which is associated with the upper troposphere (200 hPa) geopotential height and meridional velocity (Saeed et al. 2014), and that to a certain degree synchronizes surface air temperatures in Europe and central Asia (Büntgen et al. 2016). Our independent proxy confirmation of a widely debated climatological phenomenon emphasizes the value of continuously updating high-resolution proxy archives (Jones et al. 2009; Tegel et al. 2010). It should be noted in this regard that the Lake Gerber sampling site in the Spanish central Pyrenees still provides some potential to extend the composite chronology further back in time. There still remain ample subfossil trunks preserved in several smaller and larger lakes between 2100 and 2400 m MSL although to recover them will be costly, calling for professional diving campaigns, as well as numerous high-precision 14C dates in addition to conventional tree-ring dating (Sookdeo et al. 2017).
Despite its ability of tracking year-to-year and longer-term summer temperatures since ~1950 (Fig. 5), the PMR chronology somewhat underestimates the 2003 warmth as measured by the instrumental stations. This offset partly relates to the fact that the heatwave mainly occurred in July and August, and chiefly affected central Europe (Schär et al. 2004), whereas in the Pyrenees May and June were characterized by average temperatures. In addition, increased frequency of summer droughts, which have been reported to weaken the growth–temperature relationship even at high elevations in the Pyrenees (Galván et al. 2015), as well as possible changes in the diurnal temperature range (DTR; Easterling et al. 1997) and thus cloudiness (Büntgen et al. 2013), cannot be completely ruled out as complicating factors (Stine and Huybers 2014). The daily difference between maximum and minimum temperatures, the DTR, represents a spatially reliable measure of cloudiness that constitutes an ideal hydroclimatic target for the calibration of tree-ring proxy records (Büntgen et al. 2013). Since the reconstructed MJ&AS temperature range during the modern period is within the envelope of preindustrial variability since at least 1186 CE (Fig. 7a), a possible weakening of temperature sensitivity due to increased drought stress could indeed have occurred whenever past summers have been exceptionally warm and dry.
d. Amplitude preservation
The reconstructed long-term variability exceeds the preindustrial multidecadal–centennial variability in four state-of-the-art climate model simulations (Fig. 7b). This mismatch between the proxy reconstructions and the four model simulations is in line with a general tendency of state-of-the-art climate model simulations to underestimate the amplitude of reconstructed natural low-frequency temperature variability during the last millennium (Bothe et al. 2013; Fernández-Donado et al. 2013; Phipps et al. 2013; Luterbacher et al. 2016; Ljungqvist et al. 2016). Such disagreement might indicate that the role of internal unforced variability is greater than expected (Goosse 2017; Matsikaris et al. 2016), and/or that the climate sensitivity to the prescribed forcings needs adjustment. Moreover, it stresses our limited understanding of the appropriate scaling of solar versus volcanic forcing on different spatiotemporal scales (i.e., the different forcing agents themselves may need adjustments). Although there are some suggestions that past ranges of low-frequency solar variability have been larger (Shapiro et al. 2011), state-of-the-art climate models are usually driven by a relatively dampened amplitude of long-term changes in solar activity (Schmidt et al. 2011, 2012).
Although the spatial resolution of the model output is too coarse to reflect the complex topography of the Pyrenees, we assume that the simulations of near-surface temperature variability are able to capture the basic dynamical background conditions of the region. All four CMIP models reveal positive summer temperature anomalies when high pressure is located over mid-Europe and low pressure is found over the eastern North Atlantic (Fig. 10), leading to southerly winds. Generally cooler summers can be associated with the opposite pattern, a blocking-like situation over the northeastern North Atlantic consistent with negative pressure anomalies over central and eastern Europe leading to northerly winds. Despite its persistence, this pattern reveals some modification due to the position and strength of the pressure anomalies and its seasonal fluctuation.
e. Volcanic forcing
Large volcanic eruptions can perturb hemispheric-to-global-scale climate, principally via influencing atmospheric radiation and dynamics. The spatial and temporal patterns of climate forcing depend on a variety of factors, among them the quantity of sulfur gases released on eruption, the eruption intensity and plume height (stratospheric injection strongly favoring more pronounced impacts), source latitude, season and phases of other phenomena such as El Niño–Southern Oscillation (ENSO) and the quasi-biennial oscillation (QBO; Timmreck 2012; Oppenheimer et al. 2014). The strongest impacts on the global climate system are associated with large, explosive eruptions of sulfur-rich magmas from tropical volcanoes. While the 1991 eruption of Mount Pinatubo (Philippines) has yielded the most detailed insights into the impacts of volcanism on the Earth system, reflecting in part its occurrence in the era of satellite observation, wider perspectives of volcanic forcing on climate have been sought through studies of past eruptions. Although sometimes based on a combination of climate proxy evidence from ice core glacio-chemistry, historical documents, and geological records (Sigl et al. 2015; Nooren et al. 2017), tree-ring chronologies represent the dominate archive (Briffa et al. 1998, 2004; Anchukaitis et al. 2012; Esper et al. 2013a,b, 2015; Büntgen et al. 2016; Schneider et al. 2015; Stoffel et al. 2015).
All studies confirm the essential pattern of cooling over continental regions in the Northern Hemisphere in the first years following climatically significant eruptions. Observational measurements and tree-ring chronologies indicate global-to-hemispheric-scale posteruption temperature anomalies as large as −0.5°C following historical eruptions (Mass and Portman 1989; Briffa et al. 1998; D’Arrigo et al. 2009). However, large-scale atmospheric dynamics can modulate climate response (Briffa et al. 2004; Raible et al. 2016), leading to substantial latitudinal differences (Fischer et al. 2007; Hegerl et al. 2011). Consequently, there is no blueprint for the spatiotemporal response patterns (Wagner and Zorita 2005; Hegerl et al. 2011). The ability of coupled global climate models to accurately simulate the dynamical response to volcanic eruptions is also not straightforward (Stenchikov et al. 2006; Anchukaitis et al. 2012; Driscoll et al. 2012; Zanchettin et al. 2013a,b). Volcanic fingerprints on European hydroclimate (Fischer et al. 2007; Büntgen et al. 2015c), including the effects of stratospheric sulfate aerosols on tropospheric cloud formation, total precipitation, and the diurnal temperature range (Auchmann et al. 2012, 2013; Wegmann et al. 2014; Brugnara et al. 2015), are even less well understood.
Although substantial improvements have been achieved recently in ice core dating (Sigl et al. 2015), temporally varying uncertainties still exist for various episodes throughout the past 2500 years (Büntgen et al. 2017). This putative inconsistency implies that dates for volcanic sulfate layers in ice cores cannot be relied on entirely to provide accurate information on eruption years, even allowing for the time lag between eruption and deposition of aerosol over the polar regions. Mismatches between the timing and amplitude of volcanic sulfate peaks in ice cores from Greenland and Antarctica can also reflect the geographical location and the seasonal timing of an eruption. The fact that most of the events in the Common Era identified from ice core evidence cannot be matched to known historical eruptions also hampers intercomparison of ice core and dendrochronological records (Esper et al. 2013a,b, 2015). Consideration of past changes in the dominant modes of internal climate variability, such as ENSO and the southern annular mode (SAM), which can modulate the impacts of volcanic forcing among different spatiotemporal scales (Emile-Geay et al. 2008), is essentially impossible. Separate analyses of volcanic forcing effects during preindustrial and industrial times may help illuminate the mitigating or amplifying role greenhouse gases and tropospheric aerosols might play.
An individual response assessment of the estimated eruption dates of 17 volcanic events before 1850 and three afterward reveals a wide range of summer temperature anomalies (Table 4). Application of the SEA summarizes the mean response of the new MJ&AS temperature reconstruction to all 20 putative eruption dates (Fig. 8c). A marked 1-yr cooling of ~2°C is obtained from all 20 as well as the 17 eruptions before 1850. The three eruptions afterward though yield a much smaller temperature downturn. The average warm-season cooling of ~2°C, however, most likely represents the lower end of the response range to volcanic forcing in this region, because the 1961–90 baseline period is relatively high, several eruptions occurred during solar minima, and temporal adjustments of the SEA on the basis of MXD depressions might amplify the obtained pattern.
A rapid recovery from the abrupt summer cooling is a typical feature of MXD-based temperature reconstructions (Schneider et al. 2015; Stoffel et al. 2015). A less distinct and more prolonged response pattern is obtained from the ring-width–based OWDA (Fig. 8d). In evaluating spatial response patterns of OWDA to volcanic eruptions (Fig. 9), it should be noted that the persistent anomalous wet conditions over large parts of northern Fennoscandia have possibly been inferred from the inverse relationships with temperature sensitive conifer growth. This assumption is based on the fact that warm-season temperature means in this part of the world are negatively correlated with precipitation totals during the same seasonal window.
A closer look at the response of the new MJ&AS temperature to the putative 1257 Samalas eruption (Lavigne et al. 2013; Vidal et al. 2015) and the 1815 Tambora eruption (Stothers 1984; Oppenheimer 2003) reveals distinct summer cooling (−4.4°C) in 1258 with respect to the 1961–90 baseline period (Table 4), but a less pronounced depression in 1816 (−2.51°C). The summer of 1258 is the coldest between 1186 and 2014, whereas 1816, the so-called year without a summer (Stothers 1984; Luterbacher and Pfister 2015; Raible et al. 2016), is only the 17th coldest summer. In line with the temperature responses, the OWDA reflects a more abrupt precipitation surplus in 1258, compared to a more prolonged pluvial from 1812 to 1821 (not shown).
A detailed view at the responses of the ensemble simulations MPI-ESM-P and CCSM4 (Bothe et al. 2013; Otto-Bliesner et al. 2016) to the Samalas and Tambora eruptions in 1257 and 1815 CE, respectively, provides a better understanding of the role the volcanic forcing data play in the different models (Schmidt et al. 2011, 2012). While the MPI model implements the forcing dataset of Crowley (2000), the CCSM4 model uses the forcing dataset of Gao et al. (2008). The degree of post-Samalas cooling over the western Mediterranean basin varies considerably between the individual ensemble members. A wide range of spatial variability is also evident in the individual simulations at the European scale. As expected, and matched by the new MJ&AS temperature reconstruction, the modeled response to the Tambora eruption shows a less pronounced summer cooling compared with the Samalas event (Trigo et al. 2009; Stoffel et al. 2015; Raible et al. 2016), and some realizations even generated positive temperature anomalies over the Pyrenees. Inconsistency between the simulations points to the importance of internal climate variability, as well as the specifics of how the volcanic forcing was implemented in the respective model. The role of different forcings and choices of climate models on the simulated response to volcanic eruptions needs better understanding (Zanchettin et al. 2016). Reasons for mismatches in the timing and amplitude of reconstructed and simulated postvolcanic summer cooling are manifold, thus calling for an even stronger continuation of interdisciplinary research efforts.
f. Concluding remarks
This study demonstrates the relevance of updating MXD-based climate reconstructions not only back in time but also toward the present. We also advocate that new research be carried out aimed at achieving a better understanding of the biological mechanisms responsible of the formation of the widely used dendroclimatic MXD parameter. Moreover, the new paleoclimatic insight from the western Mediterranean region emphasizes the importance of comparing different proxy archives of various meteorological parameters with model simulations for better understanding regional climate dynamics.
The National Park d’Aigüestortes I Estany de Sant Maurici kindly provided sampling permission. U. Büntgen received funding from the Ministry of Education, Youth and Sports of the Czech Republic within the National Sustainability Program I (NPU I Grant LO1415). F. Reinig was supported by the Swiss National Science Foundation (SNF Grant 200021L_157187/1). Prof. Edward R. Cook and two anonymous referees kindly commented on an earlier version of this manuscript.
Current affiliation: Department of Geography, Masaryk University, Brno, Czech Republic.