Anthropogenic Forcing Decreased Concurrent Soil Drought and Atmospheric Aridity in the Historical Period 1850–2013

Abstract

As awareness of the importance of anthropogenic climate change has increased, attention is being focused on “compound extremes,” such as co-occurring soil drought (low amounts of soil moisture, SM) and atmospheric aridity (high vapor-pressure deficits, VPD), because of the disproportionate impacts of such extreme conditions on natural and societal systems. Few advances, however, have been made in isolating the net effect of anthropogenic forcing on the occurrence of such compound extremes. Nine Earth System models (ESMs) under natural-only and fully-forced simulations from Phase 6 of the Coupled Model Intercomparison Project (CMIP6) indicated that the presence of anthropogenic forcing in the historical period 1850–2013 weakened the coupling between SM and VPD and thus decreased the frequency of compound drought and aridity globally. The anthropogenically induced decrease in the strength of land–atmospheric feedbacks, which offset the drying trend effects of global warming by countering the expected natural correlation between SM and VPD, appears to have been responsible for the relative rarity of concurrent drought and aridity in the historical period. We also modeled independently the anthropogenic forcings of aerosols and greenhouse gases (GHGs) to further highlight that the widespread weakening of land–atmospheric feedbacks may have been associated primarily with the cooling effects induced by increases in anthropogenic emissions of aerosols, because the increase in intensity and frequency of compound drought and aridity that might have been expected from rising GHG concentrations was widely counter-balanced until very recently by an aerosol-driven cooling, particularly for the middle and high latitudes of the Northern Hemisphere. This finding indicates that the trade-off between these two main anthropogenic forcings may determine future patterns of concurrent drought and aridity in a changing climate. If global emissions of anthropogenic aerosols decrease in the future, as expected, our results imply a renewed strengthening of land–atmospheric feedbacks, and thus an intensification of concurrent drought and aridity.

Introduction

Drought exerts an increasing pressure on terrestrial ecosystems as we confront anthropogenic climate change (Humphrey et al., 2021; Reichstein et al., 2013; M. Zhao & Running, 2010). A low amount of soil moisture (SM, an indicator of soil drought) and a high level of atmospheric vapor-pressure deficit (VPD, an indicator of atmospheric aridity) are two main drivers of moisture stress on vegetation that can substantially reduce both agricultural production (Rigden et al., 2020) and the uptake of terrestrial carbon by non-agricultural plants (Green et al., 2019; Yuan et al., 2019), and can even drive widespread increases in tree mortality (Anderegg et al., 2013). In response to high VPD, plants are inclined to increase their stomatal closure to minimize water loss; decreased SM can further reduce stomatal conductance and impair hydraulic transfer from soil to leaves (Grossiord et al., 2020; Sperry et al., 2002). Recent studies have tried to separate the individual roles of SM and VPD in controlling vegetation productivity but with conflicting results (Fu et al., 2022; Liu et al., 2020; Yuan et al., 2019), mainly because low SM and high VPD typically co-occur (S. Zhou, Williams, et al., 2019). During SM drought, low amounts of SM limit evapotranspiration and thus reduce evaporative cooling, causing higher surface temperature and VPD (Bateni & Entekhabi, 2012). The enhanced VPD, however, in turn exacerbates SM depletion by increasing atmospheric water demand. Feedback of this sort between SM and VPD hinders the separation of their relative impacts; individually, drought and aridity would each be expected to have far weaker negative impact on vegetation productivity than when extreme high VPD and extreme low SM occur simultaneously (S. Zhou, Zhang, et al., 2019). Understanding mechanisms that would tend to increase or decrease concurrent SM and VPD extremes is thus of particular importance for managing and minimizing the emerging environmental risk of compound drought and aridity associated with anthropogenic climate change. The net effects of anthropogenic forcing on the coupling between SM and VPD, especially its implications for compound drought and aridity, however, remain largely unknown.

Previous studies have indicated that the occurrence rate of compound drought and aridity is expected to increase in the 21st century compared to the historical period, mainly due to strengthened land–atmospheric feedbacks (S. Zhou, Williams, et al., 2019; S. Zhou, Zhang, et al., 2019). However, anthropogenic forcings via greenhouse gases (GHGs) and aerosols can be expected to differ between the two periods, as the global mean emissions of anthropogenic GHGs and aerosols (e.g., sulfate, organic carbon, and black carbon) increased synchronously in the historical period, but the opposite trend has begun to occur in the 21st century. That is, GHG emissions are very likely to continue to increase while anthropogenic aerosol (AER) emissions can be expected to decrease due to the success of socioeconomic development and the pursuit of a cleaner atmosphere, both of which will contribute to continued global warming (S. Zhao et al., 2017). Comparison of the two periods thus not only highlights the role of warming, but suggests that it will be difficult to determine whether decreasing aerosol emissions or increasing GHG emissions will dominate in determining the rate of increase in instances of simultaneous soil drought and atmospheric aridity.

It is worth noting that although the Earth's global mean temperature over land has increased by 1.53°C since the preindustrial period, this rise would most probably have been even greater if GHGs-induced increase in temperature had not been largely offset by aerosol-driven cooling (Shukla et al., 2019). Does this mean that recent changes in anthropogenic forcing have already strengthened the feedback between low SM and high VPD (i.e., by promoting global warming), and thus we should expect increased occurrences of concurrent drought and aridity, compared to their frequency in the last several centuries?

