2.1 The isotope network

In this study, we investigate the first two dominant modes of variability of 26 δ¹⁸O_cel records, distributed over Europe, and their relationships with regional and large-scale climate anomalies. A total of 22 of the 26 δ¹⁸O_cel records were generated within the EU ISONET (Annual Reconstructions of European Climate Variability using a High-Resolution Isotopic network) project (Treydte et al., 2007a, b). Furthermore, four additional sites from Bulgaria, Turkey, southwestern Germany, and Slovenia were added to the ISONET network for the current study (Heinrich et al., 2013; Hafner et al., 2014). In total, the isotope network contains eight broadleaf tree sites (Quercus) and 18 coniferous tree sites (Pinus, Juniper, Larix, and Cedrus) from altitudes varying for each location from 10 m up to 2200 m above sea level (Fig. 1 and Table S1 in the Supplement). A total of 24 of the 26 sites are distributed over the European continent, whereas 2 additional sites are located in the Atlas Mountains of Morocco and in the Taurus Mountains of Turkey.

Figure 1Spatial distribution of sample sites combined with the corresponding altitude. The highest density of sample sites exists in central and western Europe. The colour indicates the tree type: Cedrus – dark blue, Juniper – yellow, Larix – grey, Pinus – red, and Quercus – light blue. The corresponding elevation (10–2200 m) is shown by the size of the circles.

The stable isotopes of oxygen in tree-ring cellulose, reported as δ¹⁸O_cel vs. Vienna Standard Mean Ocean Water (Craig, 1957), for each site were determined as described by Treydte et al. (2007a, b). At least four dominant trees were chosen per site, and two increment cores were taken per tree. After the standard dendrochronological dating following Fritts (1976), the individual tree rings were dissected from the cores. According to Treydte et al. (2007a, b), all tree rings from the same year were pooled prior to cellulose extraction for the majority of sites. For oak, however, only the latewood was used for the analyses. This procedure assumed that climate signals of the current year were predominantly applied, as the early wood of oaks frequently contains climate information of the preceding year. This is based on the fact that the proportion of the reserves of deciduous trees of the isotope network is higher at the beginning of tree-ring formation compared with the conifers of the isotope network (because they are evergreen). The isotope records have an annual temporal resolution.

The first 100 years of data from the network as well as a general description have already been published (Treydte et al., 2007a, b). Data from individual sites or regional groups of sites have also been published elsewhere (Saurer et al., 2008; Vitas, 2008; Etien et al., 2009; Hilasvuori et al., 2009; Haupt et al., 2011; Saurer et al., 2012; Rinne et al., 2013; Helama et al., 2014; Labuhn et al., 2014; Saurer et al., 2014; Labuhn et al., 2016; Andreu-Hayles et al., 2017). Here, we use the extended ISONET+ product where the longest chronologies cover a period from 1600 to 2005. The highest data density is available for the period from 1850 to 1998: a total of 26 time series available for further analysis. A total of 12 time series cover the entire period of 400 years.

2.2 Climate data

For the gridded climate information, we used the gridded fields of monthly temperature averages and monthly precipitation sums from the Climatic Research Unit (CRU TS) version 4.04 (Harris et al., 2020). Both quantities are derived from the interpolation of monthly climate anomalies from extensive networks of weather station observations. The CRU TS dataset has a spatial resolution of $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ and covers the period from 1901 to 2019.

