Significance Soil drought and atmospheric aridity can be disastrous for ecosystems and society. This study demonstrates the critical role of land–atmosphere feedbacks in driving cooccurring soil drought and atmospheric aridity. The frequency and intensity of atmospheric aridity are greatly reduced without the feedback of soil moisture to atmospheric temperature and humidity. Soil moisture can also impact precipitation to amplify soil moisture deficits under dry conditions. These land–atmosphere processes lead to high probability of concurrent soil drought and atmospheric aridity. Compared to the historical period, models project future frequency and intensity of concurrent soil drought and atmospheric aridity to be further enhanced by land–atmosphere feedbacks, which may pose large risks to ecosystem services and human well-being in the future.

Abstract Compound extremes such as cooccurring soil drought (low soil moisture) and atmospheric aridity (high vapor pressure deficit) can be disastrous for natural and societal systems. Soil drought and atmospheric aridity are 2 main physiological stressors driving widespread vegetation mortality and reduced terrestrial carbon uptake. Here, we empirically demonstrate that strong negative coupling between soil moisture and vapor pressure deficit occurs globally, indicating high probability of cooccurring soil drought and atmospheric aridity. Using the Global Land Atmosphere Coupling Experiment (GLACE)-CMIP5 experiment, we further show that concurrent soil drought and atmospheric aridity are greatly exacerbated by land–atmosphere feedbacks. The feedback of soil drought on the atmosphere is largely responsible for enabling atmospheric aridity extremes. In addition, the soil moisture–precipitation feedback acts to amplify precipitation and soil moisture deficits in most regions. CMIP5 models further show that the frequency of concurrent soil drought and atmospheric aridity enhanced by land–atmosphere feedbacks is projected to increase in the 21st century. Importantly, land–atmosphere feedbacks will greatly increase the intensity of both soil drought and atmospheric aridity beyond that expected from changes in mean climate alone.

Compound drought, heat and aridity events have recently received more attention because of their devastating impacts on the environment, economy and society (1, 2). Low soil moisture (SM) and high atmospheric vapor pressure deficit (VPD), an indicator of atmospheric aridity, have been recognized as 2 main stresses on ecosystem productivity during droughts (3). These stresses can substantially reduce terrestrial carbon uptake and food production (4⇓–6) and can drive widespread tree mortality (7, 8). Recently, it has been shown that low SM and high VPD typically cooccur, and the compound low-SM and high-VPD events are projected to be more frequent and more extreme in the future, which would strongly limit the capacity of continents to act as a carbon sink (9). Understanding the mechanisms underpinning compound drought and aridity events, specifically compound SM and VPD extremes, is of great importance to manage and minimize the risks associated with climate change (10).

Climate extremes are often interpreted as the result of large-scale atmospheric circulation and sea surface temperature anomalies (11, 12). For example, tropospheric warming induced by the El Niño Southern Oscillation and atmospheric blocking are key processes for drought and heatwave initiations (13, 14). In addition to large-scale atmospheric circulation anomalies forced by the ocean, land–atmosphere (LA) feedbacks initiated by SM anomalies can strongly modulate near-surface heat and aridity (15, 16), and also promote large-scale atmospheric circulation anomalies remotely (17, 18). Sustained SM deficit induced by negative precipitation anomaly and the related strong sensible heat flux indeed can cause extremely high temperature (16, 19). Lower evapotranspiration due to soil drying may also reduce atmospheric moisture, further increasing VPD. In turn, high VPD enhances the atmosphere’s evaporative demand, which can exacerbate SM depletion (20). Additionally, SM anomalies can induce precipitation anomalies, which may act to amplify or alleviate SM droughts (21). Given the potential importance of LA feedbacks in the occurrence of soil droughts and atmospheric aridity, it is crucial to systematically investigate the role of LA feedbacks in the occurrence of compound drought and aridity events.

