Introduction

Droughts refer to a period of abnormally dry weather that persists long enough to cause a serious hydrological imbalance1. Droughts can adversely impact various components of natural systems and economic sectors, including but not limited to the deterioration of the ecological environment and wildlife habitat2, reductions of crop yield3,4, water resources5 and hydropower generation6, as well as increasing risk of forest fires7,8. The latter can further undermine the ability of carbon uptake of terrestrial ecosystems and result in more emissions of carbon dioxide into the atmosphere, exacerbating climate change9,10.

Human activities and the resultant global warming have been identified as the main causes of the increasing frequency and intensity of drought events11,12,13. The observed change, however, is the response of drought to all anthropogenic activities. One such activity is the land use and land cover change (LULCC), which manifests primarily as the expansion of croplands and pasturelands at the expense of natural forested and non-forested lands throughout history. Despite the dominant role of greenhouse gases (GHGs) in the global climate, LULCC has been reported to have a greater impact than GHGs do on the regional hydrological cycle14. Biophysically, forests can boost evapotranspiration (ET), and thereby increasing precipitation (P) locally and in downwind regions15,16. Conversely, deforestation usually leads to reductions in P17.

While the important role of LULCC in modifying surface temperature18 and water resources19,20 has been recognized, how LULCC impacts drought characteristics is far from clear. Multiple reasons are responsible for this gap. First, there are no direct observational records of LULCC impact, as our observed signal is the climate response to the total anthropogenic and natural forcing agents, which is further complicated by internal variability from interannual to multidecadal timescales. It is nearly impossible to separate the impact of LULCC from pure observations alone. Second, in the limited observational studies, the impact of LULCC is usually quantified by comparing P observed at two adjacent land types, assuming the same background climate at both locations17,21. By equating the difference from observations made at the same time to changes of land use over time, these studies imply a space-for-time approach. A similar approach is to compare the hydroclimate conditions before and after LULCC at a fixed site22. The main limitation of this kind of study is that the results reveal direct effects of LULCC only, but miss the indirect effects arising from changes in circulation and climate feedback, as these effects were canceled out in the comparison. Besides, the changes in the mean state of P do not necessarily reflect the changes in the extreme climatic events, as changes in variability alone (e.g., longer dry spells) without changes in the mean state can substantially increase the occurrence of extreme hydrological events23,24,25. Third, most published studies only considered the P deficit when analyzing LULCC-induced droughts26,27, neglecting the crucial role of temperature, ET, or potential evapotranspiration (PET), which have been shown to be equally important as P in modifying water balance24,28,29. Lastly, the potential impact of LULCC-induced drought on human society and agriculture (e.g., exposure) is rarely measured. Drought events impact human society mainly by affecting human health30, crop yield4, and natural ecosystem31. An estimation of exposure change could provide more information to evaluate drought impact and help to reduce the exposure of populations and croplands to drought risks32.

A global assessment of LULCC-induced changes on drought is still lacking, as most previous studies are focusing on regional scales, such as the Amazon or a river basin33,34,35. However, the LULCC impact has been reported to be ___location dependent18,36. It is, therefore, not known whether the conclusions drawn at the regional scale are applicable to other regions. In addition, reforestation and forest-protection programs have been pledged worldwide (e.g., the Bonn Challenge) to combat water stress, desertification, and ecosystem degradation37. The impact of deforestation/reforestation on droughts, if any, can provide scientific guidance for policymakers planning reforestation projects.

In this study, we investigate the impact of historical LULCC since 1850 on meteorological drought from the perspective of its frequency, duration, severity, and intensity (see “Methods” section), aiming to bridge these knowledge gaps. To quantify the drought characteristics, the 3-month standardized precipitation-evapotranspiration index (SPEI), which considers roles of both P and PET in determining droughts, is employed in this study38. Our results show that historical LULCC has increased the frequency, duration, severity, and intensity of drought over half of the global land, and the exacerbation of drought events tends to occur in regions with stronger deforestation. The increasing drought frequency and intensity, however, could be substantially alleviated by reforestation and/or avoided deforestation. These findings provide scientific evidence to forest-related policymakers that, in addition to carbon uptake, reforestation can also provide hydrological benefits by alleviating drought events.

Results

Historical LULCC

Two sets of model simulations are employed in this study. The first simulation is the standard historical simulation (hist), which includes all anthropogenic (e.g., greenhouse gases and aerosols) and natural (solar and volcanic activity) forcings, as well as the evolving LULCC from 1850 to 2014. The second simulation is identical to the hist simulation but with LULCC fixed at the 1850 level (hist-noLu). The LULCC impact was quantified as the difference between the two simulations (hist minus hist-noLu; “Methods” section). Both simulations were run in coupled mode. Thus, the LULCC impact quantified here includes both the direct effects arising from changes in the surface energy balance and the indirect effects arising from climate feedback to LULCC. The period of 1950–2014 (65 years) was analyzed in this study. It is noted that historical LULCC does not modify the concentration of atmospheric carbon dioxide in both simulations; thus, only the biophysical effects related to changes in surface properties (e.g., albedo, aerodynamic roughness, and evapotranspiration efficiency) and associated atmospheric feedbacks are identified here.