To address these questions, we first examined the strong SM–VPD coupling under the natural-only (volcanic and solar radiative) forced (NAT-only) and the fully (natural and anthropogenic) forced (Fully-forced) simulations resulting from the nine climatic models in Phase 6 of the Coupled Model Intercomparison Project (CMIP6); we then examined differences among the nine models and characterized the net impacts of anthropogenic forcing on their correlations between SM and VPD, and their estimates of the occurrences of compound drought and aridity in the historical period 1850–2013. Next, to investigate how anthropogenic forcings influence the frequency of concurrent drought and aridity, we established an attribution model for SM–VPD coupling, mainly by altering the long-term trend of SM and VPD or the strength of land–atmospheric feedbacks (see Section 2). Finally, we quantified the intensity and frequency of concurrent SM and VPD extremes in the GHG-only and aerosol-only simulations and analyzed the individual impacts of the GHG and AER forcings on such compound extremes. This investigation of the role of various anthropogenic forcings in driving changes in concurrent SM and VPD extremes has led to results that can serve to broaden our understanding of how to warning, monitor, and perhaps respond to important consequences of anthropogenic climate change.

Materials and Methods Data CMIP6 Data

We used monthly data for surface and total SM, relative humidity and temperature from nine Earth System models (ESMs) (Table S1 in Supporting Information S1) participating in CMIP6. These models included the all-forcings (Fully-forced), natural-only (NAT-only), aerosol-only (AER-only) and greenhouse-gas-only (GHG-only) simulations in the historical period 1850–2014. The Fully-forced simulation included all prominently anthropogenic forcings (emissions of GHGs and AERs) and natural forcings (volcanic and solar radiative), and the NAT-only, AER-only, and GHG-only simulations represented the historical period with natural forcings only, AER emissions only and GHG emissions only, respectively. The differences between the Fully-forced (GHG-only and AER-only) and the NAT-only simulations allowed us to understand the impacts of some major anthropogenic forcings (individual forcings) on the dependence between low SM and high VPD. Vapor-pressure deficit was calculated from the relative humidity and temperature (Yuan et al., 2019).

All analyses concentrated on the warm season, which was defined as the hottest three months for each grid from 1850 to 1979 in the NAT-only simulation (Figure S1 in Supporting Information S1). The warm season corresponded to the period when loss of terrestrial carbon induced by drought and aridity mostly occurred (S. Zhou, Zhang, et al., 2019), and its duration of 3 months was a compromise between potentially longer seasonal droughts and aridity, which usually occur on shorter time scales (Zscheischler & Seneviratne, 2017). All model outputs were interpolated to 1° × 1° grid using a bilinear function for terrestrial areas to create multi-model ensembles of each scenario.

Observation Data

Monthly temperature, relative humidity, and cloud cover data with a spatial resolution of 0.5° × 0.5° were obtained from Climate Research Unit at the University of East Anglie (CRU TS 4.06) to produce the observation-based changes in VPD during 1900–2013 (Figure 1b) and to verified the global cycle of “dimming and brightening” (Figure S2 in Supporting Information S1). To verify the effectiveness of CMIP6 in simulating SM, monthly SM data with a spatial resolution of 0.25° × 0.25° from Global Land Evaporation Amsterdam Model V3.5a (GLEAM) were used to correlate with the SM output from CMIP6 models during their common period (i.e., 1980–2013; they are positively correlated globally as shown in Figure S3 in Supporting Information S1). All these data sets were bilinearly interpolated to a spatial resolution of 1° × 1° prior to analysis.

Figure 1. Anomalies of global mean (a) soil moisture (SM) and (b) vapor-pressure deficit (VPD) in the fully forced (Fully-forced) and the natural-only forced (NAT-only) simulations in the historical period. Bold lines and shaded areas illustrate the mean and standard deviation of VPD and SM simulated by nine Phase 6 of the Coupled Model Intercomparison Project models. The blue line in (b) exhibited the changes in VPD anomalies calculated from CRU TS observation data.

Methods Joint Probability of Drought and Aridity Using Copulas

The dependence between SM and VPD in our study was formulated by the bivariate copulas (Nelsen, 2006), which are extensively used as an effective statistical method to describe the dependence between random variables and to derive joint probability of a compound extreme (Wahl et al., 2015; S. Zhou, Zhang, et al., 2019). Mathematically, assuming two random variables, A (e.g., SM) and B (e.g., VPD), with marginal cumulative distribution functions $${F}_{A}(a)={P}_{\mathrm{r}}(A\le a)$$ and $${F}_{B}(b)={P}_{\mathrm{r}}(B\le b)$$, respectively, their joint probability distribution can be obtained by the copula ($$C$$): [Image Omitted. See PDF]where $$F(a,b)$$ is the joint distribution function of $$A$$ and $$B$$ (Salvadori & De Michele, 2004). The two marginal distribution functions $${F}_{A}(a)={P}_{\mathrm{r}}(A\le a)\,\mathrm{a}\mathrm{n}\mathrm{d}\,{F}_{B}(b)={P}_{\mathrm{r}}(B\le b)$$ were converted into two uniformly distributed random variables $$W$$ and $$Z$$ that ranges from 0 to 1, so the joint probability distribution of A and B could be described as: [Image Omitted. See PDF]

The probability of an extreme compound event, where the variables were below or exceeded a given threshold, is thus expressed as: [Image Omitted. See PDF]

In our study, the concurrent soil drought and atmospheric aridity event was defined as SM below its 10th percentile and VPD above its 90th percentile, and the probability of such co-occurrence extremes could thus be calculated as: [Image Omitted. See PDF]