For the large-scale atmospheric circulation, we make use of the gridded fields of geopotential height 500 mbar (Z500) from the 20th Century Reanalysis Project (20CR) version V2c (Compo et al., 2011). The 20CR reanalysis dataset has a temporal resolution of 6 h and a meridional and zonal resolution of 2^∘. In our study, the ensemble mean of 20CR was used, which was computed from 56 ensemble members. The climate variables are available for the period from 1851 to 2014, and they are provided by NOAA/OAR/ESRL PSL, Boulder, Colorado, USA (downloadable from their website: https://www.psl.noaa.gov/data/gridded/data.20thC_ReanV2c.html, last access: 2 March 2020).

As precipitation and temperature are both important for the δ¹⁸O_cel ratio, we test the relation between the isotope ratio and drought conditions at the European scale. For this purpose, the Standardized Precipitation Evapotranspiration Index dataset from Vicente-Serrano et al. (2010) is used in this study with an aggregation time of 3 months (SPEI3). The SPEI3 index is suitable for this analysis because we take the climate conditions of the preseason into account.

Furthermore, the Extended Reconstructed Sea Surface Temperature version 5 (ERSST5; Huang et al., 2017) is included in the study to investigate the correlation between the global sea surface temperature (SST) and the climate signals of δ¹⁸O_cel. The ERSST5 dataset was created by Huang et al. (2017), and it is derived from the International Comprehensive Ocean–Atmosphere Data Set. The monthly SST fields have a spatial resolution of $\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$, and they are available for the time range from 1854 to present. The gridded fields of the ERSST5 dataset are provided by NOAA/OAR/ESRL PSL, Boulder, Colorado, USA, and they can be downloaded from their website (https://www1.ncdc.noaa.gov/pub/data/cmb/ersst/v5/netcdf/, last access: 15 October 2020). All of the above-mentioned climate data were seasonally averaged: DJF (December to February), MAM (March to May), JJA (June to August), and SON (September to November). Furthermore, the linear trend from each grid cell was removed.

As our study is focused on the relation between ENSO variability and the δ¹⁸O_cel network, we also used the anomaly (1981–2010 mean removed) of the December Niño 3.4 index (HadISST1; Rayner et al., 2003). The index represents the averaged SST from 5^∘ S to 5^∘ N and from 170 to 120^∘ W (downloadable from https://psl.noaa.gov/gcos_wgsp/Timeseries/Nino34/, last access: 3 March 2020). To compare the relation between ENSO variability and our isotope network variability (also in the past), three different reconstruction of ENSO activity have been used in this study: (i) the Dätwyler et al. (2019) reconstruction of the annual Niño 3.4 index for the last millennium (based on a multi-proxy network); (ii) the reconstruction of the NDJ (November–December–January) Niño 3.4 index for the last 700 years by Li et al. (2013) (based on a tree-ring network in the tropics and mid-latitudes), and (iii) an annual reconstruction of ENSO variability based on the North American Drought Atlas (Cook et al., 2004; Li et al., 2011).

2.3 Modelled δ¹⁸O in precipitation and soil water

In addition to the observational and reanalysis-based climate data, we also investigate the relation between δ¹⁸O_cel and δ¹⁸O_Pandbetweenδ¹⁸O_celandδ¹⁸O_SW, as simulated by the water-isotope-enabled ECHAM5-wiso model (Werner et al., 2011), in order to gain further insights into the correlation with other seasons and to identify a multi-seasonal climate signal. The ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{P}}/{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{SW}}$ dataset was created by Butzin et al. (2014), who used the isotope-enabled version of the ECHAM5 atmospheric general circulation model (Roeckner et al., 2003; Hagemann et al., 2006; Roeckner et al., 2006) which is called ECHAM5-wiso (Werner et al., 2011). In Butzin et al. (2014), values of present-day insolation and greenhouse gas concentrations (Intergovernmental Panel on Climate Change, 2000) and monthly varying fields of sea surface temperatures and sea ice concentrations according to ERA-40 and ERA-Interim reanalysis data (Uppala et al., 2005; Berrisford et al., 2011; Dee et al., 2011) are used to force the model. To represent the climate conditions of the period from 1960 to 2010, Butzin et al. (2014) used an implicit nudging technique (Krishnamurti et al., 1991; the implementation in ECHAM is described by Rast et al., 2013). The nudging technique is a part of climate modelling sciences, where modelled fields of climate variables are relaxed to observations or data from reanalyses. In the study of Butzin et al. (2014), the modelled fields of surface pressure, temperature, divergence, and vorticity are coupled to ERA-40 and ERA-Interim reanalysis fields (Uppala et al., 2005; Berrisford et al., 2011; Dee et al., 2011) for the period from 1960 to 2010. The monthly grids of ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{P}}/{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{SW}}$ have a horizontal grid size of approximately $\mathrm{1.9}{}^{\circ }\times \mathrm{1.9}{}^{\circ }$. As our study is focused on seasonal variability, we have computed the seasonal averages for ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{P}}/{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{SW}}$ (DJF, MAM, JJA, and SON).

2.4 Data analysis

As a first step, the characteristics of each time series of the δ¹⁸O_cel network and their relation to altitude and latitude are investigated. For a better comparison, the linear trend of each δ¹⁸O_cel time series is removed and the time series are standardized (z values).

To combine the signals of the isotope network, we use principal component analysis (PCA) and empirical orthogonal functions (EOFs). These techniques have been described by Pearson (1902) and Hotelling (1935) and were used for the first time by Lorenz (1956) for climatological studies (Storch and Zwiers, 1999). By applying PCA, it is possible to extract a common climate signal from the δ¹⁸O_cel network that explains the highest part of the variability of the input dataset. This is done by rotating the initial data onto axes that are orthogonal to each other (Schönwiese, 2013) by the corresponding eigenvectors. Therefore, the eigenvectors are used as a transformation matrix. The corresponding analysis of the eigenvector values is known as empirical orthogonal function (EOF) analysis. The goal of this analysis is to identify the most dominant patterns of the δ¹⁸O_cel tree network variability, which explain a significant part of the variance for a specific region. The largest part of the variance can be explained by the pattern of the leading EOF. The temporal perspective on these patterns is given by the principal components (PCs), which describe the phase and the amplitude.

The resulting components are further checked to establish if they fulfil the requirements of the rule of North et al. (1982). This rule states that the pattern of the eigenvectors of one component is strongly contaminated by other EOFs that correspond to the closest eigenvalues (Storch and Zwiers, 1999). To determine whether two consecutive patterns can be interpreted as distinct patterns, it is necessary to calculate the standard error in the estimation of the eigenvalues; according to North et al. (1982), this is approximately

where λ is the eigenvalue, and n is the number of degrees of freedom of the dataset. In case that the eigenvalues of two EOFs are not more separated from each other than this standard error, it is unlikely that the two consecutive patterns can be interpreted as distinct from each other, as any linear combination of the two eigenvectors is equally significant. In our study, we only consider the eigenvectors where the successive eigenvalues are distinguishable.

