Introduction

Recent examples of high-impact summertime weather and climate extremes such as extensive forest fires1,2, record-shattering heat waves3,4,5,6,7, and intense drought6,7,8, have become a hallmark of early twenty-first century climate change, causing US$ billions of insured losses every year9. A connecting element of these high-impact events are extremely large values of vapor pressure deficit (VPD), also referred to as air dryness or atmospheric water demand, which measures the difference between saturation (es) and actual water vapor pressure (ea) in air10, i.e., VPD = es − ea.

Extreme heat often relates to large VPD11, due to the exponential dependence of es on temperature12, and meteorological drought typically reduces the atmospheric moisture availability13,14,15. Thus VPD combines the meteorological variables temperature (T) and atmospheric moisture (q) in a way that is particularly relevant for ecosystems, e.g., forest carbon uptake and mortality6,16,17,18,19,20, wildfires21,22,23, and agricultural crop yields24,25,26. In the context of plant functioning and vegetation dynamics, longer-term extremes, e.g., of seasonal mean VPD in summer, are of central importance due to the sessile and partly long-lived nature of plant organisms27. Moreover, on the seasonal timescale, positive land-vegetation-atmosphere feedbacks can further amplify large VPD as a result of soil desiccation and reduced plant evapotranspiration28,29,30. A summer VPD extreme, i.e., an unusually large seasonal VPD anomaly in summer, can therefore arise from exceptionally large T and/or low q, which are linked to land-atmosphere coupling28,31 and circulation anomalies14.

Previous studies attributed tremendous socioeconomic consequences to recent summer VPD extremes (see affected regions in Fig. 1a): Two extensive VPD extremes in North America in 1988 and 2012 were associated with nation-wide agricultural disasters of $30–40 billion in economic losses32,33. The summer VPD extreme in 2011 was the most intense in the northern mid-latitudes (mean VPD anomaly of 0.78 kPa as per our definition, see “Methods”; Fig. 1b) and was directly linked to a record-extensive wildfire season in the southern U. S.21. Furthermore, cattle herds were diminished or forced to migrate north, where they were eventually hit by extreme VPD in the subsequent year of 2012, resulting in a total loss of 3.5 million head in Oklahoma and Texas33. Also vast regions of Europe were struck by summer VPD extremes in the last two decades (Fig. 1a), which impaired vegetation vitality (measured by vegetation indices) and reduced carbon uptake (e.g., in 200334, 201535, 201820,36, and 20226). Also, these extremes were associated with unseen wildfires (e.g., in 200737, 201838, and 20222), and substantial crop losses (e.g., in 201839). The so-called “angry summer” 2018/2019 in Australia marked the hottest on record40, during which the most intense summer VPD extreme of the southern mid-latitudes was observed (mean VPD anomaly of 1.16 kPa; Fig. 1b). Also several regions of the Global South were struck by summer VPD extremes, with particularly devastating consequences: In Argentina, the summer VPD extreme in 1988/1989 reduced the cultivated area by 35%41, and that in 2022/2023 diminished annual GDP by an estimated 3% due to lower crop yields, livestock, and forestry activity42. Moreover, summer VPD extremes have been suggested to be at the root of a water conflict between Turkey and Syria in 200043 and allegedly led to high mortality rates and migration of wildlife in several South African national parks in summer 1982/198344. In summary, the socio-ecological relevance of summer VPD extremes is obvious, however, it is currently not well understood to what extent and through which mechanisms these events are affected by climate warming. In this study we thus identify summer VPD extremes (i.e., the extremely large values) in reanalysis and climate model data as spatial objects (see “Methods”) and examine changes in their intensity as the climate warms.

Fig. 1: Selected summer vapor pressure deficit (VPD) extremes in 1979–2023 and the meteorological contributions to their intensity.
figure 1

a Selected summer VPD extreme objects, see “Methods” for details of the event identification. Each object is colored according to its year of occurrence and labeled accordingly. Some objects are additionally labeled with an identifier. Regions with gray hatching are not considered in this study. b Decomposition of the intensity (I) for the summer VPD extremes labeled A to I in (a). NH and SH refer the Northern and Southern Hemisphere. Red, dark blue, light blue, and gray bars denote contributions to I from positive temperature (T) anomalies (IT), T-related specific humidity (q) anomalies (IqT), dynamically induced q anomalies (Iqd), and a residual (res), respectively (see “Methods” for detailed description of the individual contributions).

Full size image

Climatological mean VPD has been increasing in the past four decades and is expected to rise further throughout the twenty-first century16. That is, global warming causes a nonlinear increase of es of 7% per Kelvin12, while the response of ea is in many regions less than 7%/K due to limited atmospheric moisture transport and land evapotranspiration13,14,45. The past increase in mean VPD and its impact on wildfire has recently been attributed directly to the emissions of major carbon producers46. Also for heat and drought47, both strongly coupled to large VPD, as well as their compound occurrence48 significantly increased risks of extreme events to occur have been investigated in detail and attributed to global warming. However, changes of summer VPD extremes, which add to already considerable changes in mean VPD, have not been systematically assessed yet. In light of the devastating impacts of summer VPD extremes this research gap should urgently be addressed.