Model Fitting and Calculation of LMF

We obtained the uniform distributions of SM and VPD by first converting the marginal distributions of SM and VPD to a normalized rank that ranges from 0 to 1, and this is a general operation when we use copulas (Zscheischler & Seneviratne, 2017). We then fitted all possible families of bivariate copulas (40 in all and listed in Table S2 in Supporting Information S1) and chose the one of best fit in each pixel using the Bayesian Information Criterion implemented in the R package VineCopula (Schepsmeier et al., 2017). We used the selected copula (i.e., the best fit one) to calculate the probability of occurrence of compound extremes using Equation 4 and “BiCopCDF” function in VineCopula.

If SM and VPD were independent, the joint probability of an extreme low SM (below its 10th percentile) and a high VPD (above its 90th percentile) would be 0.1 × 0.1 = 0.01. We therefore defined the likelihood multiplication factor (LMF) as the ratio of the joint probability calculated by copula to that assuming they were independent (i.e., 0.01) to measure the strength of dependence of bivariate extremes (S. Zhou, Williams, et al., 2019). An LMF >1 thus indicates an increase in the occurrences rate of concurrent drought and aridity relative to the frequency if SM and VPD extremes were assumed independent. We also calculated the joint probability by directly counting concurrent SM and VPD extreme events to assess the robustness of the results retrieved by the copula. By comparation, LMF and its changes between Fully-forced and NAT-only derived by the copula or by counting were highly consistent (Figure S4 in Supporting Information S1), suggesting that the copulas could well identify the occurrence rate of compound SM and VPD extremes.

Attribution Method for SM–VPD Correlations

We determined whether anthropogenic forcing influenced the SM–VPD coupling mainly by altering their long-term trends or the strength of land–atmospheric feedbacks. Following the S. Zhou, Williams, et al. (2019) and we decomposed the time series of SM and VPD in the warm season into long-term trends and seasonality (i.e., $${\mathrm{S}\mathrm{M}}_{\mathrm{t}}$$ and $${\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{t}}$$) and the sub-seasonal and inter-annual variations (i.e., $${\mathrm{S}\mathrm{M}}_{\mathrm{v}}$$ and $${\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}$$). $${\mathrm{S}\mathrm{M}}_{\mathrm{t}}$$ and $${\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{t}}$$ were calculated as their 30-year running monthly means, and $${\mathrm{S}\mathrm{M}}_{\mathrm{v}}$$ and $${\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}$$ were defined as the difference between their monthly values and the 30-year running means for each month in the historical period. Specifically, time series of monthly SM and VPD could be expressed as: [Image Omitted. See PDF] [Image Omitted. See PDF]

The correlation between SM and VPD could therefore be written as: [Image Omitted. See PDF]where $$\mathrm{c}\mathrm{o}\mathrm{v}()$$ represents the covariance between two variables and $$\sigma $$ is the standard deviation. The relationship between $${\mathrm{S}\mathrm{M}}_{\mathrm{t}}$$ and $${\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}$$ and between $${\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{t}}$$ and $${\mathrm{S}\mathrm{M}}_{\mathrm{v}}$$ were independent, so the values of $$\frac{\mathrm{c}\mathrm{o}\mathrm{v}\left({\mathrm{S}\mathrm{M}}_{\mathrm{t}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}\right)}{{\sigma }_{\mathrm{S}\mathrm{M}}\cdot {\sigma }_{\mathrm{V}\mathrm{P}\mathrm{D}}}$$ and $$\frac{\mathrm{c}\mathrm{o}\mathrm{v}\left({\mathrm{S}\mathrm{M}}_{\mathrm{v}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{t}}\right)}{{\sigma }_{\mathrm{S}\mathrm{M}}\cdot {\sigma }_{\mathrm{V}\mathrm{P}\mathrm{D}}}$$ were very small and could thus be omitted in the calculation (S. Zhou, Williams, et al., 2019). $$r(\mathrm{S}\mathrm{M},\mathrm{V}\mathrm{P}\mathrm{D})$$ could therefore be abbreviated as: [Image Omitted. See PDF]where $$R()$$ is the normalized correlation coefficient, for example, $$R\left({\mathrm{S}\mathrm{M}}_{\mathrm{v}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}\right)$$ = $$r\left({\mathrm{S}\mathrm{M}}_{\mathrm{v}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}\right)\times \frac{{\sigma }_{{\mathrm{S}\mathrm{M}}_{\mathrm{v}}}\times {\sigma }_{{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}}}{{\sigma }_{\mathrm{S}\mathrm{M}}\cdot {\sigma }_{\mathrm{V}\mathrm{P}\mathrm{D}}}$$.

Attribution Method for Compound Soil Drought and Atmospheric Aridity