The ISONET network consists of a multi-site and multi-species tree-ring network more or less covering the period from 1600 to 2003. However, some tree-ring series cover the whole period, whereas other series cover only a shorter period. In order to be able to have a long-term perspective, one needs to find a statistically meaningful way to extend the shorter records to make use of the whole 400 years of data. As most multilinear principal component analysis algorithms do not work with gaps in the initial matrix, we make use of an algorithm developed by Josse and Husson (2016) that is able to fill the temporal gaps without a change in the PC. In the first step, we place the mean in the gap and execute a PCA with this dataset. Afterwards, the dataset is projected onto the new component axis so that the values are rotated and the value of the gap changes. The new value of the gap is placed into the initial dataset. With this new dataset, a PCA is again carried out. This process is repeated until convergence is reached. The result is a gap-free dataset that can be used for PCA. To quantify if our results are influenced by the gap-filling method, we tested the correlation between PC1 based on the ISONET network and the first four PCs of the Old World Drought Atlas (OWDA; Cook et al., 2015) for the period from 1850 to 2005 and from 1600 to 2005. If the filling algorithm altered the representation of climate signals over a longer time period, we would expect that the strength of correlation is changing.

Composite maps of average precipitation, air temperature, geopotential height 500 mbar (Z500), and SST are computed using years where the network's principal components are above or below a certain threshold. In our study, we choose maxima events above ($>\mathit{\mu }+\mathit{\sigma }$) and minima events below ($<\mathit{\mu }-\mathit{\sigma }$) 1 standard deviation (σ) with respect to the mean (μ). The composite maps allow us to analyse the general climate state occurring at times of minima (low) or maxima (high) separately. In addition, both composites of climate conditions can be combined under the assumption that the minima and maxima events show the opposite climate state. For this purpose, the minimum composite is subtracted from the maximum composite at each grid point (high−low). In addition to the composite maps, we extract the values of PC1 for those years for which ENSO values are higher than the average plus 1 standard deviation and lower than the average minus 1 standard deviation. The difference from the former distribution for the values of minima and maxima years is tested with a t test. To better understand if El Niño or La Niña events co-evolve with extremes in the δ¹⁸O_cel time series, event coincidence analysis (ECA) (Donges et al., 2016; Siegmund et al., 2017) using PC1 and a December Niño 3.4 index is applied (HadISST1; Rayner et al., 2003). ECA quantifies the simultaneity of events contained in two series of observations that can be computed with the “CoinCalc” R package (Siegmund et al., 2017). Furthermore, the CoinCalc package provides functions to test if the coincidences are significant. In our study, we analyse whether the years in which the Niño 3.4 index is above (below) the 75th (25th) percentile match the 75th (25th) percentile in PC1. In general, a significance level of α=0.05 was used in all analyses.

Finally, we analyse the relation between seasonally averaged δ¹⁸O_P and δ¹⁸O_SW from winter, spring, and summer, based on nudged ECHAM5-wiso simulations (Butzin et al., 2014), with PC1 based on the δ¹⁸O_cel values from the ISONET network. Our goal is to test the correlation between δ¹⁸O_cel and modelled ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{P}}/{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{SW}}$ to identify if the water, which is used within the photosynthesis processes, has a multi-seasonal isotopic signature.

3.1 Characteristics of the δ¹⁸O_cel network

The highest mean δ¹⁸O_cel values occur in the southern locations (i.e. in Turkey and Cazorla in Spain). For oaks and pines, the lowest mean values are found for the northern sites (Fig. 2a). Moreover, generally lower δ¹⁸O_cel values are identified for Quercus compared with Pinus. Angiosperm wood tissue contains vessels (i.e. specialized water-conducting cells) that are generally larger in diameter and, therefore, more conductive to water than conifer wood cells (Sperry et al., 2006; Carnicer et al., 2013). The overall variance in the datasets is not dependent on the type of tree species. The δ¹⁸O_cel time series from Pöllau in Austria is characterized by the highest standard deviation, whereas the time series from Lochwood in Great Britain shows the lowest standard deviation.

Figure 2The characteristics of the European δ¹⁸O time series/network. Panel (a) describes each time series using boxplots, which are firstly arranged by tree type and secondly by the average value. Additionally, the relation between the average value of each time series is plotted against their latitudinal (b) and their altitudinal position (c).

Download

As the δ¹⁸O source values and the fractionation processes are temperature dependent, it is necessary to evaluate the influence of altitude and latitude on the oxygen isotope ratio. The relation of the δ¹⁸O_cel values with respect to the altitude and latitude of each site is shown in Fig. 2b and c respectively. In Fig. 2c, the linear relationship between the average δ¹⁸O_cel values of the locations and the corresponding latitudes are plotted. As shown by the boxplots in Fig. 2a, the southern sample sites are characterized by the highest average δ¹⁸O_cel values, whereas the northern sites show the lowest average δ¹⁸O_cel values. The relation between the δ¹⁸O_cel values and latitude yields a significant linear regression.