LA feedbacks are expected to be enhanced by climate warming, as surface fluxes will be more sensitive to SM variability (22). The enhanced LA feedbacks may lead to stronger correlation between SM and VPD and hence more frequent compound drought and aridity events. As a result of anthropogenic warming, the intensities of soil droughts and atmospheric aridity are projected to increase (23, 24). The magnitude of intensity changes may also depend on LA feedbacks, which could impact the intensities of soil droughts and atmospheric aridity simultaneously, in addition to the changes driven by the long-term climate warming. Therefore, LA feedbacks might be critical for the changes in the frequency and intensity of compound drought and aridity events under future climate change.

In this study, we aim to assess the contributions of LA feedbacks to the compound drought and aridity events. As soil droughts can last for months to years, while atmospheric aridity typically occurs for weeks to months, we assess the compound drought and aridity events at the monthly scale. We first show that strong SM–VPD coupling occurs globally, using monthly VPD from the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) (25) and observationally constrained SM from the Global Land Evaporation Amsterdam Model (GLEAM) dataset (26), with important implications for the occurrence of compound drought and aridity events. We then use the Global Land-Atmosphere Climate Experiment-Coupled Model Intercomparison Project phase 5 (GLACE-CMIP5) experiment (27) to systematically assess the impacts of LA processes on concurrent soil drought and atmospheric aridity. This experiment includes pairs of transient simulations with and without interactive SM, while the boundary conditions are the same for each model (Materials and Methods). It thus offers a unique opportunity to isolate the role of SM–atmosphere interactions from the driving role of ocean–atmosphere variability. We further investigate the contributions of LA feedbacks to the frequency and intensity of compound drought and aridity events between historical and future periods in 26 CMIP5 models. Our findings demonstrate the fundamental role of LA feedbacks in the occurrence and exacerbation of compound drought and aridity events, both in the present and future climate.

SM–VPD Coupling and Cooccurring Low SM and High VPD in Observational Datasets We use the Pearson’s correlation coefficient between monthly SM and VPD [r(SM,VPD)] to measure the SM–VPD correlation in the warm season (the 3 mo with highest mean temperature; SI Appendix, Fig. S1). We find negative r(SM,VPD) globally, for both root-zone SM and surface SM, indicating that lower SM tends to be accompanied by higher VPD in the warm season (Fig. 1A and SI Appendix, Fig. S2A). In addition, surface SM is more strongly coupled with VPD than root-zone SM in most regions. The importance of the strong SM–VPD coupling for cooccurring SM and VPD extremes is supported by the bimodal distribution of SM and VPD toward both extremes: low SM–high VPD and high SM–low VPD (Fig. 1B and SI Appendix, Fig. S2B). The high probability in the top-left bin in Fig. 1B and SI Appendix, Fig. S2B indicates that the frequency of compound low SM and high VPD is much higher than that expected if SM and VPD were uncoupled. Therefore, identifying the mechanisms of the SM–VPD coupling, especially the bimodal SM–VPD distribution, is crucial for our understanding of the occurrence of the compound drought and aridity events. Fig. 1. Relationship between GLEAM root-zone SM and MERRA-2 VPD. (A) Correlation coefficient between root-zone SM and VPD during the period 1980–2017 (warm season in SI Appendix, Fig. S1). (B) Mean probability of each percentile bin of root-zone SM and VPD across all land grid cells.

