A past and present perspective on the European summer vapor pressure deficit

The response of evapotranspiration to anthropogenic warming is of critical importance for the water and carbon cycle. Contradictory conclusions about evapotranspiration changes are caused primarily by their brevity in time and sparsity in space, as well as the strong influence of internal variability. Here, we present the first gridded reconstruction of the summer (June, July, and August) vapor pressure deficit (VPD) for the past 4 centuries at the European level. This gridded reconstruction is based on 26 European tree ring oxygen isotope records and is obtained using a random forest approach. According to validation scores obtained with the Nash–Sutcliffe model efficiency, our reconstruction is robust over large parts of Europe since 1600, in particular for the westernmost and northernmost regions, where most tree ring records are located. Based on our reconstruction, we show that from the mid-1700s a trend towards higher summer VPD occurred in central Europe and the Mediterranean region that is related to a simultaneous increase in temperature and decrease in precipitation. This increasing summer VPD trend continues throughout the observational period and in recent times. Moreover, our summer VPD reconstruction helps to visualize the local and regional impacts of the current climate change, as well as to minimize statistical uncertainties of historical VPD variability. This paper provides also new insights into the relationship between summer VPD and large-scale atmospheric circulation, and we show that summer VPD has two preferred modes of variability, namely a NW–SE dipole-like mode and a N–S dipole-like mode. Furthermore, the interdisciplinary use of the data should be emphasized, as summer VPD is a crucial parameter for many climatological feedback processes in the Earth's surface system. The reconstructed summer VPD gridded data over the last 400 years are available at the following link: https://doi.org/10.5281/zenodo.5958836 (Balting et al., 2022).

1 Introduction

Evapotranspiration is a critical factor for understanding the links and feedback between atmospheric CO₂ and global climate (Good et al., 2015; IPCC, 2021b, a). Within the terrestrial water fluxes, vegetation-produced transpiration represents the dominant factor (Good et al., 2015; IPCC, 2021b, a). One key driver for such vegetation resources and dynamics is vapor pressure deficit (VPD), defined by the difference between the water vapor pressure at saturation and the actual water vapor pressure (Good et al., 2015; IPCC, 2021b, a). The VPD represents the atmospheric evaporative demand, which has an influence on the leaf-level transpiration of terrestrial vegetation and the corresponding stomatal conductance (Grossiord et al., 2020). With increasing VPD, the stomata are closing to minimize water loss (Running, 1976) due to the high atmospheric evaporative demand. Consequently, a minimal stomata opening decreases stomatal conductance and photosynthetic activity (Fletcher et al., 2007).

Vapor pressure deficit (VPD), a fundamental meteorological parameter, exerts a profound influence on diverse aspects of agriculture, ecosystems, and the environment as a whole (Novick et al., 2016). VPD plays a pivotal role in regulating plant growth and health by influencing transpiration rates, stomatal conductance, and overall water uptake mechanisms, which in turn dictate plant morphology, physiology, and productivity (Restaino et al., 2016). Elevated VPD is associated with increased forest mortality (Williams et al., 2013) and enhanced drought susceptibility (Dai, 2013). Furthermore, VPD plays a critical role in crop production and quality. Maintaining an optimal VPD range is essential for maximizing photosynthetic efficiency, nutrient uptake, and ultimately crop yield. Deviating from this optimal range can lead to stress-induced physiological disorders, reduced yields, and inferior produce quality (Zhao et al., 2017). Moreover, high VPD values are linked to elevated wildfire risk (Seager et al., 2015). Understanding VPD dynamics is crucial for wildfire prediction, management, and mitigation strategies aimed at minimizing the destructive impacts of wildfires on both natural ecosystems and human settlements. Changes in VPD driven by climate variability can induce shifts in vegetation patterns, alter species distributions, and affect biodiversity (Zhang et al., 2019). Comprehending the complexities of VPD and its interactions with various environmental factors is essential for sustainable agriculture, ecosystem conservation, and effective fire management.

Since VPD is a function of temperature (Lawrence, 2005), the effects of climate change and the associated increase in the mean global temperature lead to positive trends in the regional and global VPD (Grossiord et al., 2020; IPCC, 2021b). For instance, studies have shown that the VPD has been increasing sharply at a global scale since the year 2000 (Simmons et al., 2010; Willett et al., 2014; Yuan et al., 2019). Spatially explicit VPD records derived from remote sensing data cover only the last ∼50 years and vary in quality, so long-term perspectives of VPD variability are lacking. However, long-term perspectives can help to put recently observed trends of VPD in a long-term context relevant for estimating the significance and robustness of these changes at both local and continental scales. Furthermore, it is essential to investigate the independent physiological effects of VPD on large-scale vegetation dynamics, a topic that is far less explored (Grossiord et al., 2020). So far, the first local reconstruction studies have shown the potential for a long-term perspective on VPD (Liu et al., 2017; Roibu et al., 2022). Churakova Sidorova et al. (2020) have shown that the recent VPD increase does not yet exceed the maximum values reconstructed during the Medieval Warm Anomaly in Siberia. Nevertheless, most studies lack a wider spatial perspective as they only reconstruct VPD time series for a single location. Therefore, it appears necessary for high-resolution VPD reconstructions over extended spatial regions (e.g., the European region).

Tree ring sequences are one of the most used proxy archives for paleoclimatic and environmental reconstructions because of their annual resolution, precise dating, and widespread spatial extent (Leonelli et al., 2017). Moreover, they allow us to create chronologies of thousands of years and give us the possibility to explore the climate through different tree-ring-based climate-sensitive parameters, such as tree ring width, maximum density, and/or stable isotopes (Nagavciuc et al., 2022). Of all physical or chemical parameters to be measured of annual rings, the stable oxygen isotope ratio of tree ring cellulose (δ¹⁸O) has the advantages that the isotopic fractions controlled by physiological processes are reasonably well understood and the statistical relationship with climatic and environmental quantities is rather robust (McCarroll and Loader, 2004; Siegwolf et al., 2022). Combining these advantages with the advantages of annual tree rings, stable oxygen isotopes in tree ring cellulose may be one of the most valuable proxies for the paleoclimatic reconstructions. Nonetheless, there are still questions to be answered about the variability of δ¹⁸O in the arboreal system (Gagen et al., 2022).

The δ¹⁸O in tree ring cellulose is mainly influenced (i) by δ¹⁸O signature of precipitation (through water uptake by the roots from the soil) (McCarroll and Loader, 2004; Siegwolf et al., 2022); (ii) by biochemical fractionation including partial isotopic exchange of cellulose precursors with stem water during cellulose biosynthesis (Barbour, 2007; Roden et al., 2000; Saurer et al., 1997), which is considered to be largely constant at 27±4 ‰ (Sternberg and DeNiro, 1983); and (iii) most importantly by the evaporative ¹⁸O enrichment of leaf or needle water via transpiration of water vapor to the atmosphere (Kahmen et al., 2011; Roden et al., 2000; Saurer et al., 1997). Trees absorb water from the soil, which reflects the δ¹⁸O values of local precipitation. As δ¹⁸O values of precipitation change, this is reflected in the δ¹⁸O values of the trees' cellulose, providing a valuable proxy for past climate conditions. Balting et al. (2021) analyzed the relationship between intra-annual tree ring δ¹⁸O data from the ISONET network and δ¹⁸O in precipitation water and found significant correlations, results which are supported by numerous local-scale studies (Jones et al., 2016; Locosselli et al., 2020). The transpiration process is controlled by leaf-to-air VPD modified by the aperture of stomata controlling the conductance of water vapor (Buckley, 2019; McCarroll and Loader, 2004). The δ¹⁸O values of leaf water are typically enriched in ¹⁸O compared to the plant's parent water because evaporative losses are greater for the lighter ¹⁶O than for ¹⁸O (Roden et al., 2000). However, the δ¹⁸O signature of chloroplast water is not in isotopic equilibrium with leaf water at the actual sites of transpiration (stomata), i.e., the higher the transpiration rates, the lower the rate of enrichment of the chloroplast water, the so-called Péclet effect (Barbour et al., 2004). To a good approximation, tree ring δ¹⁸O from European sites can be considered a combined signal of the largely temperature-dependent δ¹⁸O of the water source (precipitation, soil water) and the evaporative ¹⁸O enrichment of leaf water controlled by leaf-to-air VPD, meaning that tree ring δ¹⁸O can be used as a proxy for variations in VPD (Ferrio and Voltas, 2005; Kahmen et al., 2011).