In addition to the latitudinal effect, altitude also influences the oxygen isotope ratios (as shown in Fig. 2b) which can likewise be described by a significant linear regression. It should be noted that the southern sites are found at higher altitudes than the northern sites; thus, the latitudinal and altitudinal gradients may have confounding effects on δ¹⁸O_cel. Therefore, we show that the δ¹⁸O_cel network is influenced by a latitudinal and an altitudinal effect.

3.2 Characteristics of the principal components

Based on the EOF analysis, the first component of the isotope network explains 16.2 % of the variance, the second explains 9.1 %, the third explains 6.4 %, the fourth explains 5.5 %, and the fifth explains 5.2 %. Therefore, the first five components explain a cumulative variability of around 43 %. As the first two components also fulfilled the requirements of the rule of North et al. (1982), they are investigated with respect to their temporal and spatial characteristics.

The dominant pattern (EOF1), which describes 16.2 % of the total variance, shows a spatially homogeneous structure (Fig. 3a). The majority of time series in Europe are characterized by negative eigenvector values. The pole of this EOF pattern is centred over France and Germany. In contrast, tree sites close to the Mediterranean Sea and the northernmost site in Finland are characterized by eigenvector values close to zero. Therefore, these locations contribute to a lesser extent to the first component's time series (PC1) shown in Fig. 3c. The time series of PC1 is characterized by an underlying negative trend from the 17th to the 18th century. From the beginning of the 18th century, a positive trend is observed until the beginning of the 19th century, where the highest values are reached. For the last 150 years, no clear trend is visible.

Figure 3Spatial and temporal variability of the first two δ¹⁸O components and EOFs: (a) EOF for the first δ¹⁸O component (16.2 % explained variance); (b) EOF for the second δ¹⁸O component (9.5 % explained variance). Panels (c) and (d) show the time series for the first and second δ¹⁸O components. The dashed red lines indicate the standard deviation for the years 1600–2005.

The second EOF (EOF2) is characterized by a completely different spatial pattern (Fig. 3b). Negative eigenvector values are found around the North and Baltic Sea, with the smallest eigenvector located in Norway, while positive eigenvector values are identified over the southern/southeastern Europe and the Alpine region. The highest eigenvector values are recorded for the Italian site. In summary, this component highlights a dipole-like structure between northern and southern/southeastern Europe. This dipole-like structure is a well-known feature of the European hydroclimate (Ionita et al., 2015). The time series of the second component (PC2) is characterized by an underlying positive trend from the middle of the 19th century onwards. Furthermore, the highest interannual variability is found for the beginning of the 18th and 19th centuries. The highest values of the second component (PC2) are reached at the beginning of the 18th century, whereas the smallest values are identified for the beginning of the 19th century (Fig. 3d).

Furthermore, we tested the correlation between PC1 and the first four PCs for the summer drought reconstruction based on the Old World Drought Atlas (OWDA; Cook et al., 2015) for the period from 1850 to 2005, where we have the highest sample density. The highest correlation (R=0.43, p value $=\mathrm{3.1}\times {\mathrm{10}}^{-\mathrm{08}}$) is computed with PC2 of the OWDA (EOF plot and time series are shown in Fig. S1 in the Supplement), which explains 16.1 % of the total variance of the OWDA. If the filling algorithm altered the representation of climate signals over a longer time period, we would expect that the strength of the correlation would also change. However, the correlation changes only slightly to a Pearson's R=0.39 (p value $=\mathrm{4.2}\times {\mathrm{10}}^{-\mathrm{16}}$) for the entire period from 1600 to 2005, which indicates that the filling algorithm does not impact the results. As a result, the climate signals obtained are robust and are presented in a similar manner as those for the period with a high sample coverage (i.e. 1850–2005). This argument is also supported by the comparison of the two time series in Fig. S1. We note that, when using the gap-filling method, uncertainties might arise for the first century where the sample density is low. It is, therefore, important to mention that the interpretation for the first century needs to be handled with care, and statements should be regarded as less robust.

3.3 Climate signals of PC1

For the winter season, high values of PC1 are associated with significant warm SSTs in the equatorial Pacific and on the west coast of North and South America (Fig. 4a), which indicates that the high values of PC1 are co-evolving with the occurrence of the El Niño conditions. The strong signal in the Pacific persists throughout spring (Fig. 4c) and reduces in summer (Fig. 4e). A significant warming of the tropical and North Atlantic is visible in spring and summer (Fig. 4c, e). According to Latif and Grötzner (2000), the lagged warming of the equatorial Atlantic can be observed up to 6 months after an El Niño event. In contrast, the northeastern Atlantic and the Mediterranean Sea are characterized by significantly colder SSTs in all three seasons. The low composite maps for PC1 and SST (Fig. 4b, d, f) show features of La Niña conditions associated with colder than average SSTs. This pattern is particularly prominent in the tropical Pacific during winter and spring.

Figure 4Composite maps (high and low) of SSTs related to the first δ¹⁸O component for the DJF, MAM, and JJA seasons. Panels (a) and (b) show the characteristics of the climate in DJF, panels (c) and (d) show the characteristics of the climate in MAM, and panels (e) and (f) show the characteristics of the climate in JJA. The results for maxima (a, c, e) and minima (b, d, f) events of PC1 are also given. The SSTs in winter, spring, and summer are characterized by ENSO activity. Furthermore, the significance is shown using a black grid overlay. The SST dataset from ERSST (Huang et al., 2017) is included in this figure for the period from 1854 to 2005.