Historical LULCC since 1850 has manifested mainly as the loss of primary and secondary land (psl, both forested and non-forested) and the concomitant expansion of cropland, pasture, and rangeland. Specifically, psl decreased by ~21.3% of the global land area (28.5 Mkm2) in 2014 relative to that in 1850, including roughly 6.6 Mkm2 forest loss. In the meantime, cropland (+6.4%), pasture (+3.5%), and rangeland (+10.7%) expanded extensively by the indicated fractions (Fig. 1 and Supplementary Fig. 1).

Fig. 1: Historical land use and land cover change.
figure 1

a Fractional change of psl in 1995–2014 relative to 1850. b Same as (a), but for agricultural fields, which are combinations of crops, pastures, and rangelands. c Time series of the global land area fractions for each land use type in 1850–2014.

Full size image

Global responses of drought events to LULCC

In this study, drought events are identified based on the SPEI time series of each grid. The SPEI is a standardized index formulated based on P and PET (see “Methods” section), with negative (positive) values indicating surface dryness (wetness). For a given grid, a drought event was identified when the SPEI value was smaller than –1 for at least 3 consecutive months. The duration is defined as the number of months of each drought event. Severity is defined as the cumulative summation of the SPEI in absolute values in each month within a drought event. The intensity is defined as the mean severity of each month, which is estimated by dividing the severity by the duration. The frequency at each grid is the total number of drought events within the 65-year period. The drought events were identified at each grid in each model for both simulations. By taking their differences, the impact of LULCC was isolated.

Figure 2 shows the changes in the drought characteristics due to historical LULCC (model consistency is shown in Supplementary Fig. 2). About half of the global land (without Antarctica) shows increasing drought frequency, longer duration, enhanced severity and intensity in drought events (Fig. 2). Averaged globally, all four characteristics show positive changes, denoting an increase in drought frequency and the exacerbation of drought events. A substantial fraction of such positive changes is contributed by the changes from the deforested areas (grids with psl changes <−5%). On the other hand, the reforested areas (grids with psl changes ≥5%) consistently show negative changes, indicating an alleviation of drought events. When spatially aggregated, the distributions also shifted significantly (p < 0.001) in all four aspects of droughts. Taking frequency as an example (Supplementary Fig. 3), only 5% of the global land area has a mean frequency greater than 22 within the 65-year period in the hist-noLu experiment, whereas in the hist experiment, approximately 31% of the global land area has a frequency greater than 22. Similar results are also observed for intensity, duration, and severity, but are less pronounced for the latter two. When the drought events are binned by duration (Supplementary Fig. 4a), the frequency of drought in the hist simulation is higher than that in the hist-noLu simulation across all bins and statistically significant for events lasting 6 months or longer. Similarly, the drought frequency increased for all severity and intensity levels, except for those with a severity between 4 and 5 (Supplementary Fig. 4b, c). These analyses demonstrate that historical LULCC has increased the frequency of drought events at all timescales globally and exacerbated drought events at most severity levels.

Fig. 2: Changes in drought characteristics.
figure 2

a Spatial patterns of multi-model mean (MMM) changes in drought frequency (left) and ___domain-averaged bar plot for deforestation grids, grids with no significant LULCC, reforestation grids, and all grids (right), with bars denoting the MMM values and error bars denoting the 95% confidence intervals (CIs). b–d Same as (a), but for duration (b), severity (c), and intensity (d). Gray dots in the spatial maps denote significant changes at the p = 0.05 level. Dots are shown every 5 columns and rows for clarity. The values in the lower-left corner denote the land fractions with (significant) positive change, (significant) negative change, and no change, along with the 95% CIs. Asterisks in the bar plot indicate significant changes at the p = 0.05 level. Red boxes in (a) delineate the 9 selected regions.

Full size image

Regional responses of drought events to LULCC

To better evaluate the drought responses to LULCC, a closer examination is performed on a regional scale. 9 regions with substantial increases in the drought frequency are selected here (red boxes in Fig. 2a; the boundaries are listed in Supplementary Table 3). Supplementary Fig. 5 shows the spatially aggregated whisker plots for all four characteristics of the 9 regions. 5 out of these 9 regions, including South America, North Europe, India, Southeast Asia, and Australia, show significant changes in all four characteristics of droughts. Meanwhile, a few regions, such as East Asia, show significant changes only in frequency and limited changes in other aspects.