Here we apply an identification algorithm for seasonal VPD extremes to both reanalysis and climate model data. This algorithm, which has previously been employed to identify seasonal temperature, wind and precipitation extremes49,50, consists of the following steps: First, an appropriate statistical model is fitted to the detrended seasonal mean VPD values at each grid point. Then, locally very rare (i.e., extreme) seasonal VPD values are identified as values exceeding the local 40-year return level. Finally, the algorithm produces spatial extreme event objects by connecting grid points that feature locally extreme VPD values in a particular season and year (see Fig. 1 for examples of spatial seasonal VPD extremes objects and “Methods” for more details on the identification algorithm). We first use data from ERA5 reanalyses51 as well as from the Community Earth System Model version 2 large ensemble (CESM2-LENS)52 for an “evaluation period” (1979–2023, the respective CESM2 data is hereafter referred to as CESM2eval) to show that the CESM2 model53 realistically reproduces the T and q contributions to the identified summer VPD extremes. We then address the effect of climate warming on VPD extremes by contrasting CESM2-LENS summer VPD extremes identified in the 1991–2020 and 2071–2100 periods (hereafter referred to as CESM2hist and CESM2eoc datasets, respectively). Overall, we tackle the following three research goals: (1) analyze spatial patterns of the T and q contributions to the intensity of summer VPD extremes (see also Fig. 1b), (2) quantify by how much the intensity of very rare summer VPD extremes changes under global warming on top of the changes in summer mean VPD, and (3) quantify the contributions of the dominant physical processes to these changes.

Results

Historical summer VPD extremes

We identify 139 and 29 summer VPD extremes (VPDS+) that affect at least 105 km2 land area with a center of mass in the northern (20–60°N; NH) and southern mid-latitudes (20–60°S; SH) in the ERA5 reanalysis, respectively. Figure 2a summarizes the average intensity (I; land average detrended VPD anomaly) over all these VPDS+, and the shares of I attributed to three potential sources of positive VPD anomalies (also see Eq. (13) in the “Methods”): warm temperature anomalies (\(T^{\prime} > 0\); contribution IT), T-related deficits in atmospheric moisture (\(q^{{\prime} }_{T} < 0\); contribution IqT), and negative anomalies in atmospheric moisture arising from atmospheric dynamics (\(q^{{\prime} }_{d} < 0\); contribution Iqd):

$$I={I}_{T}+{I}_{qT}+{I}_{qd}+res$$
(1)
Fig. 2: Meteorological decomposition of intensity (I) in the evaluation period (1979–2023).
figure 2

a The decomposition of I averaged over all summer vapor pressure deficit extremes (VPDS+), separately for the extratropical northern and southern hemisphere (NH and SH), for ERA5 (orange framed bars) and CESM2eval (black framed bars). CESM2eval values are computed as the mean over 1000 samples of equal size (number of VPDS+ in ERA5), with whiskers showing the resulting 95% confidence interval. Values above the bars indicate the respective share of I, with asterisks marking statistically significant differences between ERA5 and CESM2eval (5% significance level). The mean intensity (I; b), its contribution from temperature anomalies \(T^{\prime}\) (IT; c), from the climatological co-variability between temperature and specific humidity (IqT; d), and dynamically induced humidity anomalies \(q^{\prime}\) (Iqd; e) of all VPDS+ in CESM2eval. Thereby the mean value at every grid cell is calculated over all VPDS+ that comprise the respective grid cell.

Full size image

The \(q^{{\prime} }_{T}\) are estimated from an exponential regression of q with predictor T, hence quantifying the climatologically expected response of q to a given T (see “Methods”). In contrast, \(q^{{\prime} }_{d}\) are not related to T but are the result of anomalous horizontal/vertical moisture transport and other sources of anomalous moisture. Such sources include, for example, legacies of previous rain deficits/surpluses that are mediated via soil moisture anomalies. The variability in these quantities across VPDS+ is depicted in Supplementary Fig. 1a.

In both hemisphere’s mid-latitudes, the largest share of the positive VPD anomalies during VPDS+ is due to positive temperature anomalies. In the NH, IT contributes on average 75% to the intensity I = 0.30 kPa. Moreover, warm anomalies indirectly offset 9% of I via the IqT term. This implies that, on average in the NH, a moistening occurs during warmer than usual summers. Hence, dynamically induced moisture anomalies are another important source of 31% of I in ERA5 VPDS+. The VPDS+ in the SH are by 0.15 kPa more intense than those in the NH due to larger contributions IT and IqT (Fig. 2a). That is, warm anomalies directly account for IT = 65% of I and indirectly for another IqT = 14%, which results from the negative q-T-correlation in many SH regions in summer (see Supplementary Fig. 2a). Note that there is very large variability in the relative contributions of these terms to I across individual events, with IT ranging from 26% to 203% and Iqd ranging from  −14% to 89%.

The intensity of the 13,008 (NH) and 3625 (SH) VPDS+ identified in the 4500-year CESM2eval dataset, as well as its contributions, largely agree with the above results from ERA5 (Fig. 2a). Good agreement is also found for the spatial pattern and magnitude of the 40-year return level of summer VPD anomalies in both datasets (Supplementary Fig. 3). Using a bootstrapping test, we find that in the NH, the overestimation of I resulting from more positive IT and Iqd in CESM2eval is significant at the 5% level, but the overestimation is less than 17% of the respective value in ERA5. In the SH, summer VPD extremes feature slightly smaller I, smaller IT and IqT, and larger Iqd in CESM2eval compared to ERA5, whereby only the biases in IqT and Iqd are significant. These CESM2 biases are likely related to biases in the climatological q-T-correlations (Supplementary Fig. 2a, b; see “Methods”). Nevertheless, the level of agreement in the magnitude and composition of VPDS+ is remarkable, in particular given (1) that ERA5 does not have dynamic vegetation, while the Community Land Model version 5 used in CESM2 does, and (2) given the significant challenge of adequately simulating the climatological T-q-correlation, which is used here to differentiate IqT from Iqd (Supplementary Fig. 2). Overall, this analysis thus lends credibility to simulated changes in the magnitude of VPDS+ with global warming.