The related large-scale atmospheric circulation is shown in the composite map of Z500 (Fig. 5). During high values of PC1 which are co-evolving with the occurrence of El Niño conditions, the atmospheric circulation over Europe is characterized by a low-pressure regime in winter, whereas high-pressure regimes can be identified over the northwestern Atlantic as well as east of the low-pressure system (Fig. 5a). This fits well with the composite map of air temperature in Fig. 6b which shows significant cold conditions over northern Europe and significant warm conditions over southeastern Europe. The described temperature pattern shows similarities to the effects of a negative phase of the winter NAO. The opposite climate state can be observed for low values of PC1 in winter (Figs. 5b, 6b).

Figure 5Composite maps (high and low) for Z500 related to the first δ¹⁸O component for the DJF, MAM, and JJA seasons. Panels (a) and (b) show the characteristics of the climate in DJF, panels (c) and (d) show the characteristics of the climate in MAM, and panels (e) and (f) show the characteristics of the climate in JJA. The results for maxima (a, c, e) and minima (b, d, f) events of PC1 are also given. The Z500 maps show similar characteristics in winter and spring, whereas a pressure regime is directly located over central Europe in summer. The Z500 dataset from 20CRv2c (Compo et al., 2011) is used in this figure for the period from 1901 to 2005. Furthermore, the significance is shown using a black contour line.

Figure 6Composite maps (high−low) for precipitation and air temperature related to the first δ¹⁸O component for the DJF, MAM, and JJA seasons. Panels (a) and (b) show the characteristics of the climate in DJF, panels (c) and (d) show the characteristics of the climate in MAM, and panels (e) and (f) show the characteristics of the climate in JJA. The results for precipitation (a, c, e) and air temperature (b, d, f) are also given. The precipitation and air temperature dataset from CRU TS is included in this figure for the period from 1901 to 2005. The significance is shown using a black grid overlay.

The atmospheric circulation remains in a similar configuration over Europe and the North Atlantic for low and high composite maps for spring (Fig. 5c, d) compared with the winter season (Fig. 5a, b). One important difference from the preseason is that a pressure belt is visible between Europe and the Gulf of Maine in both maps.

The composite maps of Z500 for summer are characterized by a pressure regime centred over France and Germany (Fig. 5e, f). In both maps, it is visible that the low-/high-pressure regime in central and western Europe is surrounded by opposite pressure regimes. In case of a high-pressure system in the centre, it leads to a blocked zonal flow which is shown in Europe in Fig. 5f. The composite maps of precipitation and temperature also support the analysis with the Z500 data. The climate of central and western Europe is characterized by significantly higher (lower) precipitation in central Europe as well as significantly lower (higher) surface air temperatures in summer corresponding to low (high) δ¹⁸O_cel values. The significant relation between PC1 and the summer hydroclimate is also resampled by the correlation between the SPEI3 index for JJA (Vicente-Serrano et al., 2010; longitude from −5 to 10^∘ and latitude from 46 to 52^∘) and PC1. The correlation is significant (R=0.49; p<0.01) for the period from 1901 to 2005 which suggests that PC1 can capture the hydroclimate variability in summer.

To further investigate the relation between ENSO variability and the PC1, we apply two different statistical approaches. The first approach is to analyse if El Niño and La Niña events are separated in the probability density plots of PC1. During El Niño years, the distribution of PC1 is shifted towards higher values, whereas the opposite occurs during La Niña years (see the Supplement; Fig. S2). According to a t test, both shifts are significantly different compared with the distribution of PC1 (p<0.05). The second statistical approach investigates if extremes in δ¹⁸O_cel time series co-occur with El Niño or La Niña events. For this purpose, we apply event coincidence analysis (Siegmund et al., 2017; Donges et al., 2016) using PC1 and the December Niño 3.4 index (HadISST1; Rayner et al., 2003). Over the period from 1871 to 2005, 41.2 % of high and low extremes in the Niño 3.4 index coincided significantly during winter, with high and low extremes of PC1 (p<0.01). By extending the analysis period from 1750 to 1850, coincidence rates (28 % of the NDJ Niño 3.4 (Li et al., 2013) high and low extremes coincided (not significantly; p>0.1) during winter with high and low extremes of PC1) and the correlations (Fig. S2) between PC1 and ENSO reconstructions Dätwyler et al. (2019) and Li et al. (2011, 2013) are shown to weaken.

3.4 Comparison of δ¹⁸O_cel with modelled δ¹⁸O in precipitation and soil water

By employing nudged climate simulations with ECHAM5-wiso (Butzin et al., 2014), we evaluate how the δ¹⁸O_cel tree signature is related to the modelled δ¹⁸O in precipitation and soil water. A significant correlation between PC1 and the modelled δ¹⁸O_P is shown in the correlation maps for winter, spring, and summer, where central Europe is characterized by a moderate correlation (Fig. 7a–c). A similar pattern can be identified for the correlation between δ¹⁸O_SW and PC1. Compared with our previous analysis, the correlation between these quantities increases from winter to summer where it reaches the maximum correlation (Fig. 7g). As the δ¹⁸O_cel ratio is largely dependent on δ¹⁸O_SW, the relation with δ¹⁸O_SW is stronger compared with δ¹⁸O_P. Overall, the results from Fig. 5 indicate that significant correlations for both quantities can be computed for the entire European region except for the eastern parts.