Figure 3 displays the joint probability distributions constructed based on duration and intensity for these regions. For each region, a bivariate kernel density is generated based on the MMM values of duration and intensity for all grids within the region. This analysis allows us to observe shifts in the distributions of drought characteristics that are not clearly visible in the mean change. We found that, according to a two-sample Kolmogorov–Smirnov test, all 9 regions exhibited significant (p < 0.001) shifts in the joint distributions. Specifically, all 9 regions show wider distributions of the drought intensity in the hist simulations, with changes found at both tails of the distributions. The shifts in duration, however, are not obvious in Central and Eastern North America, East Europe, and two Asian regions. Jointly, the distributions of all regions in the hist simulation (red) show a larger spread than those of the hist-noLu simulation (blue). Besides, the distributions in all regions except East Europe and East Asia shifted to the upper-right, indicating that the drought events in these 7 regions are getting both longer-lasting and more intense. At the same time, some regions, such as Central and Eastern North America, South America, and East Asia, also experience shorter and less intense droughts. The above analyses demonstrate that historical LULCC not only increased drought frequency in these regions but also altered the concurrent nature of duration and intensity, despite the negligible changes in the mean value.

Fig. 3: Joint probability distributions with and without LULCC.
figure 3

Joint probability distributions of drought duration and intensity from hist-noLu (blue) and hist (red) simulations for Central and Eastern North America (a), South America (b), North Europe (c), East Europe (d), Tropical Africa (e), India (f), East Asia (g), Southeast Asia (h), and Australia (i). The four contours denote the 20%, 40%, 60%, and 80% quantiles, respectively, in each simulation. p < 0.001 indicates that the two joint distributions from hist-noLu and hist are significantly different at the 0.001 level based on a two-sample Kolmogorov–Smirnov test. Univariate distributions of duration and intensity are displayed on the outer axes of each panel.

Full size image

Relationships between the drought response and LULCC

To investigate the possible linkage between LULCC activity and changes in drought characteristics, we binned the changes in frequency as a function of psl fraction changes at intervals of 10% (Fig. 4a). Regions with stronger deforestation activity (more negative psl fraction changes) tend to have larger increases in drought frequencies. On average, the change of drought frequency increases by 0.27 in response to every 10% loss in the psl fraction. Similar results are also observed for other characteristics -- deforestation generally results in a longer duration and stronger severity and intensity of drought events locally. Specifically, the changes in mean duration, severity, and intensity increase by 0.012 months, 0.035, and 0.003, respectively, in response to every 10% loss in the psl fraction (Fig. 4b–d), demonstrating that deforestation tends to locally increase drought risk and exacerbate drought events.

Fig. 4: Changes in drought characteristics as a function of psl fraction changes.
figure 4

Changes in drought frequency (a), duration (b), severity (c), and intensity (d) binned by the fraction changes in primary and secondary land (psl) at intervals of 10%. Negative (positive) psl fraction changes denote deforestation (reforestation). Black dots are the multi-model mean (MMM) values of each bin, while vertical lines denote the 95% CI of each bin. The red lines are least-squares linear regressions based on the MMM values of each bin, with the slope and correlation coefficient values (r) shown at the upper-right corner of each panel. The red shaded regions are the 95% regression estimates based on a bootstrap technique. The inserted gray bar plot in (a) is land fractions of each bin relative to the global land (without Antarctica) on a logarithmic scale.

Full size image

The next question is how LULCC or deforestation activity increases drought frequency? The drought events in the current study are defined based on the standardized water deficit (P − PET) of the surface. Supplementary Fig. 6 shows the spatial changes in P, PET, and P − PET in MMM values. It is seen that all the regions featured with increasing frequency show strong decreases in P − PET (Supplementary Fig. 6c), with a pattern correlation of −0.70 (p < 0.001). In some regions, the negative Δ(P − PET) is mainly due to negative ΔP, such as C.E. North America, whereas in other regions, ΔPET plays a more dominant role, such as Tropical Africa and India. In general, deforestation tends to reduce P and increase PET, thereby decreasing P − PET. Such impact is more pronounced with increasing LULCC magnitude. Thus, the increased drought frequency reported here could be partially explained by the shift in the mean P − PET. It is noted that the mechanisms behind droughts are rather complex and intertwined with many components of the climate system. A lengthened interval between precipitation events without change in the mean value could possibly modify drought frequency and intensity24,25. The detailed mechanisms of drought changes are beyond the scope of this study.

Changes in exposure to droughts

Given the increasing frequency and duration, it is expected that human and terrestrial ecosystems will experience more drought events and, consequently, more disturbances from droughts. In this section, we estimate the exposure of populations, croplands, and forests to drought events (“Methods” section), providing insights and implications for evaluating potential drought impacts on human health, agriculture, and terrestrial ecosystems, as well as the socio-economy.