In the following, we assess the spatial variability of the simulated VPDS+ in CESMeval. The CESMeval offers a basis to systematically contrast characteristics of the by definition rare VPDS+ across space (Fig. 2b–d), which is not possible based on reanalysis data alone due to too few events per grid cell. Moreover, note that the spatial variability of VPDS+ in CESM2eval is almost identical to that in CESM2hist used in the next section, due to the large overlap of the considered periods (not shown). The most intense VPDS+ in CESM2eval with I > 0.7 kPa occur in the continental regions of North America and in the (semi-)arid regions of Australia and southwestern Asia (Fig. 2b). This is in accordance with the hemispherically most intense observed (i.e., according to reanalysis) VPDS+ that occurred in the U. S. in 2011 (I = 0.78 kPa) and in Australia in 2018/19 (I = 1.16 kPa; Fig. 1b). In those regions, all three meteorological anomalies contribute positively to I of VPDS+ in CESM2eval with the largest shares from IT (Fig. 2c–e), which was also the case for the highlighted VPDS+ in ERA5 (Fig. 1b). At higher latitudes, e.g., in parts of Canada, South America, and Eurasia, IqT of VPDS+ in CESM2eval is usually negative as T and q in summer correlate positively (Fig. 2d and Supplementary Fig. 2b). Overlapping with these regions, a rare group of VPDS+ (accounting for 1% of all simulated VPDS+) occurs in the Tibetan Plateau and along the lower-latitude east coasts of the continents during colder than usual summers. They hence exhibit positive IqT, which together with dynamical dry anomalies (positive Iqd) over-compensate the negative IT (not shown). Such a combination of meteorological contributions has not been observed in ERA5, but the typically negative IqT and the pronounced contributions from Iqd were also apparent in the case studies B-E and H (Fig. 1b). Generally speaking, in CESM2eval, IT tend to be smaller over eastern portions of the continents compared to the central and western portions. Accordingly, VPDS+ in these eastern continental regions are comparably less intense (I ≈ 0.2 kPa) than their counterparts in the central and western continental regions. All in all, VPDS+ in many regions acquire most of their intensity through warm anomalies but moisture-related contributions are often non-negligible.

Intensification of VPDS+ with global warming

Next we use the 100 member CESM2-LENS data to assess intensity changes of VPDS+ under global warming, by comparing their I and its composition between historical (hist) and end of the century (eoc) 30-year periods (1991–2020 and 2071–2100). The CESM2-LENS simulations are forced with historical forcing up to 2014 and SSP3-7.0 forcing thereafter52. Note that VPDS+ are defined by rarity in their respective local climate and period. That is, a comparable number of events is identified in both the hist and eoc period. On average, VPDS+ intensify by +33% and +28% in the NH and SH, respectively, and exhibit, on average, an intensity of I = 0.49 kPa across both mid-latitudes in the end-of-century period (Fig. 3a, the variability across events is depicted in Supplementary Fig. 4a). The intensification is slightly stronger in the NH showing an intensification in excess of 40% in central to northern Europe and parts of North America (ΔI ≈ +0.2 kPa; Fig. 3a, b). Recall that this intensification of extremes occurs on top of the changes in the climatological mean (VPDc), which on average equals 0.54 kPa and 0.41 kPa in the NH and SH, respectively (shown in Supplementary Fig. 5b). The absolute 40-year return level, i.e., 40-year anomaly plus the respective VPDc, increases by around 22% more strongly than VPDc alone (averaged over all regions considered here). Hence, summer VPDc will increase considerably, but seasonal VPD variability will too, with drastic consequences for VPDS+. The intensification is dominated by changing temperature effects, i.e., increasing IT in both the NH and SH (Fig. 3a, c), although in some regions it is slightly offset or further exacerbated by changes in the IqT or Iqd terms. Next, we thus investigate the physical causes of changes in I (hereafter ΔI) in greater detail, with a focus on changes in IT (ΔIT).

Fig. 3: Changes in intensity (I) between CESM2hist and CESM2eoc.
figure 3

a The decomposition of I averaged over all summer vapor pressure deficit extremes (VPDS+) in the northern and southern hemisphere (NH and SH) for CESM2hist (black framed bars) and CESM2eoc (purple framed bars). Values above the bars indicate the respective share of I. The changes in mean I (b), in the contribution related to temperature anomalies (IT) (c), the contribution related to the climatological co-variability between temperature and specific humidity (IqT) (d), and the contribution from dynamically induced humidity anomalies (Iqd) (e) of VPDS+ in CESM2eoc minus in CESM2hist, whereby the mean value at every grid cell is calculated over all VPDS+ that comprise the respective grid cell.

Full size image