$$r(\mathrm{S}\mathrm{M},\mathrm{V}\mathrm{P}\mathrm{D})$$ and the logarithm of LMF were linearly correlated across all pixels in both the Fully-forced and NAT-only simulations (Figures 2g and 2h): [Image Omitted. See PDF]where $$a$$ and $$b$$ represent the slope and intercept, respectively. The intercept $$b$$ measures the impact of tail (in)dependence on $$\mathrm{L}\mathrm{M}\mathrm{F}$$. Based on a linear relationship established between $$r(\mathrm{S}\mathrm{M},\mathrm{V}\mathrm{P}\mathrm{D})$$ and $$\mathrm{L}\mathrm{M}\mathrm{F}$$, $$\mathrm{L}\mathrm{M}\mathrm{F}$$ could be further separated into $${\mathrm{L}\mathrm{P}\mathrm{F}}_{\mathrm{t}}$$ and $${\mathrm{L}\mathrm{P}\mathrm{F}}_{\mathrm{v}}$$ by substituting Equation 8 in Equation 9. We could then obtain: [Image Omitted. See PDF]and thus [Image Omitted. See PDF] [Image Omitted. See PDF]where $${\mathrm{L}\mathrm{M}\mathrm{F}}_{\mathrm{t}}$$ represents the changes in $$\mathrm{L}\mathrm{M}\mathrm{F}$$ induced by the long-term trends and seasonalities of SM and VPD, and $${\mathrm{L}\mathrm{P}\mathrm{F}}_{\mathrm{v}}$$ represents the changes in $$\mathrm{L}\mathrm{M}\mathrm{F}$$ induced by large-scale atmospheric dynamics and the SM–VPD feedback.

Figure 2. Model mean SM–VPD correlations (r(SM, VPD); a, b) and the likelihood of multiplication factor (LMF; d, e) in the fully forced (Fully-forced; a, d) and the natural-only forced (NAT-only; b, e) simulations in the historical period. Shifts in mean r(SM, VPD) (c) and likelihood multiplication factor (LMF) (f) between the Fully-forced and the NAT-only simulations (Fully-forced minus NAT-only). Relationship between r(SM, VPD) and LMP in the Fully-forced (g) and the NAT-only (h) simulations across terrestrial grids. (i) Relationship between shifts in r(SM, VPD) and shifts in LMF across terrestrial grids.

Individual Impacts of the GHGs and Aerosol Forcings on SM–VPD Coupling

Based on the same process, we also decomposed the SM and VPD series into $${\mathrm{S}\mathrm{M}}_{\mathrm{t}}$$, $${\mathrm{S}\mathrm{M}}_{\mathrm{v}}$$, $${\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{t}}$$, and $${\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}$$ and calculated $$R\left({\mathrm{S}\mathrm{M}}_{\mathrm{v}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}\right)$$ in both GHG-only and AER-only simulations. The difference of $$R\left({\mathrm{S}\mathrm{M}}_{\mathrm{v}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}\right)$$ between the GHG-only simulation and NAT-only simulation thus mainly measures the changes in the contribution of land–atmospheric feedbacks to SM–VPD couplings due to the presence of GHGs emissions. Similarly, the difference of $$\mathit{R}\left({\mathrm{S}\mathrm{M}}_{\mathrm{v}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}\right)$$ between the AER-only simulation and NAT-only simulation mainly measures the changes in the contribution of land–atmospheric feedbacks to SM–VPD couplings due to the presence of AER emissions. In addition, we also compared the thresholds of extreme low SM and extreme high VPD in the GHG-only and NAT-only simulations to the NAT-only simulation to analyze the individual impacts of the GHGs and aerosol emissions on the intensities of soil drought and atmospheric aridity.

Results Forced Warm-Season Trends in SM and VPD in the Historical Period

Figure 1 presents long-term trends of warm-season (the 3 months with highest mean temperature; Figure S1 in Supporting Information S1) global mean terrestrial SM and VPD in the historical fully forced (Fully-forced) and natural-only (NAT-only) simulations. SM was less influenced by anthropogenic forcing, showing no apparent trend during the whole study period in the Fully-forced simulation (Figure 1a). Similarly, global mean VPD had no obvious major trends before the 1990s but increased sharply in recent decades, and this trend was highly consistent with that calculated from CRU TS observations (Figure 1b). Compared to the NAT-only simulation, VPD decreased abnormally from the 1950s and later recovered in the last two decades in the Fully-forced simulation, coinciding with a global cycle of “dimming and brightening” (Wild, 2012). Global dimming refers to the tendency, from 1960 to 1990, of the scattering effect of aerosols to reduce solar radiation at the Earth's surface; since 1990, however, primarily in the Northern Hemisphere, this scattering effect diminished, a phenomenon referred to as global brightening (Wild et al., 2005; Figure S2 in Supporting Information S1). Due to the simultaneous increase in relative humidity during the period of global dimming (Figure S5 in Supporting Information S1), global mean VPD clearly did not follow the rapid increase in temperature, especially in most regions of Northern Hemisphere, such as Europe, central North America, and northwestern Asia (Figure S6 in Supporting Information S1). Abnormal increases in VPD coincided particularly well with abnormal decreases in SM in the NAT-only simulation, but that coincidences were weaker in the Fully-forced simulation, especially during the global dimming period (Figure 1).

Impacts of Anthropogenic Forcing on SM–VPD Coupling

We specifically explored the impacts of anthropogenic forcing on SM–VPD coupling by determining the correlation between SM and VPD and modeling their warm season dependence using copulas in the Fully-forced and NAT-only simulations during the historical period 1850–2013 (Figure 2). SM and VPD were negatively correlated globally in both simulations (Figures 2a–2b), confirming the expectation that higher VPD is generally accompanied by lower SM in warm seasons. The spatial distribution of the correlations between SM and VPD were almost identical between the Fully-forced and NAT-only simulations, both exhibiting strong negative correlations (i.e., correlation coefficients <−0.8) mainly in southwestern North America, south-central Africa, southern Asia, and much of South America and Australia. By comparison, the presence of anthropogenic forcing weakens the correlation between SM and VPD in most (85% of) terrestrial areas, except for a few regions at high latitudes and in southwestern Asia (Figure 2c). This widespread weakening of the SM–VPD correlation is almost identical in different historical periods (Figure S7 in Supporting Information S1); even in recent decades (i.e., 1985–2004), when global mean VPD exhibited an apparent increasing trend, SM–VPD correlations were still weaker in most regions under the Fully-forced simulation relative to the NAT-only. A comparison of the differences of the SM–VPD correlations between 1985 and 2004 and previous decades in the Fully-forced simulation (1985–2004 minus previous decades), however, indicates lower absolute correlation coefficient values during 1985–2004 in most terrestrial areas, especially in the Northern Hemisphere (Figure S8 in Supporting Information S1), suggesting that the tendency of anthropogenic forcing to weaken the dependence between low SM and high VPD has diminished in recent decades.