Drought events can increase the risk of water-borne, vector-borne diseases, such as Diarrhea, due to the low availability of clean water39, and increase the risk of respiratory diseases through reduced air quality due to increased wildfires, erosion, and associated particulate matter30,40, posing serious health impacts on human society41. Using the gridded population data42 and the simulated response of drought frequency to LULCC, we estimated the population exposure to increasing drought frequency driven by historical LULCC in each grid cell. The results were then aggregated to the country/regional scale and are shown in Fig. 5a. Globally, the population exposure increased by ~11422 million person-months (mpm) per decade due to LULCC. Among the 220 countries/regions analyzed here, 105 countries/regions, which are home to 84% of the world’s total population, show increasing exposure. Specifically, population exposure increased by 4400 mpm decade−1 and 2333 mpm decade−1 in India and China, respectively, making them the top two countries on the list. Nigeria, Indonesia, Brazil, and the USA also show substantial increases. The decreasing exposures are mainly seen in West Asia and South Europe regions, where the drought frequency mostly decreases in response to LULCC (Figs. 5a and 2a). When normalizing the total exposure by global population, a person experienced extra 1.6 months (~48 days) of drought per decade due to historical LULCC (Fig. 5b). At the country level, Paraguay shows the largest increase per capita, adding 9.2 months of drought per person every decade, which is followed by Uruguay (5.2 months person−1 decade−1), Gambia (4.2 months person−1 decade−1), and Nicaragua (4.1 months person−1 decade−1).

Fig. 5: Changes in exposure to drought events due to historical LULCC.
figure 5

a Changes in total population exposure in each country per decade (left), along with the Top 10 countries ranked by changes in exposure. The dots are the MMM values, and the error bars are the 95% confidence intervals (CIs). b–d Same as (a) but for population per capita (b), crop (c), and forests (d). The numbers in the parentheses after the country/region names in the bar plots denote the rank of this country/region in the LULCC magnitude, with smaller numbers representing stronger LULCC activity. Note the logarithmic scale on the y-axis in (a), (c), and (d). In (b), the red horizontal dash-dotted line represents the global mean change, while the red shades denote the 95% CI.

Full size image

Multiple lines of evidence suggest that droughts can seriously reduce crop yield43,44, mainly via lowered soil water45, higher temperature46, pest outbreaks and subsequent chemical control47, soil degradation and nutrition loss48, etc. The reduced crop yield may further cause malnutrition, food insecurity (e.g., inflation of food price), and psychological distress49,50,51. Therefore, it is informative to examine the changes in crop exposure to increasing drought frequency. Altogether, the cropland exposure to droughts increased by 26.6 Mkm2 months decade−1, with India (5.5 Mkm2 months decade−1) and Brazil (3.2 Mkm2 months decade−1) leading the increases, followed by Nigeria, Argentina, China, and Australia (Fig. 5c).

Another ecosystem that is sensitive to drought disturbances is the forest. Drought events can severely reduce carbon uptake31,52 and increase the risk of forest fires8,53, both of which can increase the concentration of CO2 in the atmosphere and further amplify global warming. Globally, the forest exposure increases by 28.5 Mkm2 months decade−1 (Fig. 5d), with the largest increases found in Russia (5.4 Mkm2 months decade−1) and Brazil (5.2 Mkm2 months decade−1). Forests in the USA, China, Indonesia, and Congo (Kinshasa) are also increasingly exposed to drought events.

In addition to the total exposure changes, we also analyzed the relative changes in exposure, in which the total exposure changes were scaled by the total exposure in the hist-noLu simulation (Fig. 6). The global mean relative increases in exposure are 13.6%, 14.3%, and 6.5% for population, crop, and forests, respectively. Overall, the countries in South America and tropical Africa are more sensitive to historical LULCC than the countries in other regions. Paraguay has the highest relative change in all three respects. Uruguay, Argentina, Nigeria, Thailand, and Australia also show high relative changes in exposure.

Fig. 6: Relative changes in exposure to drought events due to historical LULCC.
figure 6

a Relative changes in total population exposure in each country per decade (left), along with the Top 10 countries ranked by relative changes in exposure. The dots are the MMM values, and the error bars are the 95% CIs. b, c Same as (a) but for crop (b), and forests (c). In the right column, the numeric values are the global mean changes and the 95% CIs, denoted by the horizontal dash-dotted lines and the shades.

Full size image