In order to have a long-term overview of the observed changes and trends in the regional and global VPD, one needs to look back in the past by employing different proxy indicators (e.g., from tree rings). In this respect, the availability of the ISONET network (Balting et al., 2021; ISONET Project Members et al., 2023; Treydte et al., 2007a, b) gave us the opportunity to develop the first gridded summer VPD reconstruction at the European level by using multiple stable oxygen isotope records derived from tree rings. The reconstruction of our VPD dataset, which covers the last 400 years, has been obtained by applying a random forest (RF) regression method (Breiman, 2001). According to Yang et al. (2020), RF has become one of the most successful machine learning algorithms for practical applications over the last 2 decades due to its proven accuracy, stability, speed of processing, and ease of use (Bair et al., 2018; Belgiu and Drăgu, 2016; Maxwell et al., 2018; Qu et al., 2019; Reichstein et al., 2019; Rodriguez-Galiano et al., 2012; Tyralis et al., 2019). The term RF describes a nonlinear and robust technique in which several decision trees are built and aggregated at the end to make predictions or perform reconstructions (Breiman, 2001). The RF approach is increasingly applied in climate and environmental sciences and has been used for the prediction of snow depth (Yang et al., 2020), solar radiation (Prasad et al., 2019), daily ozone (Zhan et al., 2018), and precipitation (Ali et al., 2020); for reconstructions of last millennium North Atlantic Oscillation (NAO; Michel et al., 2020), Atlantic Multidecadal Variability (Michel et al., 2022), El Niño–Southern Oscillation (Delcroix et al., 2022), streamflow since 1485 CE (Li et al., 2019), and vegetation cover during the mid-Holocene and the Last Glacial Maximum (Lindgren et al., 2021); and for calculating tree growth response to earthquakes (Mohr et al., 2021). Although RF models have proved to be a useful method in geosciences, studies on the spatial–temporal reconstruction of climate variables based on δ¹⁸O are relatively rare due to the low availability of δ¹⁸O time series. Here, the RF model is used independently for the reconstruction of each grid point, meaning that we build several one-dimensional climate time series reconstruction models where the respective weights of the δ¹⁸O values of tree ring cellulose vary for each grid point (see Sect. 2). In this regard, recent studies reconstructing climate time series as done here showed that RF generally outperforms most commonly used variants of linear regression model (e.g., principal component regression, partial least squares; Delcroix et al., 2022; Michel et al., 2020, 2022).

The purpose of this study is to apply an RF approach (Breiman, 2001) for reconstructing the European summer (averaged for June–July–August) VPD for the first time on a gridded scale from a proxy network that is based on 26 series of tree ring δ¹⁸O and covers the period 1600–1994 (Balting et al., 2021; ISONET Project Members et al., 2023). The main aim is to present the summer VPD reconstruction dataset for the last 400 years over Europe and provide both a spatial and a temporal long-term perspective on the past summer VPD variability. The VPD reconstruction is analyzed subsequently from a spatiotemporal perspective for northern Europe, central Europe, and the Mediterranean region (Iturbide et al., 2020), as defined and used in the Sixth Assessment Report (AR6) of the IPCC (IPCC, 2021a). The data evaluation in this paper is structured as follows. In Sect. 2 we give a detailed description of the data and methods employed, while the main results are presented in Sect. 3. First, we present the sensitivity of δ¹⁸O for VPD variability and the validations statistics for the reconstruction. Second, we investigate the spatial variability of VPD in comparison with other studies and analyses for selected extreme years for additional validation of our reconstructions. We also present the relationship between summer VPD and large-scale atmospheric circulation in a long-term context. Statistical uncertainties are discussed in Sect. 4, while conclusions and outlooks are presented in Sect. 5. In this way, the past and present VPD conditions are presented and evaluated to assess the statistical significance of observed trends to understand VPD variability on a local, regional, and continental scale.

2 Data and methods

2.1 The stable isotope network

To reconstruct the European summer VPD, we used 26 time series of stable oxygen isotopes in tree ring cellulose (δ¹⁸O) (Fig. 1, Table 1). From the 26 time series used in this study, 21 time series were obtained from the dataset generated by the EU project ISONET (EVK2-CT-2002-00147) (Balting et al., 2021; ISONET Project Members et al., 2023; Treydte et al., 2007a, b). In addition to the ISONET dataset, we added five new δ¹⁸O time series from Bulgaria, Turkey, southwestern Germany, Romania, and Slovenia (Balting et al., 2021; Hafner et al., 2014; Heinrich et al., 2013; Nagavciuc et al., 2019) (Fig. 1).

Figure 1 The site distribution of the δ¹⁸O network composed of 21 ISONET and 5 additional datasets (see Table 1). Colors indicate the past temporal availability of each tree ring dataset: starts in 1600 or earlier (purple), starts between 1601 and 1700 (blue), starts between 1701 and 1800 (orange), and starts between 1801 and 1900 (red). Symbols indicate tree species used for each record: Juniper (diamond), Larix (triangle), Pinus (circle), and Quercus (square).

The detailed measurement methodologies used within the ISONET project and the other sites are described in the description file of the dataset (ISONET Project Members et al., 2023) and also described in more detail in previously published papers (Boettger et al., 2007; Hafner et al., 2014; Heinrich et al., 2013; Nagavciuc et al., 2019; Treydte et al., 2007b, a). At least four dominant trees were selected at each site and two increment cores were taken per tree for the ISONET project; 15 cores of 5 living trees were taken from the site in Turkey (Heinrich et al., 2013), 12 trees were sampled for the Slovenian time series (Hafner et al., 2014), and 1 core per tree was taken from 9 trees in Romania (Nagavciuc et al., 2019). A standard dendrochronological dating method was performed (Fritts, 1976), and subsequent individual growth rings were dissected from the cores. All tree rings from the same year were pooled for most sites prior to cellulose extraction for the ISONET sites (Treydte et al., 2007b, a) and the Romanian site (Nagavciuc et al., 2019). The dissected tree rings from the Slovenian and Turkish sites were measured individually and not pooled (Hafner et al., 2014; Heinrich et al., 2013). For oak, only the latewood was used for the analyses. This approach was assumed to use predominantly climate signals from the current year, as the earlywood of oaks can contain climate information from the previous year (Davies and Loader, 2020; González-González et al., 2015). All δ¹⁸O time series were measured according to standard methodology, and the inter-lab differences are assumed negligible for future analysis. The results are expressed using the conventional δ (delta) notation, i.e., in per mil (‰) relative to the Vienna Standard Mean Ocean Water (VSMOW) (Craig, 1957).

Figure 2 The number of available δ¹⁸O time series within the network composed of 21 ISONET and 5 additional data (see Table 1). The Quercus sites are represented with green color, Pinus sites with yellow, Larix with red, and Juniperus with grey. For site codes, please see Table 1.