The spatial patterns of the likelihood multiplication factor (LMF, an LMF >1 indicated an increase in the cooccurrence probability of SM and VPD extremes and vice versa, see Section 2) are highly consistent with the distribution of the SM–VPD correlations, exhibiting high LMF mainly in regions with the more highly-negative correlations between SM and VPD (Figures 2a–2f). This finding suggested that SM–VPD correlations may be a good indicator of the likelihood of concurrent SM and VPD extremes in both the Fully-forced and NAT-only simulations.

We verified this possibility by identifying a linear relationship between the SM–VPD correlation coefficients and the natural logarithm of LMF over all the pixels of the terrestrial area (Figures 2g and 2h). The results indicate that the correlation between SM and VPD provides high explanatory power to LMF in both the NAT-only (R² = 0.83) and Fully-forced (R² = 0.81) simulations. Contrary to the general expectation that human activities have increased the likelihood of co-occurrence of SM and VPD extremes, we found that the frequency of such compound extremes decreased almost globally in the historical period in the Fully-forced simulation compared to the NAT-only simulation (Figure 2f). This anthropogenically induced decrease in LMF could also be especially well signaled by weak SM–VPD correlations (R² = 0.30, Figure 2i). To confirm the robustness of our results, we also calculated LMF directly by counting the rate of the joint occurrence of extreme events; the results were highly consistent with those derived using copulas (Figure S4 in Supporting Information S1).

To further clarify the pathway of the effect of anthropogenic forcing on the concurrent SM and VPD extremes, we also investigated the SM–T and SM–RH couplings in the NAT-only and Fully-forced simulations and their changes caused by anthropogenic forcing (Figure S9 in Supporting Information S1). Results indicated that SM–RH correlations were extensively stronger than SM–T correlations in both the NAT-only and Fully-forced simulations, which is consistent with the findings of Berg et al. (2016), suggesting that the association between soil dryness and atmospheric aridity is tighter than that between soil dryness and atmospheric heating. But by comparison, we interestingly found that the SM–RH correlations and the LMF of concurrent SM and RH extremes were almost globally decreased in the Fully-forced simulation compared to the NAT-only simulation, while the SM–T correlations and the LMF of concurrent SM and T extremes were widely increased. This indicated that extreme high temperatures typically lead to SM depletion in the Fully-forced simulation, but not necessarily imply extreme atmospheric aridity due to the concomitant increase in RH, thus causing a net decoupling between SM drought and atmospheric aridity.

An assessment of the impacts of anthropogenic forcing on the intensities of SM and VPD extremes indicated that in the Fully-forced simulation compared to the NAT-only simulation the thresholds of atmospheric aridity (i.e., extreme high VPD) were much higher in the Southern Hemisphere, especially for arid regions, but were lower in middle and high latitudes of the Northern Hemisphere (Figures 3b and 3c). The thresholds of soil drought (i.e., extreme low SM) were also much higher in central North America and at high Eurasian latitudes, with slight increases in middle and southern South America, southern Africa, Australia, and southern Asia, suggesting decreased intensities of SM extremes for these regions in the Fully-forced simulation (Figures 3a and 3c). This widespread alleviation of SM droughts and the intensification of atmospheric aridity further confirmed the weakened negative SM–VPD coupling in the Fully-forced simulation. In particular, the probabilities of co-occurring soil drought and atmospheric aridity were greatly reduced in southern Africa, eastern Australia, and most regions between 30°N and 60°N when the NAT-only thresholds were used to define soil drought and atmospheric aridity (i.e., extreme low SM and high VPD) in the Fully-forced simulation (Figures 3d–3f). These findings indicated that, in the historical period, anthropogenic forcing both decreased the frequency of co-occurring soil drought and atmospheric aridity and reduced the intensity of such compound extremes in most regions, especially at mid-high latitudes of the Northern Hemisphere (Figure 3f).

Figure 3. Changes in the thresholds of mean extreme low soil moisture (SM) (a) and extreme high vapor-pressure deficit (VPD) (b) between the fully forced (Fully-forced) and the natural-only forced (NAT-only) simulations (Fully-forced minus NAT-only). (c) Latitudinal comparison of mean changes in the thresholds of extreme low SM (red) and extreme high VPD (blue). (d) Model mean likelihood of multiplication factor (likelihood multiplication factor, LMF) in the Fully-forced simulation, but the thresholds were based on the NAT-only simulation. (e) Changes in LMF between the Fully-forced and the NAT-only conditions (Fully-forced minus NAT-only); thresholds used to define extreme high VPD and low SM are both based on the NAT-only condition. (f) Latitudinal comparison of mean changes in LMF between the Fully-forced and the NAT-only simulations, and the thresholds used to define extreme high VPD and low SM are both based on the NAT-only condition.

Mechanisms of Anthropogenically Forced Changes in Co-Occurrence of SM and VPD Extremes