Figure 7Links between the first δ¹⁸O component and the modelled δ¹⁸O in soil water and precipitation from nudged climate simulations with ECHAM5-wiso (Butzin et al., 2014). Panels (a–d) show the correlation between the first δ¹⁸O component (PC1) and δ¹⁸O in precipitation for winter (a), spring (b), summer (c), and autumn (d). Panels (e–h) are the correlation maps for PC1 and δ¹⁸O in soil water for winter, spring, summer, and autumn. For all maps, the significant grid cells are coloured.

3.5 Further climate signals in δ¹⁸O_cel

In addition to the multi-seasonal signal, the second component of the δ¹⁸O_cel values significantly relates to the summer climate (Figs. 8, 9). A positive (negative) geopotential height anomaly in northern Europe with its centre over the North Sea co-occurs with a negative (positive) Z500 anomaly in southeastern Europe (Fig. 9). This coincides with low (high) temperature in central and north Europe, whereas northeastern Europe is characterized by high (low) temperature (Fig. 8a). The same pattern is also shown in the composite maps for precipitation where a similar pattern is presented (Fig. 8b). Based on these patterns, the temporal distribution of extremes in the PC2 time series indicates that the 19th century has experienced increased dryness in northern Europe and enhanced precipitation in the Adriatic region (Fig. 3d).

Figure 8Composite maps (high−low) for the boreal summer related to the second δ¹⁸O component: (a) surface temperature, JJA; (b) precipitation JJA. The datasets are the same as in Fig. 6. Furthermore, the significance is shown using a black grid overlay.

4.1 Latitudinal and altitudinal dependence of the δ¹⁸O_cel network

The δ¹⁸O_cel ratio is affected by the isotopic composition of the source water (δ¹⁸O_P, δ¹⁸O_SW), which varies according to the latitude and altitude of the sample site (e.g. McCarroll and Loader, 2004). The latitudinal position has an influence on the δ¹⁸O in the atmosphere because of the strong correlation between the temperature and the δ¹⁸O composition of water vapour in the atmosphere (Dansgaard, 1964). In addition, Dansgaard (1964) suggested that the altitude also influences the δ¹⁸O ratio, primarily due to the cooling of air masses as they ascend a mountain which is accompanied by the rainout of the excess moisture (Gat, 2010).

The results from this study show that there is a linear relationship between δ¹⁸O_cel and site altitude (Fig. 2b) and latitude (Fig. 2c), which re-enforces the commonly inferred effects of latitude and altitude on the δ¹⁸O_cel (McCarroll and Loader, 2004). Therefore, our results are in line with other studies that have presented such results in the past, e.g. Szejner et al. (2016).

The effects of altitude and latitude on the photosynthesis process are also visible for other tree-ring-based proxies. For example, according to the studies of Körner et al. (1991), Marshall and Zhang (1994), and Dietendorf et al. (2010), there is a highly significant altitude effect on the tree carbon isotope composition (δ¹³C). As temperature decreases with increasing altitude and the partial air pressure is approximately 21 % lower at 2000 m than at sea level, Körner (2007) argued that these environmental conditions lead to a faster molecular gas diffusion at any given temperature.

However, the distribution of our sample sites across Europe indicate that the isotope network is spatially limited (Fig. 1). For instance, the time series from central and western Europe are overrepresented compared with those from southeastern and northern Europe. In addition, it is imperative to extend the isotope network by collecting more δ¹⁸O_cel records from eastern Europe to improve the validity of our results for this region. Furthermore, the samples were not taken from trees growing at the same altitude, which is critical for identifying a latitudinal effect. In fact, the sampled trees in southern Europe grew at higher altitudes (≥1600 m) compared with the other sample sites, which could also bias the δ¹⁸O_cel ratio. It is, therefore, not viable to compare two adjacent sites located at different altitudes, and it is challenging to distinguish between latitudinal and altitudinal effects.

Another point that needs to be considered is that the δ¹⁸O time series used in this study were created using a pooling approach. Overall, the pooling approach has been tested and proven successful for climate analysis (Treydte et al., 2007b), but the approach also has some weaknesses and is strongly discussed in the literature (e.g. Foroozan et al., 2019). For example, Hangartner et al. (2012) recommend avoiding pooling, as an unsuitable tree cannot be omitted without resampling the whole period. They suggest measuring the trees individually instead and only using pooling when there is a strong correlation between the trees. Nevertheless, the pooling approach was adopted for all of the ISONET sites, allowing, for the first time, for the establishment of a tree-ring stable isotope network of more than 20 sites across Europe in collaboration with several laboratories within a reasonable time frame. Without pooling, the spatial dataset, analysis, and interpretation presented here would not be available to the community.

4.2 Links between ENSO and PC1