The ΔIT (i.e., the temperature related changes in the intensity of VPDS+) can be further decomposed into contributions from two processes (see “Methods”): firstly, a process \({{\mathcal{T}}}_{{\rm{clim}}}\) that quantifies the VPD anomaly change related to the non-linear dependence of VPD on T, i.e., changes in the efficiency with which a given \(T^{\prime}\) generates VPD anomalies. Secondly, changes in \(T^{\prime}\) during VPDS+ themselves, termed \({{\mathcal{T}}}_{{\rm{var}}}\). These two processes can be quantified by a linearization of the respective terms over their values in the hist period (see “Methods”, Table 1 and Fig. 4). To first order, \({{\mathcal{T}}}_{{\rm{clim}}}\) and \({{\mathcal{T}}}_{{\rm{var}}}\) together equal ΔIT except for minor terms that are, for ease of interpretation, attributed to a residual \({\mathcal{E}}\) (see “Methods”). Finally, we introduce the symbol \({\mathcal{Q}}\), which denotes the sum of the changes in the two moisture related quantities ΔIqT and ΔIqd, to arrive at a more process based decomposition of ΔI (Table 1, Eq. (15)). Hereby, the residual \({\mathcal{E}}\) comprises comparably small terms that are less straightforward to interpret (see Supplementary Method 1).

Table 1 The processes contributing to the intensification of summer vapor pressure deficit extremes (VPDS+)
Full size table
Fig. 4: Processes contributing to the intensification of summer vapor pressure deficit extremes (VPDS+) between the hist and eoc period.
figure 4

The mean relative shares of intensity change (ΔI) contributions from the four processes in Table 1 (amber for \({{\mathcal{T}}}_{{\rm{clim}}}\), red for \({{\mathcal{T}}}_{{\rm{var}}}\), blue for \({\mathcal{Q}}\), and gray for \({\mathcal{E}}\)) in the northern (a) and southern (b) hemisphere (NH and SH). c–h Values of relevant VPDS+-averaged quantities again averaged over all VPDS+ in CESM2hist (light colors) and CESM2eoc (darker colors) in the NH and SH, respectively. Panels depict (c, f) the VPD anomaly induced by \(T^{\prime}=1\) K, (d, g) \(T^{\prime}\) of the respective VPDS+, and (e, h) \(q^{\prime}\) of the respective VPDS+. i The dominant process contributing to ΔI (i.e., largest contributor), whereby the attribution at every grid cell is calculated over all VPDS+ that comprise the respective grid cell. White stippling shows grid cells with ΔI < 0.05 kPa. Inset panels show the four processes as in (a, b) but in absolute terms at seven locations (see also Supplementary Fig. 6), which are denoted by black crosses (y-axis extends from 0 kPa to 0.15 kPa). See text and “Methods” for details.

Full size image

Figure 4 shows the process attribution of ΔI (Fig. 4a, b; see also Supplementary Fig. 6) alongside changes in the underlying physical quantities of VPDS+ (Fig. 4c–h) and the dominant process at each grid cell (Fig. 4i). About half of ΔI in both hemispheres is due to \({{\mathcal{T}}}_{{\rm{clim}}}\), i.e., the increasing efficiency of \(T^{\prime}\) to induce large VPD (Fig. 4a, b). That is, on average a \(T^{\prime}=1\) K during VPDS+ in CESM2eoc increases VPD by 0.23 kPa in the NH and SH, while it increases VPD only by 0.18 kPa during VPDS+ in CESM2hist (Fig. 4c, f). This is due to the non-linear dependence of VPD on T, which strongly increases the efficiency of \(T^{\prime}\) to induce large VPD’ as the global mean temperature increases by  +3.1 K in JJA and  +3.2 K in DJF. In many mid-latitude regions \({{\mathcal{T}}}_{{\rm{clim}}}\) is the dominant process, and in arid and continental regions it is almost exclusively responsible for ΔI (Fig. 4i and Supplementary Fig. 6).

The second T-related process, \({{\mathcal{T}}}_{{\rm{var}}}\), which captures changes in the T variability during VPDS+, does on average not contribute as much as \({{\mathcal{T}}}_{{\rm{clim}}}\), but at regional scales can contribute the major share of ΔI (Fig. 4i). Averaged over all VPDS+, \(T^{\prime}\) increases from  +1.5 K in the NH and from  +1.4 K in the SH to  +1.6 K, which translates to 15% and 24% of the intensification of I in the NH and SH, respectively (Fig. 4a, b). Increases in \(T^{\prime}\) during VPDS+ are related to an increase in the interannual T variability. Interestingly, in central and northern Europe as well as in parts of southeastern North America and of southern Africa, most of ΔI is attributed to \({{\mathcal{T}}}_{{\rm{var}}}\) (Fig. 4i and Supplementary Fig. 6).

Lastly, changes in the moisture contributions from CESM2hist to CESM2eoc (summarized as \({\mathcal{Q}}\)) contribute 22% and 15% of ΔI in the NH and SH, respectively. In many regions, the q-T-correlation becomes more negative and thus \(q^{\prime}\) is lower during warm summers in the end-of-century climate compared to the historical climate (see Fig. 4e, h and Supplementary Fig. 2d). This is evident in strong relative increases in IqT (i.e., the VPD’ resulting from T-related \(q^{\prime}\)) adjacent to (semi-)arid regions (Supplementary Fig. 4d). Regions that maintain a high q-T-correlation in the warmer end-of-century climate, e.g., the northeastern continental or high-altitude regions, exert an increase in climatological mean q, which also increases the q variability (see “Discussion”). Also there, \({\mathcal{Q}}\) can even be locally the most important contribution to ΔI (Fig. 4i). Note that some regions in southeastern Asia, South America, and New Zealand show a relatively weak or no intensification (ΔI < 0.05 kPa; Fig. 4i). In summary, at a global scale the intensification of VPDS+ with global warming arises mainly from \({{\mathcal{T}}}_{{\rm{clim}}}\), i.e., the same \(T^{\prime}\) in a warmer climate produce more intense VPD anomalies, while a more positive \(T^{\prime}\) and more negative \(q^{\prime}\) during VPDS+ are both adding to this intensification in a regionally differing manner.