LA Feedbacks in Concurrent Soil Drought and Atmospheric Aridity The GLACE-CMIP5 experiment consists of transient climate-change experiments from several climate models (Materials and Methods). In each model, a reference simulation (REF, similar to CMIP5 simulations) was performed with fully coupled SM responding to precipitation and evaporative demand. This reference simulation was then compared to an idealized simulation (expB) in which SM was prescribed, i.e., uncoupled from the atmosphere. The prescribed SM in the expB simulation preserves the long-term trend and seasonal cycle of SM in the REF simulation, but removes the subseasonal and interannual variability of SM in response to the atmosphere (see SI Appendix, Fig. S3 for an illustration of the methodology). Since the atmospheric model was driven by identical forcing (e.g., sea surface temperatures, sea ice, land use, and radiative forcing) in both simulations, the difference between the 2 simulations highlights the impact of SM variability on the atmosphere. We measure the SM–VPD correlation in the warm season from 1981 to 2080. Negative SM–VPD correlation is found globally in REF (Fig. 2A), which is consistent with the observational relationship (Fig. 1A). Further, concurrent SM and VPD extremes occur more frequently in regions with more negative SM–VPD correlations in REF (Fig. 2G). The bimodal distribution of SM and VPD is also evident in REF (Fig. 2E and SI Appendix, Fig. S4 A–C for each model), indicating that the models with fully coupled SM can realistically simulate the compound SM and VPD extremes. However, the SM–VPD bimodality is not found in expB (Fig. 2F and SI Appendix, Fig. S4 D–F). In expB, VPD variability is driven only by the atmosphere; without the feedback of SM variability to VPD, particularly under dry soil conditions, the frequency and intensity of extreme high VPD in expB are greatly reduced in most regions (Fig. 2 B and C). This result indicates that the feedback of SM on VPD is largely responsible for enabling aridity extremes. In addition, LA feedbacks enhance mean VPD in most regions (Fig. 2D), largely because the positive VPD response to low SM is larger than the negative VPD response to high SM (because of the exponential Clausius–Clapeyron relationship with temperature and lower sensitivity of evapotranspiration to SM at higher SM conditions). Fig. 2. SM-VPD relationship in GLACE-CMIP5 expB and REF simulations. (A) Correlation coefficient between SM and VPD [r(SM,VPD)] in REF during the period 1981–2080 (warm season). (B) Ratio of the frequency of extreme high VPD (above 95th percentile VPD in REF) between expB and REF. (C and D) Difference in (C) 95th percentile VPD (hPa) and (D) mean VPD (hPa) between REF and expB. (E and F) Mean probability of each percentile bin of VPD and SM across all grid cells in (E) REF and (F) expB. The percentile bins are sorted according to SM and VPD values in REF, and the bins are consistent in REF and expB. (G) Spatial relationship between r(SM,VPD) and the probability of concurrent VPD and SM extremes (above 95th percentile VPD and below 5th percentile SM) in REF simulations. We have demonstrated that large-scale atmospheric variability alone does not lead to extreme high VPD, which is enabled by local SM extremes. SM variability is induced by atmospheric dynamics, especially precipitation variability, but SM can also affect precipitation through its regulation of the water exchange and energy partitioning between the land surface and the atmosphere (21). Compared to REF, we find the frequency of extreme low precipitation (below 5th, 10th, and 15th percentiles of precipitation in REF) in expB is greatly reduced in most regions (SI Appendix, Fig. S5), especially in the Northern Hemisphere where VPD extremes are highly exacerbated by SM variability. This result indicates that the extreme low SM is not only attributed to atmospheric circulation anomalies, but also to local physical processes as SM–precipitation interactions also lead to more frequent extreme low SM and precipitation conditions (28). This is because low SM reduces evapotranspiration and increases sensible heat. Reduced evapotranspiration decreases moisture recycling and increased sensible heat, in some cases, may further enhance atmospheric stability and inhibit cloud formation and precipitation generation (21). These results highlight the importance of LA feedbacks in the occurrence of compound drought and aridity events and the key role of dry SM conditions in enhancing extreme VPD conditions and the mutual amplification of SM and precipitation deficits.