The results presented in Sect. 3.3 are an indicator that ENSO activity influences the climate signal of PC1. The reason for this is the described warm and cold SST patterns in the equatorial Pacific and on the west coast of North and South America which are associated with El Niño and La Niña events (Allan et al., 1996). Furthermore, the clear and significant separation between El Niño and La Niña events in probability density functions and the coincidence of high and low values of Niño 3.4 and PC1 support our argumentation.

Moreover, the temperature pattern described in Fig. 6b resembles the effects of the winter NAO (Fig. 6b). During the negative phase of the NAO, northern Europe experiences cold conditions and southeastern Europe experiences warm conditions in winter. Studies based on observations (Fraedrich and Müller, 1992; Fraedrich, 1994; Brönnimann et al., 2004; Pozo-Vázquez et al., 2005; Brönnimann et al., 2007) and models (Merkel and Latif, 2002; Mathieu et al., 2004) have shown that ENSO variability can influence the winter NAO. Their results showed that an El Niño event leads to a negative phase of the NAO. Based on our argument, we suggest that high values of PC1 co-evolve with the occurrence of El Niño conditions, and low values of PC1 co-evolve with the occurrence of La Niña conditions.

4.3 The stability of the ENSO signal in the isotope network

In our study, we test the correlation between ENSO and the PC1 of the δ¹⁸O_cel network with three different reconstructions (Li et al., 2011, 2013; Dätwyler et al., 2019) and for two different time periods (1750–1849, 1850–1949), as shown in Table S2. Despite the fact that the sample density of the isotope network is relatively high in these two periods, the correlation between the PC1 and the ENSO reconstructions is weaker and non-significant for the period from 1750 to 1849. However, not only does the correlation get weaker, the correlation between a set of different ENSO reconstructions also becomes weaker in the 18th century, which was shown for specific periods in Dätwyler et al. (2019). Dätwyler et al. (2019) also found a consistent teleconnection pattern during the 18th century, which is different to the known teleconnection pattern of ENSO during the instrumental period. Moreover, the 1850s mark the end of Little Ice Age (LIA), when the atmospheric circulation over Europe (Felis et al., 2018) and its teleconnections to ENSO changed significantly (Rimbu et al., 2003). In their study, Felis et al. (2018) found evidence of an abrupt reorganization of the atmospheric circulation over Europe at the end of the LIA, transitioning from predominantly negative phases of the NAO (weakening of westerly winds) to predominantly westerly flow patterns over central Europe. Furthermore, modelling studies (e.g. Henke et al., 2017) have also reported an increased frequency of El Niño during the LIA due to southern displacement of the Intertropical Convergence Zone.

A change in the ENSO characteristics would also have an influence on the teleconnection with the European climate. For instance, Rimbu et al. (2003) investigated coral time series from the northern Red Sea and identified a nonstationary relationship between the tropical Pacific and the European–Middle Eastern climate during the pre-instrumental period (see also Fig. S3). They showed that a Pacific–North-American-teleconnection-like pattern in its negative phase, which is compatible with La Niña conditions, is associated with positive δ¹⁸O anomalies in the Red Sea coral record from the mid-1930s to late 1960s. After 1970, they detected a shift in the teleconnections that leads to the fact that positive anomalies in the Red Sea coral δ¹⁸O are related to El Niño conditions. An unstable relationship between ENSO variability and the climate of Europe is also found in studies based on instrumental data (e.g. Fraedrich and Müller, 1992; Fraedrich, 1994; Pozo-Vázquez et al., 2005) or ocean–atmosphere coupled models (e.g. Raible et al., 2004; Deser et al., 2006; Brönnimann, 2007). The temporally unstable relationship between climate variables and ENSO is not only restricted to Europe and is also present in other regions of the planet (e.g. Álvarez et al., 2015).

Weak or inconclusive correlations between PC1 and ENSO reconstructions could also arise from the fact that the quality of the ENSO reconstructions decreases, which could be based on an overly low number of samples (especially in first years of the reconstruction period) and a non-stationarity of the teleconnection used (Batehup et al., 2015). Furthermore, ENSO reconstructions are mostly trained within the last 150 years and are used for time periods characterized by an absence of instrumental data. Therefore, confident statements can only be made from 1850 onwards, as instrumental measurements of different climate variables are available. Thus, the only possibility to test and analyse the teleconnection before 1850 is to use ENSO reconstructions.

Another reason for the decrease in correlation could be the change in the climate signal of δ¹⁸O_cel. It is important to consider that the climate signal is directly coupled to the limiting factor for tree growth. It is possible that the limiting factor changes over time, which would result in different responses to climate. Esper et al. (2017) showed that the climate signals in δ¹⁸O_cel and δ¹³C_cel change during warm and cold periods for trees in Switzerland, and they proposed splitting the calibration between these two periods or using corresponding transfer models. Therefore, future research is required to investigate the climate signal during warm and cold periods, as well as the influence on our PCA results. However, we can only suggest that the relationship between ENSO and the European climate may not be stable over time.

4.4 The winter climate signal in δ¹⁸O_cel