Discussion

This study highlights the socioeconomically highly relevant fact that very rare summer VPD extremes defined with a 40-year return period threshold (VPDS+) intensify drastically in a warming climate. Importantly, about half of that intensification would occur already if the variability in their meteorological contributors \(q^{\prime}\) and \(T^{\prime}\) were to remain constant. This is due to exponential dependence of VPD on T, which implies that the same \(T^{\prime}\) in a warmer climate will cause larger VPD anomalies (process \({{\mathcal{T}}}_{{\rm{clim}}}\)). Larger VPD variability in a warmer climate has been documented previously in projections of VPD in the U. S.22,54 and analyses of shorter-term VPD extremes55. Here, we substantially extend these studies by highlighting not only the role of mean warming separately from that of \(T^{\prime}\), but also by quantifying its contribution to ΔI of VPDS+ with 54% in the NH and 49% in the SH.

Of further importance for the ΔI in both mid-latitudes are increasing \(T^{\prime}\) during VPDS+ (process \({{\mathcal{T}}}_{{\rm{var}}}\)) and changes in the moisture contributions to I (process \({\mathcal{Q}}\)). The particularly strong intensification of VPDS+ in northwestern Europe relates to the increasing interannual variability in summer T (fostering \({{\mathcal{T}}}_{{\rm{var}}}\)) that has previously been related to land-atmosphere feedbacks during hot-dry summers56,57. Process \({\mathcal{Q}}\) includes mainly decreases in \(q^{\prime}\) during VPDS+. On the one hand, more negative \(q^{\prime}\) arise from less positive q-T-correlations in a warming climate, which have been associated with the expansion of soil-moisture limited regions (in the transitional climate regime)29, where positive \(T^{\prime}\) tend to occur with limited evapotranspiration and hence negative \(q^{\prime}\)30,31. Such a shift from a wet to transitional climate regime is most widespread in central Eurasia and central North America29,31, where \({\mathcal{Q}}\) is strongest (Supplementary Fig. 6). On the other hand, for a bounded variable like q its variability increases as a result of the increase in the climatological mean q, as has been observed for precipitation45,58. In summary, both changes in interannual summer T variability as well as regime-shifts in the land-atmosphere coupling also amplify the magnitude of VPDS+.

We have shown that, at large scales, CESM253 is able to realistically reproduce physical characteristics of the high-impact extreme VPD summers. Even though an accurate representation of characteristics of present-day VPDS+ does not imply that simulated changes in VPDS+ must necessarily be accurate, this result nevertheless lends credibility to these simulated changes. The ability of GCMs to inform about future high-impact VPDS+ has hitherto not been explored thoroughly. Evidently, the representation of VPDS+ in GCM simulations should be evaluated for various models, and our analyses should be repeated with other state-of-the-art GCMs and emission scenarios to corroborate our findings. However, note that even though our results were derived with only one GCM and emission scenario, we expect that the key finding of this study (VPDS+ are intensifying drastically even relative to increasing mean VPD) is robust across models. Around 50% of this intensification is attributable to the exponential dependence of VPD on T, which implies that in a warmer climate, the same \(T^{\prime}\) yields larger VPD (process \({{\mathcal{T}}}_{{\rm{clim}}}\))—this part of the intensification is essentially a consequence of basic physics and therefore most likely independent from the particular model that is used (albeit its magnitude might vary with the climate sensitivity of a particular model). Furthermore, at a regional scale (in particular in northwestern Europe), the pronounced increase in inter-annual summer T variability is robust across different global and regional climate models56,57. Hence, we expect qualitatively similar \({{\mathcal{T}}}_{{\rm{var}}}\) across GCMs, while the additional effects such as changes in summer q variability or regime-shifts in the land-atmosphere coupling might be less robust across models59,60,61.

Our results suggest that VPDS+ with equal rarity will intensify drastically as the climate warms, which has severe implications for ecosystems. Put in other words, summer VPD anomalies of a given magnitude will become far more frequent in the future, and thus ecosystems not only need to adapt to larger mean values of summer VPD16 but they also need to become resilient to larger positive anomalies of summer VPD. Given the drastic impacts of seasonal VPD extremes already in the current climate, we argue that this result should be considered by practitioners, e.g., when managing food producing systems, (inter-)national food security or forest ecosystems, or when designing agricultural reinsurance. Moreover, regarding certain impacts such as the Black Summer Wildfires in Australia in 2019/20201, a single VPDS+ might have a delayed impact or one that arises only after sequential VPDS+1,40. Therefore, future work should also examine consecutive or longer-than-seasonal VPD extremes so that practitioners can adequately incorporate the effects of VPD extremes in their preparedness strategies.

Methods

Datasets

We use hourly ERA5 reanalysis data on a horizontal 0.5° grid from December 1978 to August 202351. Seasonal mean values are computed first and then bi-linearly interpolated to the CESM2 grid. Moreover, in either of the considered datasets December to February (DJF) seasons assigned to any year contain December data from the previous year. Specific humidity at 2 m, denoted as q, is calculated from 2-m dew point temperature and surface pressure (p)62. The 100 member CESM2-LENS52 uses historical forcing in 1979–2014 and follows the SSP3-7.0 protocol from CMIP663 from 2015 onwards. It features a horizontal resolution of 1.25° longitude × ~0. 9° latitude. Individual model components as well as the initialization of the 100 members is thoroughly described in Rodgers et al.52.