The isotope network presented here consists of nine deciduous tree sites (Quercus) and 17 conifer sites (Pinus, Juniper, Larix; see Table 1). The sample sites are well distributed over Europe (Fig. 1). The elevation of the locations varies from 10 m a.s.l. (Woburn, UK) to 2120 m a.s.l. (Pedraforca, Spain). The longest chronologies cover a period from 1600 to 2005. The highest data density (i.e., 26 time series) is available for the period 1900–1994 (Fig. 2). For several sites or regional groups of sites from the ISONET datasets, the data are published within individual studies (Andreu-Hayles et al., 2017; Etien et al., 2008; Haupt et al., 2011; Helama et al., 2014; Hilasvuori et al., 2009; Labuhn et al., 2014, 2016; Rinne et al., 2013; Saurer et al., 2014, 2008, 2012; Vitas, 2008).

Table 1 Characteristics of each sample site used within our study. A total of 21 of the 26 δ¹⁸O records were obtained from the EU project ISONET, and the 5 additional sites are from Bulgaria, Turkey, southwestern Germany, Romania, and Slovenia.

2.2 (Paleo)climate data

Mean air temperature (°C) and relative humidity (%) from the near-surface (i.e., 2 m above ground level) were derived from the 20th Century Reanalysis Project (20CR) version V3 (Slivinski et al., 2019) at a monthly resolution to build the statistical model for VPD reconstruction. The 20CR reanalysis has a temporal resolution of 3 h, 28 different pressure levels, and a resolution of 1° × 1°. We use the ensemble mean derived from an 80-member ensemble. The climate variables are available for the period from 1836 to 2015, and they are provided by NOAA/OAR/ESRL PSL, Boulder, Colorado, USA (https://psl.noaa.gov/data/gridded/data.20thC_ReanV3.html, last access: 16 October 2023). To analyze the large-scale drivers of summer VPD in a long-term context, we make use of the monthly paleoclimate reanalysis (i.e., EKF400v2) (Valler et al., 2022) covering the period 1604 to 2003. This paleoclimate reanalysis dataset was generated using data assimilation techniques for the reconstruction of past climate states using both observation information and model simulations. The EKF400v2 is based on atmospheric-only general circulation model simulations, an observational network of early instrumental temperature and pressure data, documentary evidence, and tree ring width and density proxy records (Valler et al., 2022).

2.3 VPD computation and further pre-processing

The monthly VPD values were calculated based on the equation presented by Barkhordarian et al. (2019) by using the near-surface air temperature (e.g., 2 m above the ground) (T) and the dew point temperature (T_d), both in °C. Since dew point temperature is not available for each grid point of the datasets, we have used T and relative humidity (RH, %) to compute T_d as follows:

$$\begin{array}{}\text{(1)}& T_{\mathrm{d}}={\displaystyle \frac{a_{\mathrm{1}}(\ln \left(\frac{\mathrm{RH}}{\mathrm{100}}\right)+\left(\frac{T{a}_{\mathrm{2}}}{T+a_{\mathrm{1}}}\right))}{a_{\mathrm{2}}-(\ln \left(\frac{\mathrm{RH}}{\mathrm{100}}\right)+\left(\frac{T{a}_{\mathrm{2}}}{T+a_{\mathrm{1}}}\right))}},\end{array}$$

where a₁ and a₂ are defined as a₁= 243.04 and a₂= 17.625. This computation is reliable and is notably used in many climate models (Barkhordarian et al., 2019). We utilize the Clausius–Clapeyron relation by applying a term for the saturation vapor content of the air and a term for the actual vapor pressure to calculate VPD as follows (Barkhordarian et al., 2019; Behrangi et al., 2016; Marengo et al., 2008; Seager et al., 2015):

$$\begin{array}{}\text{(2)}& \mathrm{VPD}=c_{\mathrm{1}}\times e^{\left(\frac{c_{\mathrm{2}}\cdot T}{c_{\mathrm{3}}+T}\right)}-c_{\mathrm{1}}\times e^{\left(\frac{c_{\mathrm{2}}\cdot T_{\mathrm{d}}}{c_{\mathrm{3}}+T_{\mathrm{d}}}\right)},\end{array}$$

where c₁= 0.611 kPa, c₂= 17.5, c₃= 240.978 °C, and VPD is expressed in kPa (see WMO, 2018, for further information). We calculated these quantities for monthly reanalysis averages and computed seasonal VPD averages for summer (JJA June to August) upon which we will focus our analyses in the following.

Here, we use the 20CR reanalysis as it spans a large enough time frame to calibrate statistical reconstruction models, which is at the expense of increased data uncertainties. To ensure that values from 20CR are in line with other historical VPD estimates, we also computed VPD for the ERA5 reanalysis dataset (Hersbach et al., 2020) over the European area (35–75° N, 10° W–40° E), where we found good agreement between both datasets (Supplement Fig. S1).

2.4 Application of random forest

The VPD data derived from 20CRV3 (Slivinski et al., 2019), Sect. 2.3) and the δ¹⁸O network data were used as input. There is a small number of missing data in the δ¹⁸O time series (0.38 % entries in total) that are infilled using an iterative principal component analysis approach following Josse and Husson (2016).

To spatially reconstruct the European summer VPD for the last 400 years we have used the random forest (RF) method (see Supplement Sect. S1 and Fig. S2 for details) (Breiman, 2001). The RF is an ensemble bootstrapping method consisting of aggregating several predictions given by simple, and generally low-efficiency, nonlinear regression trees (Breiman, 2001). Regression trees consist of a recursive convex optimization. First, the predictor is chosen (δ¹⁸O time series here) that best splits the initial summer VPD values in time into two groups (leaves) where the Euclidean distances from group centers are minimized. This is recursively applied where the specific group giving the best split is also considered in the minimization (on top of the δ¹⁸O time series and its values best splitting summer VPD data). The algorithm stops when one group is composed of a threshold number of VPD values. In the case of using regression trees within a RF, this threshold is recommended to be set to 5, which we will do here for all RF models, as its optimization is costly and generally not worth the computation (Breiman, 2001). Regression trees generally overfit the data they are trained with since removing one single predictor (for instance the δ¹⁸O time series used at the first split of the tree) can result in drastically different modeled data (Breiman, 2001). This is the reason why RF was developed: it consists of constructing several regression trees, each based on a randomly drawn sample of the initial set of δ¹⁸O time series. The number of trees was shown to give converging predictions for a wide variety of datasets when set to 256 or higher (Oshiro et al., 2012). Hence, here we used a number of 256 trees for all RF models we have computed. Finally, the number of random predictors (i.e., δ¹⁸O time series) selected is tuned using a 10-fold cross-validation (Michel et al., 2020, 2022) for all RF models of the study in a possible range extending from one to the number of initially selected δ¹⁸O time series (see next paragraph); see Sect. S1 for further details. Once a RF model is trained and tuned from cross-validation, a reconstruction is obtained by browsing all trees from the training using the values of δ¹⁸O time series in the past. The reconstructed summer VPD value is obtained as the average of values.

2.5 General procedure of the European VPD reconstruction

The RF algorithm as used here is not directly adapted for spatial reconstruction. We thus computed several RF models tuned with cross-validation, each specifically computed for the reconstruction of a single grid point. For climate time series reconstruction as performed here, RF was shown to most generally overclass several linear regression methods, including principal component regression, which has been extensively used in the climate reconstruction literature (e.g., Ortega et al., 2015; Wang et al., 2017). For each grid point of the observations, sensitivity tests are performed prior to the reconstruction, where only δ¹⁸O data significantly correlated at the 95 % confidence level when the observed VPDs are considered. This makes the reconstruction of each grid point from the summer VPD dataset unique as only relevant δ¹⁸O time series are selected for each. These tests are carried out with a correction of the degrees of freedom using first-order autocorrelation coefficients estimated from both time series, as in Michel et al. (2020).