The numbers in the parentheses after the country names in Figs. 5 and 6 denote the rank of this country in the LULCC magnitude based on the loss of psl area, with smaller numbers representing stronger LULCC activity (Supplementary Fig. 7). Some countries only had moderate LULCC activity, but experienced much larger increases in exposure, such as Philippines and Vietnam in total population exposure (Fig. 5a) and Papua New Guinea in total forest exposure (Fig. 5d). On the other hand, Kazakhstan, Saudi Arabia, and Mexico has the 5th, 8th, and 9th largest LULCC among all countries, respectively (Supplementary Fig. 7), but they are not on any of the list, implying an inequality between the LULCC implementers and consequence sufferers.

Impact of potential reforestation

The above analyses illustrate that historical LULCC has increased drought risk and exacerbated drought events, especially in the deforested areas. The last question is whether reforestation can alleviate the drought conditions. To answer this question, we estimated the drought characteristics in two sets of model simulations with and without reforestation under the shared socioeconomic pathway (SSP) 3-7.0 scenario54 (“Methods” section). The SSP3-7.0 scenario represents a medium-to-high emission pathway for the 21st century, with a global mean temperature rise of ~2.8 K in 2100 relative to that in 202055. Globally, the land area fraction of psl will decrease from 52.3% in 2014 to 47.0% in 2100 under the standard SSP3-7.0 simulation but increase to 54.5% by 2100 for the reforestation simulation (Supplementary Fig. 8). Differences between these two simulations (SSP3-7.0 with reforestation minus the standard SSP3-7.0) isolate the impact of potential reforestation on droughts under a warming climate.

The impacts of reforestation on drought conditions are shown in Fig. 7 and the model consistencies are shown in Supplementary Fig. 9. Specifically, approximately half of the global land area shows decreasing drought frequency due to reforestation, and statistically significant changes are observed in tropical Africa, Eastern North America and India where large reforestation activities occur (Fig. 7a). On the whole, the global land shows decreasing drought frequency with statistical significance at the p = 0.05 level, which is largely contributed by reforestation grids (Fig. 7a, bar plot). Similar results are also observed for duration, severity, and intensity (Fig. 7b–d), but are not statistically significant. When binned by duration, severity and intensity, the frequency from reforestation simulations also shows decreases over most timescales and levels (Supplementary Fig. 10). Similar to historical results, the drought conditions tend to be mitigated in reforested regions and exacerbated in deforested regions (Fig. 7 and Supplementary Fig. 11). In summary, reforestation and/or avoided deforestation can substantially lower drought risk and alleviate drought events under a warming climate.

Fig. 7: Changes in drought characteristics due to potential reforestation under the SSP3-7.0 scenario.
figure 7

a Spatial patterns of multi-model mean (MMM) changes in drought frequency (left) and ___domain-averaged bar plot for deforestation grids, grids with no significant LULCC, reforestation grids and all grids (right), with bars denoting the MMM values and error bars denoting the 95% confidence intervals (CIs). b–d Same as (a), but for duration (b), severity (c), and intensity (d). Gray dots in the spatial maps denote significant changes at the p = 0.05 level. Dots are shown every 5 columns and rows for clarity. The values in the lower-left corner denote the land fractions with (significant) positive change, (significant) negative change, and no change, along with the 95% CIs. Asterisks in the bar plot indicate significant changes at the p = 0.05 level.

Full size image

Discussion

Using the latest coupled model simulations with and without historical LULCC, our study reveals that historical LULCC has increased the SPEI-based drought frequency and exacerbated drought events in duration, severity, and intensity on both global and regional scales. More longer-lasting and more intense drought events are observed in some regions, such as India and South America, owing to historical LULCC. The deterioration of drought events tends to occur in regions with strong deforestation activity. The increasing drought frequency significantly increased the exposure of populations, croplands, and forests to droughts, posing serious potential impacts on human health, agriculture, and forest ecosystems. Further analyses demonstrate that reforestation has the potential to alleviate drought conditions. This provides a scientific foundation and guidance for forest-related policymakers. The spatial patterns in the changes of drought frequency (Fig. 2) and mean P − PET (Supplementary Fig. 6) reported here are qualitatively consistent with one recent evaluation, in which large-scale idealized deforestation causes reductions in the mean SPEI and thus, a drying surface56. However, the results could not be further compared due to different metrics used. Moreover, different locations and magnitudes of deforestation, as well as different timescales of model integration, all contribute to discrepancies. Compared with the idealized deforestation simulations, our evaluation represents a more realistic LULCC impact and provides more accurate information to local policymakers.