To minimize effects of conditional bias49,64 in our statistical modeling of ERA5 and CESM2 VPD data (see section “Identification of VPDS+”), we group CESM2 data into subsets, each with the same number of years as the ERA5 data has, and then process these subsets of CESM2 data individually. That is, for identifying VPDS+ in CESM2eval, which consists of 100 times 45 years of data (1979–2023) each of the CESM2-LENS members is processed individually, but exactly as ERA5 data. Note that for this 45-year evaluation period, the CESM2 simulations combine historical forcing until 2014 and SSP3-7.0 forcing thereafter. For CESM2hist and CESM2eoc we combine one member (i.e., 30 years) with half of another member (i.e., 15 years) so that the resulting subset again contain 45 years of data. This procedure yields 66 subsets of 45 years (i.e., 2970 years in total) for both CESM2hist (in 1991–2020) and CESM2eoc (in 2071–2100). These subsets too are then processed exactly analogous to ERA5 data in order to identify the respective VPDS+. Note that keeping the number of years constant in each subset ensures comparability between the VPDS+ in CESM2 data and ERA5, as the sample size affects the skill with which extreme events with large return periods can be identified as such through our statistical modeling49. We use JJA and DJF (December of the previous year) mean values in the NH and the SH, respectively, and refer to them as summer mean values. Moreover, we provide a short model evaluation based on our results and the consideration of relevant literature in the last subsection below.

Computation and detrending of VPD

We derive a formulation of VPD [Pa] in terms of T [K], q [kg kg−1], and p [Pa] starting with the definition of VPD according to Grossiord et al.10. For es(T) we use Murray’s formulation65, which is valid over a wide range of T, and for ea(q, p) we follow the derivation by Allen et al.66:

$${\rm{VPD}}={e}_{s}-{e}_{a}=610.78\cdot \exp \left(\frac{17.2694\cdot (T-273.16)}{T-35.86}\right)-\frac{q\cdot p}{0.622+0.378\cdot q}$$
(2)

Note that VPD is computed from seasonal mean values. Next, we compute detrended VPD anomalies that are physically consistent with detrended T and q (Fig. 5). Following Sutton et al.67, we first compute a least square regression for local T with global mean summer T (GMST, computed separately in each CESM2 subset as in ERA5) as predictor (Fig. 5a):

$$T={\alpha }_{1}+{\alpha }_{2}\cdot {\rm{GMST}}+{\epsilon }_{T}={T}_{c}+{\epsilon }_{T}$$
(3)

resulting in the two regression parameters (α1, α2) and Tc, which is a forced trend estimate for T, i.e., a temporally evolving climatological reference. As physics-based scaling arguments suggest an exponential relationship between q and T 68, we estimate the forced trend in q based on an exponential fit of q on Tc (Fig. 5d). To do so we fit a linear least square regression to \(\ln (q)\) with Tc as predictor. This results in regression parameters β1 and β2, which are used to estimate the trend in q, i.e., qc:

$$\ln (q)={\beta }_{1}+{\beta }_{2}\cdot {T}_{c}+{\epsilon }_{q}$$
(4)
$${q}_{c}=\exp ({\beta }_{1}+{\beta }_{2}\cdot {T}_{c})$$
(5)
Fig. 5: Illustration of the forced trend estimation for ERA5 summer mean values of temperature, specific humidity, and vapor pressure deficit at the example grid cell at 54 °N, 67 °E.
figure 5

Scatter plots show temperature (T) vs. Global mean surface T (GMST) (a), specific humidity (q) vs. the forced trend estimate for T (Tc) (d), and vapor pressure deficit (VPD) vs. GMST (g), including the respective trend curves in black. Inset numbers in (a, d) mark the slope of the trend line, which is expressed in relative terms for q(Tc). b, e, h Evolution of the respective values (colored lines) including their trend estimate (black line) in 1979–2023. c, f, i Detrended anomalies, i.e., deviations from the trend line.

Full size image

Finally, we assume no trend in p with respect to GMST due to the weak dependence of VPD on p. Based on the time-evolving climatological values of q and T we compute the climatological VPD by using Eq. (2) (Fig. 5g):

$${{\rm{VPD}}}_{{\rm{c}}}={\rm{VPD}}({T}_{c},{q}_{c},p)$$
(6)

Deviations from the temporally evolving climatological reference values VPDc, Tc, and qc are assumed to be entirely due to internal variability. Hence, for X ∈ {VPD, T, q} the (detrended) seasonal mean anomaly (\(X^{\prime}\)) equals (Fig. 5c, f, i):

$$X^{\prime}=X-{X}_{c}$$
(7)

Identification of VPDS+