In our study, we focus on the continental area of the European region (35.5 to 74.5° N, 9.5° W–39.5° E). The learning is set from 1900 to 1994, which constitutes the longest period covered by all tree ring data from Table 1 starting from 1900 (Fig. 2). Because δ¹⁸O time series cover different periods in the past, we use a nesting approach to produce the reconstruction. Starting from the longest time frame of reconstruction (i.e., 1600–1994), a new gridded reconstruction is computed as a new δ¹⁸O time series becomes available (see Table 1). At the end of the calculations, all imbricated time series for a given grid point are put one after the other to obtain the final reconstruction (Michel et al., 2022), which is done for all cells from the grid.

Each of the reconstructions for each grid point and nested time frame is evaluated by splitting up the data 30 times through time into training and testing samples of 80 % and 20 % of the size of the initial learning period (i.e., 1900–1994), respectively. For each training sample, a RF model is built and tuned, and a reconstruction of summer VPD for the testing period is obtained using corresponding values from δ¹⁸O time series only. The evaluation over testing samples is made using the coefficient of efficiency score (S_{CE}, Nash and Sutcliffe, 1970). This score is defined from −∞ to 1, where positive values indicate better reconstruction skill than what the empirical average of the training set would give. The reconstruction for a grid point over a given time frame of the nested reconstruction is thus considered robust and reliable only if the average of the 30 S_{CE} is significantly positive at the 95 % confidence level through the sample, based on a one-tailed Student's t test. Hence, our final reconstruction extending back to 1600 is composed of gaps in space and time, where filled data are only those where robust reconstruction skills as described here were found. Once the RF model is evaluated through the 30 training and testing splits, a final reconstruction is computed as the median of the 30 individual reconstructions obtained from training and testing splits. Therefore, statistical uncertainties of the gridded reconstruction are provided as the 10 %–90 % range of these 30 individual reconstructions.

In addition, standard errors are computed from reconstructed time series with respect to summer VPD observations over the period 1900–1994. For each grid point and each time frame of the nested reconstruction, we computed the uncertainties as an envelope of twice the standard errors on each side of the reconstruction as in, e.g., Ortega et al. (2015).

2.6 Reconstruction analyses

In order to further test the feasibility of our VPD reconstruction, in the final step of the analysis we compare the regional averages and running averages (30 years) of the reconstructed VPD with temperature (Luterbacher et al., 2004), precipitation (Pauling et al., 2006), and Palmer Drought Severity Index (PDSI) reconstructions (Cook et al., 2015) for the three European regions, as defined by IPCC Sixth Assessment Report (AR6) (IPCC, 2021a; Iturbide et al., 2020): northern Europe, central Europe, and the Mediterranean. Furthermore, we show maps of the European summer VPD for the most extreme positive years (1868, 1707, 1835) and most extreme negative years (1785, 1742, 1747) to investigate the spatial variability in concordance with past historical data (Brázdil et al., 2013; Brooks and Glasspoole, 1922; Cook et al., 2015; Glaser, 2008; Ionita et al., 2021; Marusek, 2010; Pauling et al., 2006; Trigo et al., 2009). Prior to this mapping, all VPD grid cells are centered and standardized (z transformation) to present z anomalies for each grid cell.

3 Results and discussion

3.1 Vapor pressure deficit – δ¹⁸O relationship

The network of δ¹⁸O time series is significantly correlated with the European Summer VPD over the observational period (i.e., 1900–1994). In Fig. 3, we show the correlations of the individual δ¹⁸O time series with the respective time series of the grid cell from the calculated VPD of 20CRV3 (Slivinski et al., 2019). High values of significant correlations are reached between the local summer VPD and the δ¹⁸O time series (Fig. 3). We found that 23 of the 26 sites show significant correlations with local summer VPD (i.e., from the closest grid point) at the 95 % confidence level, with most of these sites being located in France, Spain, Germany, Scandinavia, and Great Britain. One time series from Turkey, one from Finland, and one from Bulgaria do not have a significant correlation with local summer VPD (Fig. 3). We note, however, that tree rings from southeastern Europe (e.g., Turkey, Bulgaria) may be more sensitive to spring VPD rather than summer VPD as demonstrated by a former study (Heinrich et al., 2013). In addition, we notice that there is no influence of the tree species used and no influence of the altitude of sampled trees on the correlations found (Fig. 3, Table 1, Balting et al., 2021). Based on these results and the relationship between δ¹⁸O and VPD variability depicted above, we can argue that δ¹⁸O in tree ring cellulose can be used as a reliable proxy to perform the reconstruction of VPD for the summer months (JJA).

Figure 3 Correlation of the δ¹⁸O time series with local summer VPD for the period 1900 to 1994. For each tree ring dataset, the correlation is computed with the time series from the closest grid point. The significance of correlations is calculated using a Student's t test for correlation with corrected degrees of freedom from time series as in McCarthy et al. (2015) and Michel et al. (2020). Significant correlations at the 95 % confidence level are highlighted in red (p < 0.05).

3.2 Validation statistics

Since we used a nesting approach in our study (Sect. 2.5), we made 13 successions of RF models over space in total (Sect. 2.5), each starting from time steps where at least one time series from the tree ring network becomes available (see Sect. 2.5). The S_{CE} scores (see Sect. 2.5) are shown in Fig. 4 for four selected time steps (1600–1699, 1700–1799, 1800–1899, and 1900–1999).

Figure 4 The median Nash–Sutcliffe S_{CE} for four different time slices. The colored areas show the regions where the model setup has a suitable quality, which was tested with a one-sided Student's t test for the average of 30 S_{CE} scores obtained for testing samples (p < 0.05, Sect. 2.5). For the sake of graphical representation, for each time slice plotted here, the median S_{CE} scores were averaged over the corresponding nested time frames. Locations of available records for each time slice are given by green crosses (see Table 1).

For the first time step 1600 and older, S_{CE} scores indicate a satisfactory quality (Sect. 2.5) for the reconstruction of summer VPD for northeastern Spain, Italy, Greece, France, Germany, and large parts of Scandinavia (Fig. 4a). These regions coincide with the locations of the 11 available time series for this time frame. For 1700 or older, when our δ¹⁸O network has increased to 16 time series, the validation scores became higher, and the spatial coverage of grid points with a significantly robust reconstruction at the 95 % confidence level now expands to some parts of eastern Europe and the whole Scandinavian region (Fig. 4b). A further improvement of the S_{CE} values can be observed for the time span until 1800 (Fig. 4c), whereby the most noticeable feature is that Great Britain is now largely covered by grid points with a significantly robust reconstruction at the 95 % confidence level. This improvement is due to the fact that another time series from Scotland has been included in the analysis for this time slice. Finally, for the last time step (i.e., 1900–1994), for which all 26 δ¹⁸O time series are now available, an improvement in the spatial coverage of the reconstruction is found, in particular, a large gap in eastern Europe is now robustly covered by our reconstruction due to an additional time series from Romania, which became available from 1900 onwards (Fig. 4d). It can also be noted that adjacent regions of Europe, such as parts of Turkey, also have a suitable quality for reconstructing past VPD variability changes since a tree ring series is available in this region. The overall spatial strength of the network is in southern, western, and northern Europe, whereas eastern Europe can only be partially covered. Therefore, only those grid cells with a satisfactory reconstruction performance (where the mean S_{CE} scores are significantly positive at the 95 % confidence level) will be included in the following analyses. In addition to the validation scores, we also computed and provided statistical uncertainties (Fig. 5) in order to compute error bar estimates of the reconstruction (see Sect. 2.3). Overall, the validation statistics passed the conventional verification tests, and the obtained statistical test results indicate that the reconstruction model was statistically sound and can be used for the reconstruction of summer VPD values at the European scale.