LMF_t and LMF_v in the two simulations based on the attribution model established with $$R\left({\mathrm{S}\mathrm{M}}_{\mathrm{t}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{t}}\right)$$ and $$R\left({\mathrm{S}\mathrm{M}}_{\mathrm{v}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}\right)$$ (see Section 2) indicated that the SM–VPD correlations and the high likelihood of co-occurrence of SM and VPD extremes were tightly associated with $$R\left({\mathrm{S}\mathrm{M}}_{\mathrm{v}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}\right)$$ and LMF_v, respectively (Figures 4g–4l), suggesting that the frequency of compound extremes was controlled mainly by land–atmospheric feedbacks in both the NAT-only and Fully-forced simulations. Interestingly, we found by comparation that anthropogenic forcing decreased the $$R\left({\mathrm{S}\mathrm{M}}_{\mathrm{v}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}\right)$$ and LMF_v almost globally in the historical period (Figures 4f and 4l). Although the anthropogenically induced atmospheric drying trend contributed to a larger LMF_t across more than half of the terrestrial area (Figures 4a–4f), the effects of such a trend were offset by a lower LMF_v (Figure 4l). These findings collectively suggested that anthropogenic forcing reduced the correlation between SM and VPD and the frequency of concurrent soil drought and atmospheric aridity in the historical period, mainly by weakening land–atmospheric feedbacks.

Figure 4. Contributions of the variations in soil moisture (SM) and vapor-pressure deficit (VPD) (RSMv,VPDv $R\left({\mathrm{S}\mathrm{M}}_{\mathrm{v}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}\right)$) and the long-term trends of SM and VPD (RSMt,VPDt $R\left({\mathrm{S}\mathrm{M}}_{\mathrm{t}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{t}}\right)$) to the SM–VPD correlation (a–f), as well as to the likelihood multiplication factor (LMFv for variations and LMFt for long-term trends) (g–l) in the fully forced (Fully-forced; a, d, g and j) and the natural-only (NAT-only; b, e, h and k) simulations. The last column exhibits the difference between the Fully-forced and the NAT-only simulations.

Examining the individual impacts of AER and GHG forcings on $$R\left({\mathrm{S}\mathrm{M}}_{\mathrm{v}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}\right)$$ indicated that the reduced contribution of land–atmospheric feedbacks to SM–VPD coupling may have been associated most closely with the increased aerosol emissions in the historical period: compared to the NAT-only simulation, $$R\left({\mathrm{S}\mathrm{M}}_{\mathrm{v}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}\right)$$ were lower in the AER-only simulation in 62% of the terrestrial area, especially in the Northern Hemisphere (Figure 5a). Interestingly, the presence of GHG forcing also slightly weakened $$R\left({\mathrm{S}\mathrm{M}}_{\mathrm{v}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}\right)$$ in 56% of the terrestrial areas, primarily in extreme dry regions (e.g., the Sahara Desert) and in vegetated regions (Figure 5b). Perhaps SM was too low to notably alter surface energy balance and induce perturbations in climatic variables in these extreme dry regions. But for the vegetated regions, GHG-induced weakened $$R\left({\mathrm{S}\mathrm{M}}_{\mathrm{v}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}\right)$$ may have been associated with the water-saving effects of vegetation in response to continuously increased CO₂ levels, which can reduce stomatal conductance and thus the evaporative loss of water, alleviating the development of soil drought but simultaneously intensifying atmospheric aridity, thereby weakening the land–atmospheric feedbacks (S. Zhou, Williams, et al., 2019; S. Zhou, Zhang, et al., 2019).

Figure 5. Changes in RSMv,VPDv $R\left({\mathrm{S}\mathrm{M}}_{\mathrm{v}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{v}}\right)$ (a, b), thresholds of extreme low soil moisture (SM) (c, d) and extreme high vapor-pressure deficit (VPD) (e, f) and the likelihood of multiplication (likelihood multiplication factor, LMF) of concurrent SM and VPD extremes (g, h) under the greenhouse-gas (GHG-only) and aerosol (AER-only) forcings compared to the natural-only (NAT-only) simulation (GHG-only or AER-only minus NAT-only values). Thresholds used to define extreme low SM and high VPD are both based on the NAT-only condition when calculating LMF (g, h).

The GHG-induced water-saving effect of vegetation, however, is negligible compared to its warming and drying effects. A comparison of the GHG-only and historical NAT-only climatic conditions found that the presence of GHG emissions greatly increased the thresholds of VPD extremes globally, especially in southeastern North America, the Amazon Basin, southern Africa, and the extremely arid regions (e.g., Sahara Desert and the Arabian Peninsula; Figure 5f). As for thresholds of SM extremes, in the GHG-only simulation they were much lower in the Amazon Basin and the high latitudes of Eurasia than in the NAT-only simulations (Figure 5d). The forced response to aerosols, however, produced an opposite global pattern of changes in SM and VPD threshold extremes (Figures 5c and 5e). Using the NAT-only thresholds to define soil drought and atmospheric aridity (i.e., extreme low SM and high VPD) in the AER-only and GHG-only simulations indicated that the larger decrease in LMF under the AER-only conditions occurred mainly in the regions where LMF was much higher in the GHG-only simulation (Figures 5g and 5h). This finding suggests that the trade-off between the aerosol and GHG forcings likely determines the strength of the negative coupling between SM and VPD and thus the patterns of concurrent drought and aridity in the historical period (Figure 2f).

Discussion