SM–VPD Correlation and Compound Drought and Aridity Events We further assess the frequency of compound drought and aridity events in CMIP5 models (listed in SI Appendix, Table S1) during the historical (1881–1980) and future (1981–2080) periods. The dependence of SM and VPD in the warm season is modeled with copulas (29) to derive the cooccurrence probability of extreme events. A compound drought and aridity event is defined as a month when mean VPD is above its 95th percentile and mean SM is below its 5th percentile. We use the probability multiplication factor (PMF) to measure the increase in the occurrence frequency of compound extreme events compared to the frequency if these extremes were assumed independent (P = 0.05 × 0.05 = 0.0025). PMF is defined as the ratio of the joint probability of extreme low SM and high VPD and the probability if they were assumed independent, and a PMF of 1 therefore represents no increase in the cooccurrence probability. Consistent with the GLACE-CMIP5 results, we find more frequent compound SM and VPD extremes (compared to if they were independent) in regions with more negative SM–VPD correlations in CMIP5 models (Fig. 3). The spatial patterns of the SM–VPD correlation and PMF are almost identical between the historical and future periods (Fig. 3 A–D). By comparison, negative SM–VPD correlation (PMF) becomes stronger (higher) in most mid- and high-latitude regions, but weaker (lower) in tropical regions in future simulations (SI Appendix, Fig. S6). To predict the dependence of PMF on the SM–VPD correlation, we use a linear relationship between r(SM,VPD) and the natural logarithm of PMF [exponentially increasing with r(SM,VPD)] across all grid cells of the land area (excluding Greenland and Antarctica in all analyses), with high explanatory power (R2 = 0.77 for historical simulations and R2 = 0.79 for future simulations) (Fig. 3 E and F). Comparing the functions between the 2 periods, we find future PMF becomes less sensitive to r(SM,VPD) (i.e., linear dependency between SM and VPD), but more strongly impacted by tail dependence in extreme deviations (SI Appendix, Eq. S7 and Table S2), which is modeled by copulas (30). The high explanatory power of the PMF-r(SM,VPD) function demonstrates that the frequency of compound extreme events can be well captured by the SM–VPD correlation in both periods. Fig. 3. SM-VPD correlation and probability multiplication factor (PMF) in CMIP5. (A and B) Mean correlation coefficient [r(SM,VPD)] between SM and VPD across 61 pairs of historical and future simulations. (C and D) Model mean PMF of concurrent extreme VPD above its 95th percentile and extreme SM below its 5th percentile. (E and F) Relationship between r(SM,VPD) and PMF across the land grids (excluding Greenland and Antarctica) for historical and future simulations. The SM–VPD correlation depends on the correlation between SM and temperature as well as that between SM and relative humidity. Interestingly, we find stronger SM–relative humidity correlation than SM–temperature correlation, and higher occurrence frequency of compound SM and relative humidity extremes than that of compound SM and temperature extremes, except in the northern high latitudes (SI Appendix, Fig. S7). These results are consistent with recent results (15) that the association between soil dryness and atmospheric dryness is tighter than that between soil dryness and atmospheric heating. This is because atmospheric dryness is enhanced not only by the increased sensible heat flux, but also the associated drying of the boundary layer induced by reduced evapotranspiration and increased boundary layer dry and warm air entrainment (31).