Figure 5 Regression uncertainties for each grid cell. Uncertainties are calculated as the time-averaged 10 %–90 % range of individual reconstructions from each grid point (for the longest period with reconstruction skills, see Fig. 4).

3.3 Temporal variability of the European summer VPD reconstruction

The VPD reconstruction for the Mediterranean region shows the largest VPD values and higher variability compared to the other two regions (Fig. 6a, b). It is noticeable that at the beginning of the reconstruction period (1610–1650), the VPD values decrease ($m=-\mathrm{0.0017}$ kPa yr⁻¹; 0.05 < p value < 0.01, where m is the linear regression slope), which is followed by a VPD increase leading to an almost constant VPD level until ∼ 1715. After this point, the VPD decreased steadily until 1761 when it remained at a similar level until 1773. This is followed by a short increase until ∼ 1790, after which VPD decreases again. The minimum VPD value is also reached during this period (VPD = 0.96 kPa for the year 1814). Around 1830, the VPD values started to increase again (m= 0.0003 kPa yr⁻¹; 0.05 < p value < 0.01). Furthermore, the reconstructed VPD time series shows similar variations to those for central Europe, but they are most pronounced in the last 150 years. Low periods of variability are shown, for example, between the end of the 19th century and the beginning of the 20th century and for the 1960s and 1970s. In total, the 3 years with the lowest VPD z anomalies in this region are 1735, 1814, and 1876, whereas the years with the highest VPD are 1682, 1686, and 1945.

Figure 6 Temporal variability of the VPD reconstruction. Maps (a, c, e) give the location of regions investigated: Mediterranean (a, green), central Europe (c, blue), and northern Europe (f, orange), following the IPCC AR6 definition (Iturbide et al., 2020). Area-averaged reconstructed time series (b, d, f) for (top to bottom): VPD (this study), surface temperature (Luterbacher et al., 2004), precipitation  (Pauling et al., 2006), and PDSI (Cook et al., 2015) for the period 1600–1994, for the same three regions. For VPD time series in (b), (d), and (f), grey-shaded areas give the range of spatially averaged VPD uncertainties from Fig. 5 at the annual timescale, where values are spatiotemporally filtered based on reconstruction skills (Fig. 4).

The reconstructed VPD for Central Europe (Fig. 6c, d) shows that VPD increases from 1600 up to ∼ 1700 (m=0.0007 kPa yr⁻¹; p value < 0.01). The maximum of the 30-year moving time series is reached in 1697 (VPD = 0.6398 kPa > 2.4 × σ, where σ is the standard deviation of the time series). This increase is followed by a downward trend of VPD in central Europe that ends in 1743 (m = −0.0014 kPa yr⁻¹; 0.01 < p value < 0.05), when the minimum VPD value was recorded. The period 1740 to 1760 is characterized by very low VPD values and the lowest 30-year rolling VPD. From this time on, the rolling average VPD is characterized by a rising trend (m= 0.0002 kPa yr⁻¹; p value < 0.01), but low values are reached between 1790 and 1830. Furthermore, the VPD reconstruction is characterized by a significant 60- to 80-year variability, where the lowest values are reached during the periods 1890 to 1920 and 1960 to 1980 (Fig. S3). In total, the 3 years with the lowest VPD in this region are 1814, 1838, and 1980, whereas the years with the highest VPD are 1687, 1707, and 1835.

The average VPD reconstruction for northern Europe (Fig. 6e, f) shows significant differences in variability and in the long-term mean. Northern Europe shows the lowest VPD of all three regions, which is due to the comparably low temperatures and high humidity in this area. The increase in VPD between 1620 and 1660 is significant (m = 0.0017 kPa yr⁻¹; p value < 0.01), with the highest VPD value of the 400 years being reached in 1652 (VPD = 0.4619 kPa > 2.7 × σ). However, we note that the increase during 1675–1695 is not as long-lasting as the former increase from 1620 to 1660. The VPD values start to decrease again from the year 1700. The subsequent drop lasts until the 1720 to 1730 period, after which the time series shows further low-frequency variability and an upward tendency in the 30-year rolling average that lasts until the end of the 20th century (m= 0.00002 kPa yr⁻¹; p value < 0.01). The 3 years with the lowest VPD in this region are 1632, 1674, and 1802, whereas the years with the highest VPD are 1652, 1959, and 1973.

3.4 Validation of the VPD reconstruction by comparison with other datasets

In order to test the feasibility of our VPD reconstruction, in the following two sub-sections we compare the spatiotemporal variability of our VPD reconstruction with other datasets available for the same time span. Time series from Fig. 6b, d, and f show some interesting similarities between our summer VPD reconstruction and former ones based on different climatic variables (e.g., precipitation, temperature, drought indicators) (Cook et al., 2015; Luterbacher et al., 2004; Pauling et al., 2006). For the overlap period of these reconstructions (1600–1994), we find consistent positive correlations with area-averaged reconstructed summer temperatures at 2 m (T2m, Luterbacher et al., 2004). The obtained correlations vary between r= 0.36 (p < 0.01) for the Mediterranean region, r= 0.31 (p < 0.01) for central Europe, and r= 0.38 (p < 0.01) for northern Europe (Fig. 7). When averaged over the same three regions, an opposite (negative) relationship is found between our VPD reconstruction and a summer precipitation reconstruction (Pauling et al., 2006), where correlations for spatially averaged time series are $r=-\mathrm{0.26}$ (p < 0.01) for the Mediterranean region, $r=-\mathrm{0.35}$ (p < 0.01) for central Europe, and $r=-\mathrm{0.39}$ (p < 0.01) for northern Europe (Fig. 7). The same results are found when we use a drought indicator, namely the Palmer Drought Severity Index (PDSI). In the case of PDSI, we find significant negative correlations for the three analyzed regions, namely $r=-\mathrm{0.39}$ (p < 0.01) for the Mediterranean region, $r=-\mathrm{0.24}$ (p < 0.01) for central Europe, and $r=-\mathrm{0.34}$ (p < 0.01) for northern Europe (Fig. 7). Interestingly, the Mediterranean region has highest correlations on average for T2m, whereas it has the lowest correlations for precipitation and PDSI.

Figure 7 Comparison with former climate field reconstructions. (a–c) Correlations between the VPD reconstruction from this study and other climate reconstructions: (a) temperature at 2 m (T2m) (Luterbacher et al., 2004), (b) precipitation (precip) (Pauling et al., 2006), and (c) Palmer Drought Severity Index (PDSI) (Cook et al., 2015) for the period 1600–1994. Black crosses indicate no significance at the 95 % confidence level. (d–f) Distribution of correlations from (a–c) for the same three variables, respectively, for the three regions investigated, i.e., Mediterranean (green), central Europe (blue), northern Europe (red), and all regions (purple dashed line), following IPCC AR6 (Iturbide et al., 2020; Fig. 5).