The exact mechanism through which the δ¹⁸O_cel captures a climate signal of the preseasons is still debated. For example, Heinrich et al. (2013) mentioned that winters with very low temperatures may damage the cambium more than usual, requiring a longer recovery period. Such winters may have a negative effect on the cambial activity and on the photosynthetic process. Furthermore, Vaganov et al. (1999) showed that precipitation during winter can sustainably affect tree growth in the following year. Their findings are similar to Treydte et al. (2006), who showed that δ¹⁸O_cel can contain a winter signal. Treydte et al. (2006) further argued that, depending on the root system, winter snow fall, and the characteristics of groundwater reservoirs, it is likely that trees use precipitation from the preseasons. Nevertheless, it is possible that winter climate conditions can also be memorized in δ¹⁸O_cel through different climate feedback processes. For example, Ogi et al. (2003) highlighted that a positive NAO is frequently followed by higher pressures and warmer temperatures in Europe during the next summer. The authors suggested that SST, sea ice extensions, and snow fall anomalies capture the winter climate conditions and influence the summer climate.

However, the oxygen isotope signal in cellulose depends primarily on the corresponding oxygen signal of the soil water and precipitation. According to Saurer et al. (2012), δ¹⁸O_SW constitutes the average δ¹⁸O input to the arboreal system over several precipitation events, and it is modified by partial evaporation from the soil (depending on soil texture and porosity) and by a potential time lag, depending on the rooting depth. Here, we obtain the strongest correlations with δ¹⁸O_P in winter, spring, and summer (Fig. 7a, b, c). Because the correlations with δ¹⁸O_SW are strongest in summer and autumn (Fig. 7g, h), and δ¹⁸O_SW is the input of the arboreal system, we suggest that the isotopic signal of δ¹⁸O_cel corresponds to an average over δ¹⁸O_P events from winter, spring, and summer, transferred via the δ¹⁸O_SW. This may also explain the reason behind the strong ENSO signal that δ¹⁸O_cel is able to capture during winter. Moreover, it indicates the high potential of δ¹⁸O_cel to capture climate signals even outside of the growing season.

When considering the multi-seasonal signal, one has to account for the fact that this can lead to certain limitations in the analysis of δ¹⁸O_cel. This is because there is only one δ¹⁸O_cel value in the network used for each site and for each year; however, this value does not represent a clear seasonal signal referring to one season, instead representing a mix of signals referring to several seasons. Overall, the seasonal contributions to the climate signal can vary in strength. Due to the combinations of these influencing factors, the correlation with individual seasons is certainly weakened, as no clear seasonal signal can be represented by the δ¹⁸O_cel value. On the other hand, this enables the observation of several seasons and an understanding of which season has an influence on the biochemical processes in the tree.

Figure 9Composite map (high−low) for Z500 related to the second δ¹⁸O component (PC2) for the boreal summer. The Z500 dataset from 20CRv2c (Compo et al., 2011) is used in this figure for the period from 1901 to 2005. Furthermore, the significance is shown using a black contour line.

4.5 Links between PC2 and large-scale atmospheric modes

In the composite maps for the PC2 (Fig. 9), we obtain a dipole structure between northern Europe and the Mediterranean region. The dipole is characterized by a pressure anomaly centred on the North Sea that expands from the northeastern Atlantic to the Baltic Sea. Its counterpart is present in the northern Mediterranean region, especially in Italy, the northern parts of Greece, and the Adriatic region.

The described pattern shows similarities to the summer European blocking pattern (Barnston and Livezey, 1997; Cassou et al., 2005) which is often associated with the Summer North Atlantic Oscillation (SNAO; Hurrell and van Loon, 1997). According to Cassou et al. (2005), 17.8 % of the positive phase and 17.9 % of the negative phase of the summer European blocking pattern influence the total summer weather regimes in Europe.

The summer European blocking pattern is a surrogate indicator for storm track activities. During the positive index phase, the storm track moves further northwards (Folland et al., 2009; Lehmann and Coumou, 2015). This results in a low storm activity over northern Europe that is characterized by dry conditions, less cloudiness, high temperatures, and a blocked cyclonic flow (Lehmann and Coumou, 2015). On the other hand, the Mediterranean region is affected by lower temperatures and more precipitation. The opposite phenomenon can be identified for the negative index phase. Northern Europe experiences an enhanced storm activity through the southward movement of the storm track over northwestern Europe (Folland et al., 2009; Lehmann and Coumou, 2015), which leads to higher precipitation, higher cloudiness, and lower temperatures. In contrast, the northern Mediterranean experiences dry and warm conditions. These predominant summer European blocking pattern features are well represented in our composite maps of precipitation and temperature for PC2 (Fig. 8). Based on the temporal evolution of PC2, we suggest that there is a tendency towards a negative index phase starting at the beginning of the 20th century.

In addition to the link to summer blocking activity in Europe, the geopotential height pattern is often used in another context. For example, Sillmann and Croci-Maspoli (2009) showed that a positive geopotential height anomaly over the North Sea describes an atmospheric-blocking-like pattern that relates to climate extremes like floods and droughts in the European mid-latitudes. Moreover, this circulation anomaly pattern has been identified as the main driver for extreme dry periods over the eastern Mediterranean (Oikonomou et al., 2010; Rimbu et al., 2014; Ionita and Nagavciuc, 2020) and for summer air temperature variability in Greece (Xoplaki et al., 2003a, b).