Several implications could be drawn from this study. First, most policymakers only pay attention to carbon uptake and local temperature cooling57,58 when planting trees to mitigate climate change59. Our results demonstrate that reforestation could also provide hydrological benefits by alleviating drought conditions in climate mitigation, which should be considered in future land management policy-making. Second, the inequality between LULCC performers and sufferers (Fig. 5) indicates a non-negligible role of indirect effects of LULCC arising from changes in circulations and the background climate60,61. To put it differently, a country/region that plans to perform reforestation should consider not only the hydrological effects within its borders but also the effects it may cause in remote countries/regions19. Conversely, the consequence of reforestation within its borders is also subject to the impact of LULCC activity in remote regions. Third, the historical LULCC activities are mainly characterized by deforestation and concomitant crop expansion. However, ironically, this process can simultaneously reduce crop yield and degrade the remaining forests through the biophysical effects of droughts (Fig. 5). To compensate for the loss of crop yield, it is imperative to increase crop yield per unit area with more innovative agricultural technology and sustainable land management policies. Otherwise, we need to clear forests for more cropland, which will further exacerbate the climate.

Three limitations exist in the current study. First, the magnitude of LULCC-induced forcing is much smaller than that of GHGs and aerosols62. This leads to limited impacts on a global scale (Fig. 2), limited land fractions with significant changes, and model consistency (Supplementary Figs. 2 and 9). Statistically significant changes are only observed in some regions, such as South America and tropical Africa. It is acknowledged that such a low signal-to-noise ratio adds uncertainties to our results. We stress that future LULCC studies should focus more on regional scales where the LULCC signal is significant (e.g., regions with strong LULCC activity). The second issue is associated with irrigation, a water management practice that can change the surface energy budget (including ET), temperature, near-surface humidity, and precipitation63,64. The joint changes in these variables may further modify drought characteristics via changes in P and PET. Among the 11 models used in this study, only the CESM2 model has active irrigation. Therefore, the remaining models that do not consider the irrigation process may under- or overestimate the biophysical effect of LULCC on droughts. Third, current models are still unable to perfectly represent the biophysical effects of forest cover changes, which may bias the simulated effects on ET65,66 and, subsequently, P and P − PET16. Despite these limitations, our multi-model ensemble study still provides useful information on how LULCC modifies the water cycle biophysically. Future studies with irrigation and more realistic representations of biophysical processes are needed to provide a more complete assessment.

Materials and methods

CMIP6 model output

Eleven models (one realization per model) participating in both the Climate Model Inter-comparison Project Phase 6 (CMIP6)67 and Land Use Model Inter-comparison Project (LUMIP)68 are employed in this study. For each model, two experiments are used to quantify the impact of historical LULCC. The first experiment is a standard historical simulation that includes all anthropogenic forcings (e.g., greenhouse gases, aerosols, and LULCC) and natural forcings (e.g., solar and volcanic forcings). The second is the hist-noLu simulation, in which all forcings are identical to those in the historical experiment but with land use fixed at the 1850 level. The LULCC data in all CMIP6 models are prescribed by the Land-Use Harmonization dataset69. It is worth noting that not all models can fully adopt this LULCC data due to different classification schemes for land cover; thus, the LULCC implemented in each model has different levels of complexity, depending on the model's needs and configurations. All simulations were run in fully coupled mode from 1850 to 2014, and the period of 1950–2014, 65 years in total, is chosen for analyses in this study. The detailed information of each model is listed in Supplementary Table 1. Notably, the atmospheric concentrations of greenhouse gases in both experiments are identical and are not influenced by historical LULCC. Differences between the two experiments are solely attributed to the biophysical effects of historical LULCC.

Model performance evaluation

CMIP6 models have been widely used in the climate community70. However, it is still necessary to evaluate the performance of the model simulations for the main variables used in this study. We compared some selected variables from the hist simulation against observational data. The variables selected for comparison are those used for estimating P − PET (see below), such as temperature, humidity, and radiation. The comparison was performed based on 30-year mean values over 1985–2014 from observational records and the CMIP6 output (Supplementary Fig. 12). In general, the models can reproduce the observed climatology reasonably well for most of the variables examined here. The correlation coefficients (r) between the MMM values and observational values are larger than 0.80 for all the selected variables and larger than 0.95 for temperature and surface radiation. The r between the MMM and the observed P − PET is close to 0.90. Given this evaluation outcome, it is suitable to utilize these models in our analyses.

Calculation of the standardized precipitation-evapotranspiration index (SPEI)

In this study, we employed the SPEI to identify drought events and assess drought characteristics. The SPEI has the advantages of flexible timescales and the inclusion of the impact of background climate change on drought, which has been widely used in drought analyses12,71,72,73,74. The SPEI is a standardized time series of surface water balance, defined as P − PET, in which P is precipitation and PET is the potential evapotranspiration. The latter is the evapotranspiration under the condition without water limitation, which is estimated by the Food and Agricultural Organization Penman–Monteith method75:

$${PET}=\frac{0.408\varDelta \left({R}_{n}-G\right)+\gamma \frac{900}{T+273}u\left({e}_{s}-{e}_{a}\right)}{\Delta +\gamma \left(1+0.34u\right)}$$
(1)