Changes in the Frequency and Intensity of Compound Drought and Aridity Events Due to LA Feedbacks We further remove the component of r(SM,VPD) due to the trends and seasonal cycles (warm season) of SM and VPD, R(SM t ,VPD t ), which may be related to the regional climatology and factors such as increased greenhouse gases. In this way, we isolate the contribution of the correlated subseasonal and interannual variations in SM and VPD, R(SM v ,VPD v ), induced by the atmospheric circulation dynamics and LA feedbacks (Fig. 4 A–F). R(SM v ,VPD v ) strongly depends on SM variability, i.e., the larger the SD of SM v , the stronger the R(SM v ,VPD v ) (SI Appendix, Fig. S8). This result further confirms the GLACE-CMIP5 result that SM variability is largely responsible for the SM–VPD feedback. Fig. 4. Attribution of SM-VPD correlation and PMF in CMIP5. Separated contributions from the subseasonal and interannual variations in SM and VPD [R(SM v ,VPD v )] and from the long-term trends and seasonality of SM and VPD [R(SM t ,VPD t )] to the SM–VPD correlation (A–F), and the contributions to the PMF (PMF v for variations and PMF t for trends and seasonality) (G–L). The last column shows the difference between historical (1881–1980) and future (1981–2080) periods (future minus historical values). Given the linear function established between r(SM,VPD) and the natural logarithm of PMF (SI Appendix, Table S2), we decompose PMF into the contribution from the correlated SM and VPD variations (PMF v ) and that from the long-term trends and seasonality of SM and VPD (PMF t ) (SI Appendix, Methods). PMF v measures the extent to which the compound extreme events are exacerbated by the SM–VPD feedback, with necessary impact of atmospheric circulation dynamics on SM extremes. The CMIP5 results show that the high probability of compound extreme events is mainly associated with the correlated SM and VPD variations [R(SM v ,VPD v )], although the trends and seasonality of SM and VPD [R(SM t ,VPD t )] also play a small part (Fig. 4 G–L). This result confirms a widespread SM–VPD feedback that promotes compound drought and aridity events. We find relatively low PMF v in some monsoonal regions, such as Brazil, West and Southeast Africa. In these regions, VPD and SM are strongly influenced by moisture transport from the ocean in the warm season (32). As a result, the contribution of seasonality is relatively large on the occurrence of compound SM and VPD extremes (SI Appendix, Fig. S9). Low PMF v also exists in some very dry regions such as the Sahara Desert and the Arabian Peninsula because SM is too low to induce variations in evapotranspiration and climatic variables. Compared to the historical period, future PMF v is slightly enhanced across 72% of the land area, although R(SM v ,VPD v ) is weaker in the future (Fig. 4 A–F). This is because of reduced dependence of PMF on the SM–VPD correlation, i.e., smaller sensitivity of PMF to r(SM,VPD) (Fig. 3 E and F), and stronger tail dependence in extreme deviations induced by large-scale atmospheric dynamics and LA feedbacks. PMF t also becomes larger across 68% of the land area, mainly due to the negative trend of warm-season SM and positive trend of VPD in most regions in the future period (SI Appendix, Fig. S9). The intensities of SM and VPD extremes are greatly enhanced, with much higher (lower) thresholds and mean values of extreme high VPD (low SM) in the future (SI Appendix, Fig. S10). After subtracting the changes in mean SM and VPD due to the long-term climate change, the thresholds and mean values of SM and VPD for compound extreme events still become more extreme (Fig. 5). Specifically, the threshold and mean values of extreme high VPD increase by 18 ± 9% and 28 ± 15% (relative to mean VPD in the historical period), respectively, with the highest increases located in the Amazon and Europe (Fig. 5 C and D). The increase of extreme high VPD occurs even in the regions where the extreme low SM increases (less intense SM drought) in the future. This is because, as future temperature rises, VPD is projected to become increasingly sensitive to SM variability, especially under low SM conditions. These results indicate that LA feedbacks will greatly exacerbate future extreme aridity beyond that induced by mean climate change alone. The combination of a slight increase in the frequency and large increases in the intensity of compound SM and VPD extremes induced by LA feedbacks indicates large risks of compound drought and aridity events in the future. Fig. 5. Changes in the thresholds and mean values for compound extreme events from historical to future periods in CMIP5. (A and C) Changes (future minus historical values) in thresholds of (A) extreme low SM and (C) extreme high VPD. (B and D) Changes in (B) mean SM and (D) mean VPD for compound extreme events. We have subtracted the differences in mean SM and VPD between the 2 periods from the changes in the thresholds and mean values of SM and VPD for compound extreme events, and the remaining values are normalized by mean SM and VPD in historical simulations, respectively.