Increased VPD from rising temperatures due to anthropogenic warming has likely enhanced the atmospheric evaporation demand and thus accelerated the depletion of SM (Teuling et al., 2013), which in turn has increased temperatures of the terrestrial surface and instances of VPD as a result of increased sensible heat flux and reduced near-surface humidity (Hirschi et al., 2011). From this perspective, the frequency of compound SM and VPD extremes should, in theory, have increased with anthropogenic warming. Interestingly, the NAT-only and Fully-forced simulations in the CMIP6 models, however, indicate that the presence of anthropogenic forcing decreased the correlation between SM and VPD and the frequency of concurrent drought and aridity in the historical period (Figure 2), due primarily to anthropogenically weakened land–atmospheric feedbacks (Figure 4). This unexpected phenomenon may be tightly associated with the continuously increasing anthropogenic emissions of aerosols over the historical period considered here (Figure 5), which by multiple pathways can have altered land–atmospheric feedbacks and thus altered SM–VPD coupling.

First, AERs can affect the Earth's radiative balance (i.e., exert a radiative cooling of about −0.8 W/m² globally) by directly absorbing and scattering solar radiation and indirectly by acting as ice nuclei and cloud condensation, thus increasing cloud albedo (Mahowald, 2011). These cooling effects have tended to offset a large portion of the climatic forcing and warming induced by CO₂ and other GHGs in some areas with high aerosol emissions, leading to regional decreases in evaporation, precipitation, and VPD, and thus weakening land–atmospheric feedbacks (S. Zhao et al., 2017; H. Zhou et al., 2020).

Second, increased AERs can produce more but smaller droplets as it act as seeds for almost all formation of liquid water in the atmosphere (Rosenfeld et al., 2008). These smaller droplets, however, are less likely to merge into raindrops, and the lifetime of clouds can thus be extended (Westervelt et al., 2015). This possibility suggests that even if abnormally high temperatures reduce SM by increasing evaporation, the evaporated SM may persist in the atmosphere for longer time than it would have in the absence of the aerosols, leading to a pattern of concurrent low SM, high relative humidity, and low VPD, thus weakening the historical correlation between low SM and high VPD (Figure S9 in Supporting Information S1). This mechanism was also verified by our result that the co-occurrence of extreme high temperature and relative humidity in the historical period is substantially more frequent in the Fully-forced simulation than in the NAT-only simulation (Figure S10 in Supporting Information S1).

Third, historical changes in aerosol emissions have been neither linear nor spatially homogeneous (Undorf et al., 2018). For example, aerosol emissions in Europe and North America rose steadily from pre-industrial times until their peak around the 1970s; however it was only in the 1950s that economic development caused a rapid increase in Asian aerosol emissions, which have continued to increase since then (Hoesly et al., 2018). Such spatial differences in aerosol emissions can affect vertical and horizontal gradients of temperature, which in turn alter atmospheric circulation and thus SM–VPD coupling (Chen et al., 2018; Dong et al., 2022; Shen & Ming, 2018). In short, we find that the effects of aerosol-induced cooling may be the major explanation for the widespread weakening of land–atmospheric feedbacks in the historical period; that weakening tended to offset the expected tendency of GHG-induced warming to increase the frequency of co-occurrence soil drought and atmospheric aridity (Figure 5).

Note, however, that radiative forcing is likely to become strongly positive under scenarios of future twenty-first century climatic mitigation; aerosol emissions are expected to decrease continuously before the atmospheric CO₂ level stabilizes (Xu et al., 2018). The severe GHG-driven increase in dependence between low SM and high VPD may thus no longer be offset by the effects of aerosol-induced cooling, leading to a renewed strengthening of land–atmospheric feedbacks and thus to an intensification of concurrent drought and aridity. This may also partially account for the waning tendency, already observed in recent decades, of the effects of anthropogenic forcing to weaken the SM–VPD coupling (Figure S8 in Supporting Information S1); this return to stronger SM–VPD coupling has come as the trend of increasing aerosol emissions stalled or even reversed in most developed countries (Klimont et al., 2013).

The shifting effects of various anthropogenic forcings could have profound implications for humans and natural systems. In particular, in addition to their effects on the correlation between low SM and high VPD, AER emissions can increase the uptake of carbon through rising diffuse radiation, which would be favorable by delivering more light to shaded leaves, improving the efficiency of light use throughout the canopy; such a process is referred as the diffuse-fertilization effect (Mercado et al., 2009; Yue & Unger, 2018). Indeed, it has been suggested that an increased diffuse-radiation effect during the period of global dimming may have enhanced the terrestrial carbon sink by approximately 25% (Mercado et al., 2009). These observations suggest that if aerosol emissions drop while GHGs continue to rise in the near future, we can expect to see decreased diffuse-fertilization effects and higher risks of concurrent drought and aridity, with catastrophic impacts on agricultural production (Asseng et al., 2015) and continental carbon storage (S. Zhou, Zhang, et al., 2019), posing a great challenge to global food security and further accelerating future global warming (Mahowald, 2011).