The spatial correlations with other (paleo)climate reconstructions (the temperature at 2 m (T2m; Luterbacher et al., 2004), precipitation (Pauling et al., 2006), and PDSI (Cook et al., 2015)) are given in Fig. 7a–c. It appears that for all three variables, significant correlations are found, with a focus on western and central Europe, which are the areas with the highest tree ring data coverage (Fig. 1, Table 1) and that show the highest S_{CE} validation scores (Fig. 4). Despite the relatively low correlation coefficients (lower or higher than ±0.2), the findings present a remarkable achievement considering the vast geographical scope (encompassing the entirety of Europe) and the extensive temporal range (spanning over 400 years) of the correlated datasets. Thus, the obtained significant correlations mean that there is a strong statistical relationship between the analyzed datasets, while the modest correlations suggest that they do not fully capture the complex relationship between these datasets. Therefore, the significant correlations found for most regions highlight the important accordance and consistency between the summer VPD reconstruction and former summer climate reconstructions (Cook et al., 2015; Luterbacher et al., 2004; Pauling et al., 2006).

3.5 Comparison of the temporal variability with other reconstructions

In this section, we compare the temporal variability of VPD with existing reconstructions of summer temperature (Luterbacher et al., 2004), precipitation (Pauling et al., 2006), and drought (PDSI) reconstruction (Cook et al., 2015) in Figs. 6, 7, and 8.

Figure 8 (a) Stability maps between the gridded VPD reconstruction from this study and the reconstructed time series of VPD over southern Fennoscandia (Treydte et al., 2023). (b) The temporal evolution of the VPD time series averaged over the black box in (a) and the VPD time series over southern Fennoscandia (Treydte et al., 2023). Panel (c) is the same as (a) but for western Europe, panel (d) is the same as (b) but for western Europe, panel (e) is the same as (a) but for southeastern Europe, panel (f) is the same as (b) but for western southeastern Europe, panel (g) is the same as (a) but for the Alpine region, and panel (h) is the same as (b) but for the Alpine region. In (a), (c), (e), and (g) the regions where the correlation is stable, positive, and significant for at least 80 % of the 31-year moving window correlations are shaded with dark red (95 %), red (90 %), orange (85 %), and yellow (80 %). The corresponding regions where the correlation is stable but negative are shaded with dark blue (95 %), blue (90 %), green (85 %), and light green (80 %).

For northern Europe, the characteristics of VPD until 1700 show little resemblance with the reconstructions of precipitation and temperature (Fig. 6e, f). However, the PDSI indicates dry conditions for the time of the maximum VPD in 1652 (Fig. 6d). Furthermore, the decrease of VPD between 1700 and 1730 is shown as a wet period in the drought reconstruction, which could be a possible explanation for these low values. From 1730 to 1800, the observed strong temperature increase is only partly visible in VPD. A possible reason for this could be the concurrent increase in summer precipitation. Since the temperature and precipitation changes in this region have similar trends, it is difficult to detect the characteristics of the two variables in the VPD time series.

For central Europe, the increase of the VPD until 1700 cannot be explained by the temperature or precipitation reconstruction (Fig. 6d). The 30-year running means of both time series show diverging trends. Nevertheless, the reconstruction of the PDSI also shows a trend towards drier conditions from 1670 to 1690 (Fig. 6d). However, before the end of the 17th century, the data indicates a wetting trend. Afterwards, the periods of high and low temperature and precipitation match well with the corresponding periods of high and low VPD. For example, a simultaneous drop in temperature and precipitation between 1730 and 1745 is also present in the VPD time series. From 1800 onward, the VPD closely follows the pattern of temperature (Fig. S3). However, the variability of the VPD time series is also influenced by the precipitation variability. The precipitation shows strong interannual variability from the second half of the 19th century, but the magnitude is weaker than in the VPD or the temperature time series. The temperature, precipitation, and VPD reconstructions show similar patterns, it gets warmer, there is less precipitation, and the VPD increases in summer. The PDSI, on the other hand, shows a wetting trend from 1800 onwards. We suggest that this could be based on the accumulation processes of precipitation in the preceding months, which is represented by the PDSI (Alley, 1984). This could compensate for the increasing temperature and the decreasing precipitation in summer.

For the Mediterranean area, the summer of 1884 is an extreme event in the temperature time series, as it shows the second coldest temperatures for this region in the last 400 years (Fig. 6b). The variability of periods with high and low temperatures are almost identical to the corresponding periods of the VPD time series. During the period 1600–1730, the summer VPD reconstruction presents a high variability, with similar variability being observed for the precipitation reconstruction. After which the summer VPD reconstruction presents a slowly decreasing trend, similar to the temperature reconstruction. Since the decrease in precipitation is almost linear in the 30-year mean from 1850, we suggest that comparatively the variability of temperatures can be clearly observed in the VPD time series. A similar situation is also described by the drought reconstruction for the Mediterranean area (Fig. 6b). Thus, the reconstructed PDSI of Cook et al. (2015) shows wetter conditions towards the middle of the 18th century, which is complemented by the low VPD values for this period. Like the VPD and the temperature reconstructions, the PDSI shows drier conditions after the 1770s, intensifying after the onset around 1790. Furthermore, the year 1814 is also the year with the wettest conditions in the PDSI time series. After ∼ 1850, the anti-correlation between the temperature and VPD and the PDSI is again evident, meaning that variability is also visible.

In line with the recently published study of Treydte et al. (2023) dealing with the intensification of the atmospheric drying over Europe, in the following we also compare the temporal variability of our VPD reconstruction with the data from the aforementioned paper. The comparison is made only for particular regions since no gridded data are available for the Treydte et al. (2023) data. In the following, we analyze the relationship between the VPD reconstruction of Treydte et al. (2023) for different regions across Europe, identified based on principal component analysis (PCA), namely southern Fennoscandia, western Europe, eastern central Europe, and the Alps and Pyrenees (see Fig. 1 in Treydte et al., 2023). First, we compute the stability maps (Ionita et al., 2008; Nagavciuc et al., 2019) between the aforementioned time series and our gridded VPD reconstruction (Fig. 8, left column). As can be inferred from Fig. 8, the correlation between the time series for the four regions across Europe and our gridded summer VPD reconstruction is stable and significant throughout the analyzed period (i.e., 1600–1994). We chose the stability maps, instead of simple correlations, because we wanted to extract only the information for which the relationship between the analyzed variable stays stable in time. Based on the stability maps, we defined four indices by averaging the reconstructed VPD over the grid points where the correlation was stable and above the 90 % significance level (black squares in Fig. 8, right column). These four indices correspond to the same regions as the ones used in Treydte et al. (2023), with the only exception that ours are averaged over the gridded VPD reconstruction, while the others are based on a clustering algorithm. For southern Fennoscandia (Fennoscandia in our graph), the correlation between our reconstruction and the one from Treydte et al. (2023) is r = 0.75 (p < 0.001), indicating a very good agreement between the two time series. In addition, the interannual and the decadal variability of the two-time series is in good agreement. For western Europe the correlation between the two reconstructions is r = 0.68 (p < 0.001), also indicating a good agreement between the time series over this region.

In the case of southeastern Europe and the Alpine region, the correlation coefficients between the two VPD reconstructions are smaller compared to southern Fennoscandia and western Europe (r = 0.52, p < 0.01 for southeastern Europe and r = 0.51, p < 0.01 for the Alpine region). Nevertheless, the interannual and decadal variability is coherent in both reconstructions (e.g., prolonged periods with low or high VPD are similar in both time series). In general, the reconstruction of Treydte et al. (2023) has higher VPD values than our reconstruction, but as previously mentioned the interannual and the decadal variability is similar for both reconstructions. For example, both reconstructions indicate an increase in the VPD at the beginning of the 20th century over all analyzed regions and a period of relatively low VPD values starting in the 1970s (Fig. 8, right column). Some of the discrepancies, especially over the southeastern part of Europe and the Alpine region, might arise from the fact that our reconstruction is based on a network of 26 isotope records, while the one from Treydte et al. (2023) is based on 45 isotope records. Moreover, our indices are defined by averaging the VPD over the gridded reconstruction covering a larger region than the ones in Treydte et al. (2023).