Discussion Compound drought and aridity events can be attributable to both large-scale atmospheric dynamics and local LA feedbacks. Large-scale atmospheric circulation anomalies (e.g., a high-pressure system that blocks precipitation and promotes heat via enhanced sunshine and subsidence) may trigger climate extremes (13, 33, 34) and extreme low SM. For compound extreme events assessed at the monthly scale, our results show exacerbation by SM–climate feedbacks, as atmospheric circulation anomalies alone cannot generate substantial compound extreme events (Fig. 2). Reduced SM due to initial precipitation deficit can enhance the heating and drying of the near-surface atmosphere (16, 19), leading to high VPD values; further, it can decrease moisture recycling and the likelihood of cloud/precipitation generation in some cases (21), resulting in a positive feedback to SM itself and hence the exacerbation of compound soil drought and atmospheric aridity. We note that the impacts of atmospheric circulation dynamics and LA feedbacks on compound drought and aridity events are intertwined since SM variability induced by atmospheric circulation dynamics is necessary for the SM–climate feedback, and SM itself can also affect atmospheric circulation dynamics (17). Yet, the experimental simulations in GLACE-CMIP5 demonstrate the strong enhancement of compound drought and aridity events by LA feedbacks. The crucial role of LA feedbacks in promoting compound drought and aridity events has important implications for the occurrence of flash droughts and prolonged droughts, which occur at shorter or longer timescales than the compound extreme events (monthly scale) assessed in this study. The dependence of extreme high VPD on low SM promotes rapid drought intensification, i.e., flash droughts, which occur due to quick SM depletion and atmosphere heating and drying (35). Extreme high VPD may also precede and drive SM depletion to trigger flash droughts (35, 36). Declining vegetation health and its weakening regulation on water and energy exchanges also accelerate the development of flash droughts (37). The regions with strong LA feedbacks are also more prone to prolonged droughts: indeed, atmospheric dynamics on their own do not have much memory, but LA feedbacks induce memory that is likely important in promoting prolonged droughts, possibly including the infamous megadroughts known from proxy records to have occurred repeatedly over the past millennium in Europe (16), Asia (12), and North America (38). In addition, the heightened atmospheric aridity tends to drive more wildfires during soil droughts (39), especially in a future with strong aridity intensification. Models project LA feedbacks to increase the frequency and especially the intensity of compound drought and aridity events in the future (Figs. 4 and 5). The strength of LA feedbacks, measured by R(SM v ,VPD v ), is not projected to increase, probably because of the relatively stable magnitude of SM variability between historical and future periods (SI Appendix, Fig. S8). Thus, we find only a slight increase in the frequency of compound drought and aridity events enhanced by LA feedbacks (i.e., PMF v ) because of stronger tail dependence in extreme deviations in future simulations. We note the attribution analysis of PMF relies on the relationship between PMF and r(SM,VPD) (Fig. 3 E and F and SI Appendix, Table S2). Although r(SM,VPD) can well predict PMF in both periods, this relationship may vary slightly across space due to varying dependence structures (i.e., copulas) between SM and VPD, which may lead to some uncertainties in the attribution results. It should also be noted that the CMIP5 models simulate weaker SM–VPD correlations than those obtained from the observationally constrained data (Figs. 1 and 3). Therefore, the simulated frequencies of historical and future compound drought and aridity events may be underestimated. On the other hand, LA processes strongly enhance the responses of SM and VPD extremes to anthropogenic warming, leading to more intense compound extreme events. This is because higher temperature enhances both saturated water vapor and evaporative demand. Therefore, VPD will be more sensitive to soil drying and in turn amplify SM deficits. The increases in the frequency and especially intensity of compound drought and aridity events could have dramatic impacts on humans and natural systems (9). VPD and SM stresses can directly limit plant stomatal conductance and carbon uptake, and even trigger vegetation mortality (3, 7, 40). Future rising CO 2 may reduce plant stomatal conductance and evaporative water loss, slowing down the development of soil drought, but exacerbating atmospheric aridity (15). Thus, rising CO 2 may weaken the negative SM–VPD coupling, as shown in Fig. 4C. The CO 2 fertilization effect is expected to enhance photosynthetic carbon assimilation, but this effect may be largely offset by intensified compound drought and aridity events, increased respiration losses due to rising temperature, and vegetation mortality (9, 41). The strengthened risks of compound drought and aridity events may result in substantial decreases in agricultural production (6, 42) and continental carbon storage (9) if terrestrial ecosystems cannot adapt to future climate changes. This study highlights the importance of SM variability in enabling a series of processes and feedback loops affecting near-surface climate. There is thus a crucial need to better quantify and evaluate the representation of such processes in climate models. In particular, accurate model representation of both SM variability and associated feedbacks appears indispensable to provide reliable simulations of the frequency, duration, and intensity of compound drought and aridity events and of their changes in a warmer climate, and ultimately to mitigate future risks associated with these events.