where \({PET}\) is in mm day−1, \(\varDelta\) is the slope of the vapor pressure curve (kPa °C−1), \({R}_{n}\) is the surface net radiation (MJ m−1 day−1), \(G\) is the soil heat flux (MJ m−1 day−1), \(\gamma\) is the psychometric constant (kPa °C−1), \(T\) is the 2 m air temperature (°C), \(u\) is the surface wind speed (m s−1), \({e}_{s}\) is the near-surface saturation vapor pressure (kPa), and \({e}_{a}\) is the actual vapor pressure (kPa).

For each grid cell in each month, the 3-month cumulative P and PET were first obtained by summing all values of that month and the preceding two months. For example, the cumulative P in March 1960 is the summation of P from January 1960 to March 1960. Then, the water deficit on a monthly scale (D) is calculated as follows:

$${D}_{m}={P}_{m}-{{PET}}_{m}$$
(2)

where \({D}_{m}\) is the 3-month cumulative water deficit in month \(m\). \({P}_{m}\) \({{{\rm{and}}}}\) \({{PET}}_{m}\) are 3-month cumulative P and PET, respectively, in month m. The monthly anomalies of D were then obtained by subtracting the monthly climatology from each month at each grid. Previous analysis indicates that the best fitted probability distribution of hydrological variables varies grid by grid76. Considering the different model simulations from different scenarios used in this study, a nonparametric approach76 is used in this study instead of fitting the data into the log-logistic distribution38. This approach generates comparable standardized indices from different simulations with different background climates, and was applied by Chiang et al.12. The monthly D anomalies of each grid in the hist-noLu simulation are first ranked following Gringorten77:

$$p\left(D\right)=\frac{i-0.44}{n+0.12}$$
(3)

where p is the empirical probability of each D value, i is the rank, and n is the sample size, which is 778 in this case (65 × 12–2, there are no D values in the first two months). 0.44 and 0.12 are two empirical coefficients. The empirical probability \(p\) was then translated into the standardized index (SI) using the standard normal distribution function:

$${SI}={\phi }^{-1}\left(p\right)$$
(4)

where \({SI}\) is the SPEI that will be used for analysis. For the hist simulations, the monthly D anomalies were ranked against the hist-noLu simulation to obtain the SPEI values. We used the 3-month SPEI in the current analysis because the 3-month SPEI is capable of capturing seasonal drought events. In addition, we also repeated our analysis with the 6-month SPEI (Supplementary Figs. 13 and 14), and the results are essentially similar to those of the 3-month SPEI.

Drought identification, PDF construction, and frequency change

For a given grid, a drought event was identified when the SPEI value was smaller than –1 for at least 3 consecutive months. The duration is defined as the number of months of each drought event. Severity is defined as the cumulative sum of the SPEI absolute value during a drought event. The intensity is defined as the mean severity of each month within a drought event, which is estimated by dividing the severity by the duration. The frequency at each grid is the total number of drought events within the 65-year period. The drought events were identified at each grid in each model for both simulations. Then the duration, frequency, severity, and intensity were averaged across all events to obtain the mean value of each grid.

To compare the drought frequency for hist and hist-noLu simulations in more detail, we showed the drought frequency grouped by duration, severity and intensity, as shown in Supplementary Fig. 4. In the hist-noLu simulation, we first estimated the total number of drought events by aggregating all the drought events in all land grids in each model, in which the drought events of each grid were weighted by grid area. The frequency for each bin is then obtained by dividing the total number of drought events within that bin by the total number of drought events across all the grids. Notably, the frequency in the hist simulation was estimated against the total number of drought events in hist-noLu for comparison. This process was repeated for all the models, and the MMM values are shown.

The probability distribution functions (PDFs) of each drought feature (Supplementary Fig. 3) were constructed with the spatially aggregated MMM values to examine how historical LULCC has impacted the distribution of each drought characteristic. The PDFs were estimated by fitting a kernel density estimate to the land grids for both the hist and hist-noLu simulations separately. A similar approach was applied to the bivariate distribution constructed from the MMM values of duration and intensity (Fig. 3).

Impact of LULCC on drought

The LULCC impact, denoted as Δ or change, was quantified as differences in the same variables between the two simulations in each model. Taking the drought frequency as an example, the impact of LULCC on frequency is defined as:

$$\Delta {{{\rm{F}}}}{requency}=\overline{{F{requency}}_{h{ist}}}-\overline{{F{requency}}_{h{ist}-{noLu}}}$$
(5)

In Equation [5], Δ is the change in drought frequency due to the impact of LULCC, and the overbars represent the mean value of 1950–2014. The changes in other drought characteristics, such as duration and intensity, as well as other meteorological variables (e.g., P and PET), were estimated similarly by replacing the frequency in Equation [5] with the corresponding variables.