3.6 Large-scale drivers of summer VPD over the last 400 years

Although for other paleoclimate reconstructions (e.g., precipitation, temperature, drought) a clear link has been established between the large-scale circulation and the variability of the respective paleoclimate variables (Nagavciuc et al., 2019; Casty et al., 2005; Jones et al., 2009; Steiger et al., 2018; Cook et al., 2022, among others), in the case of VPD not many studies are available that deal with analyzing the influence of large-scale circulation on the variability and trends of VPD. For example, in their study, Treydte et al. (2023) have inferred that the variability of VPD is influenced by the position of the jet stream and the phase of the NAO, but no analysis was actually performed to show that such a relationship actually exists. Therefore, in this study, we go a step forward by using our gridded VPD reconstruction together with a long-term paleo-reanalysis dataset (EKF400v2; Valler et al., 2022), and we analyze the couple modes of variability between summer VPD and the larger-scale atmospheric circulation (i.e., the geopotential height at 500 mb) over the last 400 years. The large-scale atmospheric circulation can significantly impact vapor pressure deficit (VPD) through its role in controlling humidity and temperature patterns. Persistent high-pressure systems have long been acknowledged for their significant role in generating heat waves, soil moisture deficits, and droughts at both global and regional scales (Bakke et al., 2023; Ionita et al., 2022; Kingston et al., 2015; Schubert et al., 2014). For instance, in regions affected by descending air masses (high-pressure systems), such as subtropical zones, VPD tends to be higher due to the associated warming and drying of the air. Conversely, in areas influenced by ascending air masses (low-pressure systems), such as the Intertropical Convergence Zone, increased moisture content often leads to lower VPD. Thus, the relationship between atmospheric circulation and VPD is crucial in understanding and predicting local water availability and plant stress. To effectively examine the connection between large-scale summer atmospheric circulation and summer VPD variability in Europe over the last 400 years, we employ a canonical correlation analysis (CCA) to compute the coupled modes of variability between summer geopotential height at 500 mb (Z₅₀₀) and summer VPD. The initial 10 empirical orthogonal functions (EOFs) of Z₅₀₀ during summer (JJA) account for approximately 85 % of the total variance, while the initial 10 EOFs of reconstructed summer VPD capture around 81 % of the total variance.

The first coupled mode of variability between summer Z₅₀₀ and summer VPD reveals that an elongated high-pressure system extending from the central North Atlantic towards Fennoscandia, flanked by a cyclonic circulation over the southeastern part of Europe and over Greenland (Fig. 9a), corresponds to high VPD values (i.e., atmospheric dryness) over the western and central part of Europe, as well as over Fennoscandia, while low VPD values prevail over the southeastern part of Europe (Fig. 9c). There is a clear dipole-like variability in the VPD field between the central and western parts of Europe including Fennoscandia and the southeastern part of Europe. The anomalies in the VPD field follow the ones in the Z₅₀₀ field, meaning that positive VPD values are associated with an anticyclonic circulation, while negative VPD values are associated with cyclonic circulation. The interannual variability of the normalized temporal components of the first CCA pairs for Z₅₀₀ and VPD is shown in Fig. 9e. The two time series are significantly correlated (r = 0.84, p < 0.001). The significant positive correlation between the standardized amplitudes corresponding to CCA1 Z₅₀₀ and VPD CCA1 (Fig. 9e) confirms the significant relationship between the summer VPD variability at the European scale and the prevailing large-scale atmospheric circulation. The two canonical time series present strong interannual variability over the period 1604–1994, with the highest amplitudes (towards higher pressure and atmospheric dryness) over the latter half of the 20th century. The increase in the summer Z₅₀₀ in the second half of the 20th century is in line with previous studies (Bakke et al., 2023), which have shown that the geopotential height at 500 mb has increased significantly over large regions covering the North Sea and the western part of Europe form 1979 onwards.

Figure 9 The first coupled mode of variability (CCA1 – left column) and second coupled mode of variability (CCA2 – right column) between summer Z₅₀₀ (a, b) and summer VPD (c, d). (e) The standardized amplitudes corresponding to CCA1 Z₅₀₀ (black line) and CCA1 VPD (red line). (f) The standardized amplitudes corresponding to CCA2 Z₅₀₀ (black line) and CCA2 VPD (red line).

The second coupled mode of variability between summer Z₅₀₀ and summer VPD reveals that a high-pressure system centered over the southern and central part of Europe and a low-pressure system over the northern part of Europe (Fig. 9b) corresponds to high VPD values (i.e., atmospheric dryness) over the central and southern part of Europe, while low VPD values prevail over Fennoscandia and the British Isles (Fig. 9d). As in the case of CCA1, there is a clear dipole-like variability in the VPD field, but this time it is between the southern and northern parts of Europe. The interannual variability of the normalized temporal components of the first CCA pairs for Z₅₀₀ and VPD is shown in Fig. 9f.

The two time series are significantly correlated (r = 0.72, p < 0.001). The two canonical time series also present strong interannual variability over the period 1604–1994. Two prolonged periods of low VPD (low Z₅₀₀) can be observed over the southern part of Europe, one around the 1750s and the other one after the 1950s, while positive VPD values prevailed at the beginning of the 18th century and the beginning of the 20th century. The NW–SE dipole-like mode (CCA1) and the N–S mode (CCA) in the VPD field captured by the CCA analysis follow the patterns of drought variability over Europe over the observational records (Ionita et al., 2015; Kingston et al., 2015; Nagavciuc et al., 2022).

3.7 Extreme years in the VPD reconstruction and comparison of past or historical data

To illustrate the spatial variability of our VPD reconstruction, we selected three of the driest summers (1842, 1868, and 1835) and three of the wettest summers (1604, 1602, and 1747) when considering time frames with significantly positive S_{CE} scores at the 95 % confidence level for each grid point (Fig. 4). We compared their spatial variability (Fig. 10) with respect to other available reconstructions or documentary pieces of evidence at the European level to once again validate the occurrence of the extreme events in our reconstruction. Since VPD variability and mean state differ strongly for different latitudes, these six extreme years were selected using spatially averaged normalized VPD values (i.e., z anomalies). The year 1747 is part of a relatively wet period that corresponds to strongly decreasing PDSI (Cook et al., 2015; Ionita et al., 2021) for Europe on average, notably for the northernmost regions (Fig. 6f) as also suggested by Fig. 10f. Contrarily, the summer of 1835 is characterized by positive anomalies over southern, western, and central Europe, with the highest VPD anomalies over France, which is in concordance with the extremely hot summer conditions in France during that summer (Marusek, 2010). During the summer of 1835, 41 hot days and 11 very hot days were recorded in Paris, France (Marusek, 2010). Heat and dryness conditions were also registered in Belgium, England, and Switzerland in 1835 (Marusek, 2010), which is in concordance with our results (Fig. 10c). Summer 1868 is characterized by positive anomalies over the western, central eastern, and northern parts of Europe (Fig. 6a). Documentary data indicate a hot and dry summer in England (e.g., on 21 July 1835 a temperature of 34.1 °C was recorded in London) (Marusek, 2010).