The VPDS+ in ERA5 and all CESM2 subsets are identified using a scheme targeted to seasonal extremes, which has previously been applied to identify seasonal extremes in T, surface winds, and precipitation49,50, and consists of three steps: (1) Estimate at each grid cell the distribution of VPD’ by statistical modeling and, based on that, calculate the local return period (LRP) of each VPD’ value, (2) identify locally extreme VPD’ based on exceedances of an LRP threshold τ, and (3) form spatially coherent objects within which LRP > τ and quantify their characteristics. As VPD has a lower bound of zero, the distributions of VPD’ are typically positively skewed. To account for this skew in step (1) we fit a Yeo-Johnson transformed normal distribution (NYJ) with the three parameters μ (mean), σ2 (variance), and λ (transform parameter) to VPD’69, which has previously been used to model seasonal mean T and near-surface wind speed49,50. The three parameters are allowed to vary between grid cells and datasets, and are estimated with a maximum likelihood estimation. As for wind and T, also for seasonal VPD a goodness-of-fit test reveals no widespread departures of the data from the Yeo-Johnson transformed normal distribution in ERA5 and the CESM subsets (Supplementary Fig. 7). Furthermore, in steps (2) and (3), we use τ = 40 years, following Boettcher et al.50. Hence every grid cell experiences on average about one VPDS+ in ERA5 and in each CESM2 subset (see Supplementary Fig. 8). Further details regarding the identification of seasonal extremes are provided by the authors of the framework49. For our analysis we only consider VPDS+ with a center of mass in 20–60°S (SH) and 20–60°N (NH) and a land area larger than 105 km2. This yields 168 VPDS+ in ERA5, 16,633 in CESM2eval, 11,063 in CESM2hist, and 10,948 in CESM2eoc.

Intensity of VPDS+ and its meteorological decomposition

The central aspect of every VPDS+ studied here is its intensity (I), which we define as the average VPD’ across all land grid cells contained in the respective VPDS+. In the following, we decompose I into contributions from \(T^{\prime}\), \(q^{{\prime} }_{T}\), and \(q^{{\prime} }_{d}\). To do so, we first use a multivariable linearization of VPD from the climatological reference, i.e., we approximate any grid cell value VPD(T, q, p) tangentially from VPDc = VPD(Tc, qc, p)70:

$${\rm{VPD}}(T,q,p) \, \approx \, {\rm{VPD}}({T}_{c},{q}_{c},p)+(T-{T}_{c}) \cdot \left.\frac{\partial {\rm{VPD}}}{\partial T} \right|_{{{T}_{c}},{q}_{c},p} \\ +(q-{q}_{c})\cdot \left.\frac{\partial {\rm{VPD}}}{\partial q}\right|_{{{T}_{c}},{q}_{c},p}$$
(8)

which can be formulated as:

$${\rm{VPD}}^{\prime}=T^{\prime} \cdot \left.\frac{\partial {\rm{VPD}}}{\partial T} \right|_{c}+q^{\prime} \cdot \left.\frac{\partial {\rm{VPD}}}{\partial q} \right|_{c}+{\epsilon }_{VPD}$$
(9)

where ϵVPD marks a residual term containing higher-order terms, and the derivatives of VPD(T, q, p) wrt. T and q are evaluated at (T, q, p) = (Tc, qc, p) according to Eq. (2), for simplicity referred to with “c”. In a second step, we account for the climatological co-variability of T and q. That is, in some regions warmer than usual summers relate to a limited moisture supply (\(q^{\prime} < 0\) when \(T^{\prime} > 0\)), and in others moistening results from the increased water holding capacity (\(q^{\prime} > 0\) when \(T^{\prime} > 0\)). To estimate the climatological co-variability of T and q, we compute an exponential regression for q with predictor T:

$$\ln (q)={\gamma }_{1}+{\gamma }_{2}\cdot T+{\epsilon }_{qT}$$
(10)
$${q}_{T}=\exp ({\gamma }_{1}+{\gamma }_{2}\cdot T)$$
(11)

So the regression parameters γ1 and γ2 are used to compute the climatologically expected q given T (qT). Again the exponentiated regression structure is chosen due to the exponential dependence of q on T derived from physical scaling arguments (i.e., Clausius-Clapeyron relation). Next, we rewrite the expression for the q-anomaly (\(q^{\prime}\)) as

$$q^{\prime}=q-{q}_{c}=(q-{q}_{T})+({q}_{T}-{q}_{c})=q^{{\prime} }_{d}+q^{{\prime} }_{T}$$
(12)

that is, \(q^{\prime}\) is the sum of the deviation of the actual q from the expected q for a given temperature (qT) plus the difference of qT minus the climatological q. The latter difference is what we refer to as \(q^{{\prime} }_{T}\), i.e., the \(q^{\prime}\) resulting from the climatological co-variability with T. We assume the remaining part of \(q^{\prime}\), i.e., \(q^{{\prime} }_{d}=q-{q}_{T}\), to be the result of atmospheric dynamics and land-atmosphere interactions.

To obtain the decomposition of I for a single VPDS+, we combine Eq. (9) with Eq. (12) and average over all land grid cells contained in that VPDS+ (denoted by square brackets):

$$I=[{\rm{VPD}}^{\prime} ]=\left[T^{\prime} \cdot \left.\frac{\partial {\rm{VPD}}}{\partial T}\right| _{c}\right]+\left[q^{{\prime} }_{T}\cdot \left.\frac{\partial {\rm{VPD}}}{\partial q}\right|_{c}\right]+\left[q^{{\prime} }_{d}\cdot \left.\frac{\partial {\rm{VPD}}}{\partial q}\right|_{c}\right] \\+[{\epsilon }_{VPD}]={I}_{T}+{I}_{qT}+{I}_{qd}+res$$
(13)

where the four terms on the r.h.s. quantify the impact of each meteorological anomaly on I in [Pa]. Note that in the main text, we use the term I and its contributions when referring to the intensity of a single VPDS+ and to the average over all VPDS+ in one hemisphere. The context of the notation thus specifies whether I and its contributions refer to a single (e.g., in Fig. 1b) or multiple VPDS+ (e.g., in Fig. 2).