The presence of anthropogenic forcing decreased the occurrence rate of compound SM and VPD extremes in the historical period, but also increased the concurrent of extreme high temperature and high relative humidity (Figure S10 in Supporting Information S1). High humidity can reduce evaporative cooling of warm-blooded animals, including humans. Extremes of temperature and humidity thus pose a great threat to human health and even life (Li et al., 2020) and may already have been responsible for most of the increase in lethal heat in the historical period (Mora et al., 2017). A renewed strengthening of compound SM drought and atmospheric aridity, however, would be expected to be beneficial in alleviating deadly heat stress by reducing the levels of humidity; if GHG emissions dominate future climate change, no longer offset by the cooling effects of aerosols, this benefit of increased aridity (i.e., reducing relative humidity) may outweigh some of the detrimental consequences of increased air temperatures (Wouters et al., 2022). To summarize, instances of concurrent drought and aridity extremes may reduce lethal human and animal heat stress even as they reduce food production and carbon sequestration. It is clear that better understanding of the coordinating roles of relative humidity can thus be critical for predicting or mitigating the risks of future droughts and lethal heatwaves.

Several uncertainties in this study represent crucial challenges. First, our attribution analysis of LMF has relied on a linear relationship between LMF and coefficients of SM–VPD correlation, that is, $$r(\mathrm{S}\mathrm{M},\mathrm{V}\mathrm{P}\mathrm{D})$$. $$r(\mathrm{S}\mathrm{M},\mathrm{V}\mathrm{P}\mathrm{D})$$ predicts LMF well, both in the Fully-forced and NAT-only simulations, but this relationship can be somewhat spatially different as a result of diverse structures of dependence between SM and VPD (i.e., depends on different copulas across pixels); this may contribute to slight incertitude in the attribution results (S. Zhou, Williams, et al., 2019). Further studies that examine comprehensively the possible spatial heterogeneity and even nonlinearity of the relationship between $$r(\mathrm{S}\mathrm{M},\mathrm{V}\mathrm{P}\mathrm{D})$$ and LMF are thus needed to more precisely attribute regional changes to the strength of negative coupling between SM and VPD. Second, inclusion of anthropogenic influences other than aerosol and GHG emissions, such as irrigation, could also profoundly affect how the models display regional dependence between SM and VPD. For example, Ambika and Mishra (2021) indicated that irrigation increases the flux of latent heat and evapotranspiration, which induces surface cooling and thus decreases instances of VPD. Combined with the supply of SM, irrigation has substantially reduced the frequency of concurrent SM and VPD extremes in India. Our findings and other projections based on CMIP6 or CMIP5 that do not take into account the role of irrigation (S. Zhou, Williams, et al., 2019; S. Zhou, Zhang, et al., 2019), may thus overestimate the risk of compound drought and aridity in heavily irrigated regions. Including these additional anthropogenic factors in the CMIP6 models would therefore be conducive to a better understanding of the long-term changes in compound extremes. Lastly, we have accounted for the impacts of anthropogenic forcing on SM–VPD coupling from a directly climatic perspective, but these direct climatic impacts may be regulated primarily through vegetation. For example, we have suggested that increased aerosol emissions, through their cooling effects, may have played a major role in reducing compound SM and VPD extremes in the historical period, but such effects may have been further amplified by increased diffuse-radiation fertilization on vegetation, leading to increased absorption of atmospheric CO₂ by terrestrial ecosystems (Mercado et al., 2009). In contrast, the strengthened land–atmospheric feedbacks caused by increased GHG-driven warming may also promote the rapid onset of flash droughts (Qing et al., 2022), which can rapidly kill trees and thus further accelerate warming by releasing more CO₂ into the atmosphere (Zhang et al., 2021). These possibilities highlight the importance of understanding the biophysical mechanisms of the interactions among aerosols, GHGs, and meteorology, and of incorporating into ESMs the pathways of the various effects of aerosols and GHGs, in order to improve prediction of large-scale impacts of aerosols and GHGs and their associated climate feedbacks mediated by vegetation.

Conclusions

Our results indicate that the presence of anthropogenic forcing decreased the correlation between SM and VPD and thus the expected occurrence rate of compound drought and aridity in the historical period 1850–2013. Weakened land–atmospheric feedbacks associated with the effect of aerosol cooling appear to have offset the expected positive effect of GHG-induced warming on SM–VPD coupling, leading to the patterns of widespread decreased concurrent drought and aridity in our results. Emissions of aerosols and GHGs may follow different pathways in the future, so a stronger negative coupling between SM and VPD and thus an intensification of compound soil drought and atmospheric aridity is likely. These results call for urgent efforts to adapt to and mitigate the environmental risk associated with compound drought and aridity in the near future.

Acknowledgments

This work was jointly supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant XDA20090000), the National Natural Science Foundation of China (Grant 42271317), and the Innovation Research Team Project of Natural Science Foundation of Hainan Province (Grant 422CXTD515). J.P.’s research was funded by the Spanish Government (Grant PID2019-110521GB-I00), the Catalan Government (Grant SGR2017-1005), the Fundación Ramon Areces (project CIVP20A6621), and the European Research Council (Synergy Grant ERC-SyG-2013-610028, IMBALANCE-P).

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

Data Availability Statement

The monthly surface and total soil moisture, temperature and relative humidity data during 1850–2014 under four simulations from nine Earth System models (ESMs) participating in Phase 6 of the Coupled Model Intercomparison Project (CMIP6) can be assessed online through the Earth System Grid Federation (ESGF) system. The local node used in this study is https://esgf-node.llnl.gov/search/cmip6/. Based the information we provided in Table S1 in Supporting Information S1, CMIP6 data can be easily and publicly downloaded from this node. The Gleam SM data set is available at http://www.gleam.eu/. Enter your email address in the “Downloads” window and then you will receive an SFTP-site (including username and password to the SFTP) to download the Gleam data. The CRU TS temperature, relative humidity, and cloud cover data are available at https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.06/cruts.2205201912.v4.06/.