Figure 10 Spatial pattern of VPD anomalies for selected years. Values of each grid cell are VPD z anomalies (standardized and centered) with respect to the whole period of study (1600–1994). The chosen years were selected as they are the three with the highest VPD z anomalies (a–c, a is the highest) and the three with lowest VPD z anomalies (d–f, d is the lowest) for preindustrial data (1600–1899). Time series only include years where S_{CE} scores are significantly positive at the 95 % confidence level (see Fig. 4).

Interestingly, none of the three cases with high VPD anomalies (Fig. 10a–c) overall shows spatially consistent positive values over the continent compared to low VPD cases (Fig. 10d–f). Strong positive anomalies appeared over the western and northern parts of Europe in 1868 (Fig. 10a) and western to central Europe in 1707 and 1835 (Fig. 10b and c). In all three cases, these large positive anomalies are accompanied by negative anomalies in other regions (in general northeastern and southeastern Europe, Fig. 10) and form a significant contrast over the continent. Oppositely, the three wettest years happened mainly as a large and spatially consistent set of negative VPD z anomalies, where only scattered positive values occurred in a few regions (Fig. 10d–f).

4 Limitations of the reconstruction

The δ¹⁸O network used in our study (Table 1, Fig. 1) is characterized by specific limitations that inevitably influence the quality of the reconstruction presented here. It is likely that the reconstruction may provide some level of sample sites over-representation in central and western Europe compared to the ones from southeastern, eastern, and northern Europe as shown in Fig. 1. These spatial distribution characteristics of sites can also be seen in the validation statistics, where the regions with a good sample density show higher validation scores (Sect. 3.2). Therefore, further time series from uncovered regions may help in improving the quality and the spatial extent of our reconstruction. Besides the δ¹⁸O climate signal, the observational or reanalysis data of VPD are given by the ensemble mean of 20CRV3 (Slivinski et al., 2019) for the period 1900 to 1994. Even though the quality and quantity of instrumental data available during this period are comparatively good, the ensemble mean can only represent the variability and diversity of the reanalysis with 80 ensemble members to a limited extent.

In our study, the reconstruction is based on the application of the RF algorithm calibrated and evaluated over the period 1900 to 1994. The model is therefore trained to represent exactly this period. Thus, when this model is applied to the previous years, we assume stationarity between the VPD variability and the proxy records, as is generally the case for climate reconstructions (Cook et al., 2015). Since the climate system is not stationary, the assumption of stationarity must be included as a potential source of error. However, the RF approach can represent nonlinearities, which is not possible with the classical approaches, for example, principal component regression (Cook et al., 2015). Therefore, we compared our reconstruction with other available observational, historical, and proxy-based reconstructed data in order to assess errors. While our δ¹⁸O VPD reconstruction provides valuable insights into past climate conditions, it is important to acknowledge the existence of other sources of paleoclimate reconstructions, which could potentially complement or refine our findings and represent a great potential for future studies.

Another potential limitation of our study can be related to the fact that we only used the δ¹⁸O record to reconstruct the summer VPD, and we did not take into account the potential of δ¹³C as a candidate for reconstructing VPD. This was mainly due to the fact that at the time when we wrote this paper, the δ¹³C network was not publicly available, and at this stage it is beyond the scope of our paper to make a δ¹³C-based reconstruction. Moreover, it seems that δ¹⁸O is preferred as a candidate for long-term paleo-reconstructions, especially over Europe (Freund et al., 2023; Nagavciuc et al., 2022; Treydte et al., 2023). Nevertheless, we acknowledge that for future studies one needs to test the potential of δ¹³C as a proxy indicator for paleoclimate reconstructions.

Finally, even if the nested reconstruction approach is used by default for reconstructions (Cook et al., 2007; Freund et al., 2019; Luterbacher et al., 2004; Pauling et al., 2006) to cover the longest possible time range, it is important to check the quality of the model for the respective time range. Therefore, we recommend considering the validation scores for the considered time periods that are also uploaded to the repository.

By acknowledging these limitations, we can better understand the uncertainties associated with paleoclimate reconstructions and avoid over-interpreting the data. This will help us to develop more robust and reliable paleoclimate records in the future.

5 Conclusion and outlook

Here, we present the first gridded reconstruction of the European summer VPD over the past 400 years, a dataset that was produced to cover a large knowledge gap about past European VPD variability. The validation statistics performed to obtain our summer VPD reconstruction passed all of the conventional verification tests (e.g., Nash–Sutcliffe model efficiency coefficient, S_{CE}), and the obtained statistical test results indicate that the reconstruction model was statistically robust. Thus, we can present the first gridded reconstructions over 400 years of the summer VPD at the European level. Our spatial reconstruction of summer VPD conditions over Europe based on tree ring δ¹⁸O records provides evidence of historical changes and a unique long-term context for recent trends and variability of the summer VPD. Moreover, our reconstruction helps to analyze the spatiotemporal variability of summer VPD in a long-term context emphasizing the local and regional impacts of the current climate change.

The obtained results indicate that the past variability of VPD is different for the three regions we investigated, namely northern Europe, central Europe, and the Mediterranean. Over the Mediterranean region much higher summer VPD values were recorded compared to the other two regions for the analyzed period. The lowest summer VPD values were found for northern Europe, which can be explained by lower temperatures and higher humidity in this area. Similar results were found for the variability of the summer VPD.

The comparison of our summer VPD reconstruction with other climate datasets reveals consistent positive correlations with summer temperatures but negative correlations with summer precipitation and a drought indicator in different European regions, demonstrating agreement with previous climate reconstructions, especially in western and central Europe over a 400-year period. Our reconstruction also shows statistically significant correlations with a recently published regional VPD reconstruction (Treydte et al., 2023). The correlations are higher over the northern and western parts of Europe and smaller over the southern part of Europe and the Alpine regions. Some of the discrepancies, especially over southern Europe and the Alpine region, might arise from the fact that our reconstruction is based on a network of 26 isotope records, while the one from Treydte et al. (2023) is based on 45 isotope records, and their reconstruction has a higher number of records, especially over the mountain regions. Some differences might appear also from the different techniques used to reconstruct the summer VPD, namely random forest (this study) versus fuzzy cluster analysis (Treydte et al., 2023). Moreover, we reconstructed the VPD on a gridded field covering all of Europe and made these 400 years of gridded VPD available to the general public for further analysis, while the study of Treydte et al. (2023) reconstructed clusters (i.e., time series) of VPD over different regions in Europe.

This paper also provides new insights into the relationship between summer VPD and large-scale atmospheric circulation, a relationship that has been neglected so far or at least not treated into a direct manner. Here we show that summer VPD has two preferred modes of variability, namely a NW–SE dipole-like mode and a N–S dipole-like mode. In general, high VPD values are associated with an anticyclonic circulation, while low VPD values are associated with cyclonic-like circulations. Since most of the climate simulations from the Coupled Model Intercomparison Project phase 6 (CMIP6, (Eyring et al., 2016)) indicate an increase in the anomalously high-pressure systems over Europe (Bakke et al., 2023), this needs to be taken into account when analyzing the trends and variability of VPD for the next few decades. Our summer VPD reconstruction was also validated by examining the occurrence of the extreme years, the extremely hot and dry summers coincide with high VPD anomalies, while the extremely wet summers coincide with low VPD anomalies. Interestingly, the high VPD anomaly years show spatial variability with large positive anomalies in specific regions, while the wettest years display consistent negative VPD anomalies over the continent.

Furthermore, the interdisciplinary use of the data should be emphasized, as VPD is a crucial parameter for many climatological processes. As a logical next step, the regional and temporal boundaries of the reconstruction can be extended by using more and longer δ¹⁸O time series from tree ring cellulose. It is also possible to disentangle the influence of solar events on the VPD on local, regional, and continental scales.