Process attribution of changes in I

A common approach to attribute temporal changes in a given quantity to its components is based on a linearization of the respective relationship71. We apply this framework to the change in IT averaged over all VPDS+ in one hemisphere between the historical and end-of-century period. First, we contrast the terms in Eq. (13) between the two periods:

$$\Delta I={\langle I\rangle }^{eoc}-{\langle I\rangle }^{hist}=\Delta {I}_{T}+\Delta {I}_{qT}+\Delta {I}_{qd}+\Delta res$$
(14)

where angle brackets with superscript eoc and hist denote averages over all considered VPDS+ in the given period. The Δ-terms refer to the mean changes of a given contribution over all considered VPDS+ between the two periods. Next, we re-write all components of IT evaluated in the end-of-century period (Eq. (13)) as the sum of the respective values in the historical period plus the respective deltas, e.g., \({\langle [T^{\prime} ]\rangle }^{eoc}={\langle [T^{\prime} ]\rangle }^{hist}+\Delta T^{\prime}\). Hence each Δ-term quantifies the mean change of a quantity over all land grid cells of the considered VPDS+ between the two periods. These steps are detailed in Supplementary Method 1 and yield a refined form of Eq. (14):

$${{\Delta}} I=\underbrace{{{\Delta}} T^{\prime} {\cdot} {\left\langle \left[ \left.\frac{\partial {\rm{VPD}}}{\partial T}\right|_{c} \right] \right\rangle}^{hist}}_{{{\cal{T}}}_{\rm{var}}}+\underbrace{{{\Delta}} \left.\frac{\partial {\rm{VPD}}}{\partial T}\right|_{c} {\cdot} \langle [ T^{\prime} ] \rangle^{hist}}_{{{\cal{T}}}_{\rm{clim}}}+\underbrace{{{\Delta}} I_{qT}+{{\Delta}} I_{qd}}_{{{\cal{Q}}}}+{{\cal{E}}}$$
(15)

which attributes the total ΔI to changes from the so-called processes \({{\mathcal{T}}}_{{\rm{clim}}}\), \({{\mathcal{T}}}_{{\rm{var}}}\), \({\mathcal{Q}}\), and \({\mathcal{E}}\). Note that \({{\mathcal{T}}}_{{\rm{clim}}}+{{\mathcal{T}}}_{{\rm{var}}}\approx \Delta {I}_{T}\) with differences only arising from covariance and higher-order terms, which are contained in the residual \({\mathcal{E}}\). All such details regarding the computation and definition of the four processes are available in Supplementary Method 1. In the main part of the study the processes are introduced more conceptually in Table 1 in accordance with Eq. (15).

Significance assessment

We conduct a bootstrapping test to identify statistically significant differences in I and its contributions between VPDS+ in ERA5 and in CESM2eval. The procedure draws 1000 random samples of n VPDS+ from CESM2eval, whereby n is the number of VPDS+ under consideration in ERA5 (n = 139 in the NH, n = 29 in the SH). We test the null hypothesis that the VPDS+ in ERA5 and CESM2eval stem from the same distribution of VPDS+. The null distribution is constructed from the mean values over all random samples. The null hypothesis is rejected if the ERA5 mean value lies outside the 2.5th to 97.5th percentile range of the null distribution (significance level α = 5%). In that case, a statistically significant difference is identified for the respective characteristic of the VPDS+ in ERA5 and CESM2eval.

Model evaluation

The key finding of this study, namely the fact that very rare summer VPD extremes intensify drastically in a warming climate is likely robust across models, as it is mostly due to basic physics, i.e., the non-linear dependence of VPD on T. Nevertheless, changes in T variability and land-atmosphere coupling also affect the magnitude of VPDS+ intensity changes, and may be model-dependent. A more thorough model evaluation regarding VPD is thus warranted. Although the hemispherically aggregated I in CESM2eval is largely comparable to ERA5, and observed VPDS+ agree with the spatial patterns identified in CESM2eval, CESM2 reveals biases regarding some aspects of VPD. In the NH, \(q^{{\prime} }_{T}\) and \(q^{{\prime} }_{d}\) are influenced by regional biases in the q-T-correlation (Supplementary Fig. 2a, b) which partly compensate, hence resulting only in small discrepancies in the I decomposition between ERA5 and CESM2 (Fig. 2a). Biases in the land-atmosphere coupling have been identified in several older, CMIP5 generation models, and are inherently coupled to T biases72. In the SH, overestimated Iqd likely relate to the lower IqT and IT, since the q-T-correlation is too strong, i.e., warm summers show a weaker T-related drying (\(q^{{\prime} }_{T}\)) and hence attribute more of \(q^{\prime}\) to \(q^{{\prime} }_{d}\) instead. Lastly, model mean state biases add to the identified discrepancies. Nevertheless, the future drying trend in many mid-latitude regions (Supplementary Fig. 2d) is consistent across CMIP561 and CMIP6 models59,73, and their spatial patterns are in agreement with the identified VPDS+ changes, especially regarding the processes resulting from stronger land-atmosphere feedback and decreasing soil moisture in summer. Finally, it should be noted that, no matter which model one uses, the exact meaning of simulated two-meter T, q and VPD is not obvious in regions with tall vegetation. Here, however, we compare simulated T, q and VPD with reanalysis data (rather than direct observations), which should be more comparable to one another than climate model data and direct observations.