These differences were computed for each grid cell in each model and then averaged to obtain the multi-model mean (MMM) values, with an equal weighting factor assigned to each model. All model outputs were bilinearly resampled to a common spatial resolution of 0.94° × 1.25° in latitude and longitude before data processing.

Significance test

A bootstrap technique was used to test whether the MMM values are significantly different from zero at the grid-cell level. The 11 model values were sampled 11 times randomly with replacement to obtain a mean value. The process was repeated 1000 times to construct a 95% confidence interval (CI) (function bootci in MATLAB). The MMM changes are considered significant if zero falls outside the confidence interval. An advantage of this technique is that it does not require the results of the 11 models to be normally distributed. This approach was also applied to the binned frequency in Supplementary Fig. 4 to examine whether the frequency differences between hist and hist-noLu are statistically significant. To test whether the distributions of drought characteristics from hist-noLu and hist simulations are significantly different (Fig. 3 and Supplementary Fig. 3), a two-sample Kolmogorov–Smirnov test was used to determine the statistical significance.

Exposure to drought

To evaluate the potential impact of drought on populations, croplands, and forests due to LULCC activity, the exposure of populations, croplands, and forests to drought events was estimated. Taking population exposure as an example, the exposure (\({{Exposure}}_{{Pop}}\)) at each grid was estimated as:

$${{Exposure}}_{{Pop}}={\sum }_{1}^{N}{Pop}$$
(6)

In Equation [6], N is the total number of months in drought events within the 65-year period, and \({Pop}\) is the total population of this grid. To eliminate the change in the population over the 65-year period, the population in 2010 was used in the estimation of exposure. The exposure was estimated separately for the hist and hist-noLu simulations. The impact of LULCC was then obtained by taking the differences between the two simulations (hist minus hist-noLu). The results were divided by 6,500,000 to obtain exposure per decade (million person-months), as shown in Fig. 5. The relative changes in exposure were estimated by scaling the total exposure changes with the exposure in the hist-noLu simulations (Fig. 6). The exposures of cropland and forest are estimated similarly by replacing the population data with areas of cropland and forest in 2010 in each grid. The gridded crop and forest data are from the Land-Use Harmonization 2 (LUH2)69. The gridded results were then aggregated to the country/region level. The 95% CI was estimated with the same bootstrap technique as described above.

Impact of reforestation on drought (SSP3-7.0 data analysis)

To evaluate the impact of LULCC under a warming climate, the shared socioeconomic pathway (SSP3-7.0) simulation from the Scenario Model Intercomparison Project54 and the SSP370-126Lu simulation from the LUMIP project are used. The standard SSP3-7.0 scenario represents a medium-to-high emission path with substantial deforestation and limited climate mitigation policies. Under this scenario, the radiative forcing is approximately 7.0 W m−2 at the top of the atmosphere by the end of the 21st century relative to 1850, and the global mean temperature is about 2.8 K warmer than the current level. The simulation covers the period of 2015–2100, with all anthropogenic forcings and LULCC. The SSP370-126Lu simulation is identical to the SSP3-7.0, but LULCC is adapted from the SSP1-2.6 scenario. SSP1-2.6 represents the low-emission developing pathway with strong climate mitigation policies and aims to constrain the global mean temperature change below 2 K by 2100 relative to 1850. In this scenario, forests are substantially recovered at the expense of fraction loss in agricultural land (crop, pasture, and rangeland). Consequently, the SSP370-126Lu simulation represents the reforestation version of SSP3-7.0. The land use data is also prescribed by the Land-Use Harmonization dataset69. Again, the concentrations of greenhouse gases in both simulations are the same, and the LULCC impact identified here is also a pure biophysical effect.

By the time of this study, seven models that archived both SSP3-7.0 and SSP370-126Lu simulations are used for the analysis of future response. The availability of models is listed in Supplementary Table 2. One realization per model is employed. The output of 2036–2100 (65-year period) is used for data analyses to keep consistency with the historical analyses. The SPEI values are calculated similarly to those of historical analyses, and the monthly D anomalies of both SSP3-7.0 and SSP370-126Lu were ranked against the hist-noLu simulation to obtain comparable SPEI values as in Equation [3]. Then, the impact of reforestation (e.g., Fig. 7) is estimated as follows:

$${\Delta {{{\rm{F}}}}{requency}}_{r{eforestation}}=\overline{{{{{\rm{F}}}}{requency}}_{S{SP}370-126{Lu}}}-\overline{{{{{\rm{F}}}}{requency}}_{S{SP}3-7.0}}$$
(7)

The impact of reforestation on other drought characteristics is estimated similarly by replacing frequency in Equation [7] with duration, severity, and intensity.