Materials and Methods Datasets. We used GLEAM SM and MERRA-2 VPD in the warm season of 1980–2017 to assess the relationship between SM and VPD (SI Appendix, GLEAM Dataset and MERRA-2 Reanalysis). The experimental simulations (expB) with prescribed SM and the reference simulations (REF) with interactive SM in 3 models (i.e., ACCESS, ECHAM6, and GFDL) participating the GLACE-CMIP5 experiment were used to isolate the impacts of SM-atmosphere feedbacks on the occurrence of compound drought and aridity events (SI Appendix, GLACE-CMIP5 Experiment). We further used 26 CMIP5 models (SI Appendix, Table S1) to evaluate the contributions of LA feedbacks to the changes in the frequency and intensity of compound drought and aridity events between historical (1881–1980) and future (1981–2080, RCP8.5 scenario) simulations (SI Appendix, CMIP5 Model Simulations). PMF. The dependence structure of SM and VPD was modeled with bivariate copulas (29), which have been widely used to assess the relationship between dependent variables (43, 44). The copulas can overcome the shortcomings of counting the cooccurrence rate of extreme low SM and high VPD with few samples (44), especially in the regions where SM and VPD are weakly or not coupled. Here, we have the joint probability distribution of SM and VPD F S M , V P D ( x , y ) = P ( S M ≤ x , V P D ≤ y ) , [1] and the marginal cumulative distribution functions F S M ( x ) = P ( S M ≤ x ) and F V P D ( y ) = P ( V P D ≤ y ) . We use a bivariate copula C to describe the joint probability distribution of SM and VPD as F S M , V P D ( x , y ) = C [ F S M ( x ) , F V P D ( y ) ] = C ( u , v ) , [2] where F S M ( x ) and F V P D ( y ) are transformed into 2 uniformly distributed random variables u and v ranging from 0 to 1 (i.e., the normalized ranks of SM and VPD). As we define compound extreme events with VPD above its 95th percentile and SM below its 5th percentile, the joint probability of the compound extreme event is given by P ( u < 0.05 ∩ v > 0.95 ) = P ( u < 0.05 ) − P ( u < 0.05 ∩ v ≤ 0.95 ) = 0.05 − C ( 0.05 , 0.95 ) . [3] To model the joint probability distribution, we considered commonly used bivariate copula families (Gaussian copula, Student’s t copula, and Archimedean copulas) and tested their goodness of fit according to the Bayesian Information Criterion (45). In each grid cell, the copula with the best fit was selected to represent the joint probability distribution of SM and VPD. These analyses were performed using the R package, VineCopula (46). With the joint probability distribution of SM and VPD, we derived the occurrence probability of compound extreme events using Eq. 3. The PMF was defined as the ratio of the joint probability derived from the copula with the best fit and that assuming independent distributions, for the compound extreme events. We derive the joint probability of independent distributions from the independent copula: P = 0.05 − C ( 0.05 , 0.95 ) = 0.05 − 0.05 × 0.95 = 0.0025 . The independent copula is commonly used as a scale factor to transform a low joint probability to a ratio that measures the dependence strength of bivariate extremes (9, 43). We also counted the joint occurrence rate of extreme events based on model simulations and calculated the PMF in a similar method. We compared the PMF derived with copulas and by counting and found that the copulas could capture the occurrence frequency of compound extreme events well (SI Appendix, Fig. S11). We also calculated the PMF for the concurrent SM and temperature extremes, and the SM and relative humidity extremes using copulas.

Acknowledgments We acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups (listed in SI Appendix, Table S1 of this paper) for producing and making available their model output. For CMIP the US Department of Energy’s Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. S.Z. acknowledges support from the Lamont-Doherty Postdoctoral Fellowship and the Earth Institute Postdoctoral Research Program. P.G. acknowledges support from NASA ROSES Terrestrial hydrology (NNH17ZDA00IN-THP) and NOAA MAPP NA17OAR4310127. A.P.W. and B.I.C. acknowledge support from the NASA Modeling, Analysis, and Prediction (MAP) program (NASA 80NSSC17K0265) and Columbia University’s Center for Climate and Life. This work was also supported by the Lamont-Doherty Earth Observatory Contribution No. 8347.

Footnotes Author contributions: S.Z. designed research; S.Z., A.P.W., A.M.B., B.I.C., Y.Z., S.H., R.L., S.I.S., and P.G. performed research; S.Z. and Y.Z. analyzed data; A.P.W., A.M.B., B.I.C., Y.Z., S.H., R.L., S.I.S., and P.G. edited the paper; and S.Z. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission. P.A.D. is a guest editor invited by the Editorial Board.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1904955116/-/DCSupplemental.