Understanding meteorological, runoff, and agricultural drought propagation and their influencing factors in an ensemble of multiple datasets

Understanding the propagation of diverse drought conditions is essential for effective drought preparedness. This study evaluated the propagation of meteorological, runoff, and agricultural droughts across global land areas from 1958 to 2024 using an ensemble of reanalysis data (ERA5), land surface model simulations (GLDAS), and merged observational datasets (TerraClimate). Two distinct methodological frameworks were employed to characterize drought propagation: time-lag correlation analysis and multi-threshold run theory. Based on standardized drought indices derived from precipitation, runoff, and soil moisture, the drought propagation characteristics of response time (RT), propagation rate (PR), and lag time (LT) were examined. Moreover, the climatic and geographical factors influencing drought propagation were quantified using the SHapley Additive exPlanations (SHAP)-based attribution method. The results demonstrate the propagation pathways of meteorological-runoff-agricultural drought at the global scale, with the average RT, PR, and LT from meteorological to runoff drought at 5.0 months, 55.3 %, and 1.23 months; from meteorological to agricultural drought at 8.7 months, 30.3 %, and 2.60 months; and from runoff to agricultural drought at 5.8 months, 35.0 %, and 2.49 months, respectively. Drought propagation exhibits pronounced spatiotemporal heterogeneity, closely related to the regional climatic background. Temperature and potential evapotranspiration are the primary factors influencing the propagation of meteorological drought to runoff drought, whereas precipitation plays a decisive role in the propagation from meteorological or runoff drought to agricultural drought. The findings underscore the importance of taking climatic characteristics into account in the development and implementation of regional drought risk management.

1 Introduction

Drought is one of the most frequent natural disasters and is generally defined as a prolonged period of moisture deficits within the water cycle (Liu et al., 2020; AghaKouchak et al., 2023). Under global warming, the magnitude, frequency, and spatial extent of droughts have increased in recent decades, driven by precipitation variability and increased atmospheric evaporative demand (Chen et al., 2025; Gebrechorkos et al., 2025). The intensified droughts pose a significant threat to ecosystems and socio-economic sectors, such as agricultural production (Hendrawan et al., 2022), ecosystem productivity (Cao et al., 2022; Gu et al., 2025), and water resources (Fowler et al., 2022; Liu et al., 2025; Xie et al., 2025). Moreover, multiple climate and hydrological projections reveal a consistent drying trend across many regions during the 21st century (Cook et al., 2020; Christian et al., 2023; Li et al., 2025). The increasing frequency and intensity of droughts are expected to exacerbate water scarcity and ecosystem degradation in the foreseeable future, thus posing significant threats to both the natural environment and human society. Therefore, characterizing the spatio-temporal dynamics of droughts is crucial, as it has significant implications for adapting to and mitigating the impacts of drought-related hazards.

Drought is a complex and multifaceted natural phenomenon (Wu et al., 2022). Although a drought event typically originates from inadequate precipitation and excessive evapotranspiration (referred to as meteorological drought), its impacts on human and natural systems are closely related to subsequent development, such as diminished runoff and streamflow (hydrological drought), reduced soil moisture (agricultural drought), and declined groundwater (groundwater drought). There exists a strong interrelationship among different types of droughts, owing to the close linkage of their driving factors within the hydrological cycle. The transition processes from one type of drought (i.e., meteorological drought) to another (i.e., hydrological drought) are referred to as drought propagation (Apurv et al., 2017; Colombo et al., 2024). Understanding the characteristics of drought propagation is essential for elucidating drought evolution mechanisms and improving early warning systems. These characteristics include propagation time (the time lag between different drought types), propagation probability (the likelihood of transition between drought types), and propagation threshold (the critical level of drought conditions required to trigger a subsequent drought response). Over the past decades, numerous studies have assessed the propagation characteristics of different drought conditions, using identification methods such as time-lag correlation analysis (López-Moreno et al., 2013; Barker et al., 2016; Liu et al., 2023; Geng et al., 2024), run theory (Aryal et al., 2024; Xiong et al., 2025), and event-coincidence analysis (Baez-Villanueva et al., 2024). For example, Shi et al. (2022b) examined the response time from meteorological and hydrological droughts using the maximum correlation coefficient method, and analyzed the variations in response time across different climatic regions. Han et al. (2023) investigated the propagation pathways of various types of droughts across China and revealed the long-chain propagation mechanisms involving meteorological, hydrological, agricultural, and groundwater droughts. Aryal et al. (2024) evaluated the propagation time from meteorological drought to hydrological and agricultural drought across Australia, emphasizing the significant influence of climatic conditions and drought indices in assessing drought propagation dynamics. Among the aforementioned studies, correlation analysis and run theory are two of the most commonly used methodologies for quantifying drought propagation (Zhang et al., 2022). The time-lag correlation analysis based on standardized drought indices (SDIs) with varying accumulation periods provides a direct and effective method for assessing the response time among different drought types from a statistical perspective. However, the correlation analysis can only reflect the average linear relationship between different drought types, but cannot capture the variations among individual drought events (Zhou et al., 2024; Brunner and Chartier-Rescan, 2024). In comparison, the run theory identifies discrete drought events based on the time series of drought indices, thereby providing a more physically meaningful interpretation of the time lag relationship among various drought conditions. Although extensive research has been conducted to examine the characteristics of drought propagation, the results of these studies are heavily dependent on the datasets, evaluation methods, drought indices, and thresholds employed. A comparison of the differences and consistencies in drought propagation characteristics derived from different datasets and methods is desired to improve our understanding of drought propagation, particularly at the global scale.

The propagation characteristics of different types of drought vary depending on climatic conditions and underlying surface conditions (Apurv et al., 2017; Sattar et al., 2019; Apurv and Cai, 2020). Over the past decades, a large number of attribution studies have been conducted to quantify the impacts of climatic and geographical factors on drought propagation, using methods such as statistical analysis (Gevaert et al., 2018), clustering analysis (Liu et al., 2023), causality analysis (Shi et al., 2022a), and machine learning models (Muthuvel and Qin, 2025). Data-driven machine learning approaches are increasingly being employed in drought modeling due to their inherent advantages in capturing nonlinear patterns from complex and high-dimensional data (Sundararajan et al., 2021; Prodhan et al., 2022). Although machine learning models achieve satisfactory simulation accuracy, their reliability remains questionable due to their black-box nature and lack of physical interpretability (Rudin, 2019; Hassija et al., 2024). SHapley Additive exPlanations (SHAP) provides a unified attribution framework for explaining the machine learning outputs, enabling the interpretation of the causal relationships between independent variables and dependent variables (Antwarg et al., 2021; Nohara et al., 2022). Recently, SHAP-based attribution models have been utilized to investigate drought dynamics across various temporal and spatial scales (Xue et al., 2024; Feng et al., 2025). To the best of our knowledge, there is currently a lack of studies that analyze the factors influencing drought propagation from the perspective of SHAP-based machine learning.

Reliable drought monitoring and an in-depth understanding of the underlying mechanisms depend on datasets that accurately describe variations in drought-related hydro-meteorological variables. For large-scale and global assessments, gauge observations and gauge-based gridded datasets are often constrained by limited spatial and temporal coverage, the occurrence of missing values, and challenges in data accessibility (Wang et al., 2020; Gebrechorkos et al., 2024). Numerous satellites, reanalysis, earth system models, and merged datasets have been developed, providing long-term and spatially continuous records of hydro-meteorological variables (Abatzoglou et al., 2018; Hersbach et al., 2020). Over the past decades, extensive efforts have been made to evaluate drought dynamics using different datasets at both regional and global scales. For example, Yuan et al. (2023) assessed the global patterns of flash drought, which is characterized by the rapid depletion of soil moisture, using ERA5 reanalysis and climate model datasets. Gebrechorkos et al. (2025) revealed the critical role of atmospheric evaporative demand in accelerating global drought severity, based on an ensemble of reanalysis, gridded observation, and hydrological model datasets. Wu et al. (2025) conducted an assessment of the dynamic predictability of agricultural drought across global land areas, utilizing the gridded observation, reanalysis, and hydrological model datasets. However, inconsistent findings across studies can be attributed to the inherent uncertainties and errors within different datasets; few systematic comparisons have been conducted to quantify the discrepancies among the multiple datasets in representing drought propagation characteristics (Chen et al., 2020; Huang et al., 2025).

Although previous studies have evaluated drought propagation across various temporal and spatial scales, a comprehensive assessment of the propagation characteristics of meteorological, runoff, and agricultural droughts – derived from ensembles of multiple datasets – remains lacking, particularly at the global scale. Moreover, comparisons among different evaluation methods are also needed to fully understand the drought propagation process and its underlying influencing factors. Therefore, the objectives of this study are as follows: (1) to assess the spatial and temporal patterns of response time, propagation rate, and lag time across meteorological, runoff, and agricultural droughts – derived from an ensemble of multiple datasets; (2) to quantify the main climatic and underlying surface factors that influence drought propagation characteristics from the perspective of machine learning; (3) to compare the robustness and uncertainty associated with different methods and datasets in characterizing drought propagation.

2 Data and methodology

2.1 Datasets

Monthly precipitation, runoff, and soil moisture were derived from the ERA5, the Global Land Data Assimilation System (GLDAS), and TerraClimate datasets to calculate the drought indices. ERA5 is the fifth-generation global atmospheric reanalysis product developed by the European Centre for Medium-Range Weather Forecasts. It integrates extensive records of both in-situ and satellite observations through an ensemble-based data assimilation system (Hersbach et al., 2020). Precipitation in ERA5 was generated by the atmospheric component of the Integrated Forecasting System, whereas runoff and soil moisture were simulated by a land surface model (Boussetta et al., 2021). The soil moisture in ERA5 was aggregated to 1 meter volumetric soil water using weighted data from three layers: 0–7, 7–28, and 28–100 cm. GLDAS is a multi-model ensemble comprising three land surface models – Noah, Catchment, and the Variable Infiltration Capacity – which integrate satellite and in-situ observations through advanced land surface modeling techniques. The soil moisture in GLDAS models has different soil layer structures, all of which were weighted to the root zone depth of 1 m to be consistent with ERA5. TerraClimate integrates multiple datasets, including WorldClim, Climate Research Unit, and Japanese 55-year Reanalysis, to generate hydro-meteorological variables (Abatzoglou et al., 2018). The soil moisture in the TerraClimate refers to the plant extractable soil water based on the root zone storage capacity, as modeled by an empirical water balance model. The original spatial and temporal resolutions of the raw datasets are as follows: ERA5 (1940–present) at 0.25° × 0.25°, GLDAS v2.0 at 1° × 1° (1948–2014), GLDAS v2.0 at 0.25° × 0.25° (2015–2024), and TerraClimate (1958–2025) at $\mathrm{1}/\mathrm{24}$° × $\mathrm{1}/\mathrm{24}$°. To ensure spatial and temporal consistency, the period from 1958 to 2024 was selected as the reference period, and all datasets were uniformly interpolated onto a 1° × 1° latitude–longitude grid using bilinear interpolation.

In addition, the temperature and potential evapotranspiration (PET) were also obtained from the ensemble of ERA5, GLDAS, and TerraClimate datasets. Potential evapotranspiration in these datasets was calculated using the Penman-Monteith method (Abatzoglou et al., 2018). The Normalized Difference Vegetation Index (NDVI) was obtained directly from the Advanced Very High Resolution Radiometer instruments operated by the National Oceanic and Atmospheric Administration (NOAA) (Pinzon and Tucker, 2014). The elevation dataset was obtained from the ETOPO Global Relief Model developed by the National Centers for Environmental Information (https://www.ncei.noaa.gov/products/etopo-global-relief-model, last access: 1 May 2026). The aridity index dataset was derived from the Global Aridity Index and Potential Evapotranspiration Database – Version 3 (Zomer et al., 2022).

2.2 Drought definition and drought indices

Drought is a complex phenomenon related to multiple variables in the water cycle, and there is no universally accepted definition in the current literature (Van Loon, 2015). Drought propagation is a hierarchical top-down process in which meteorological drought, caused by insufficient precipitation, can propagate to other hydrological variables (Teutschbein et al., 2025). A large number of drought indices and datasets have been developed to characterize varying drought conditions at different spatial and temporal scales (AghaKouchak et al., 2023; Gebrechorkos et al., 2025). To provide a consistent and comparable assessment of drought propagation, standardized drought indices (SDIs) derived from precipitation, runoff, and soil moisture were used to define meteorological, hydrological, and agricultural droughts. The SDI time series were obtained by fitting the drought variables of interest to a suitable probability distribution and subsequently normalizing the probabilities to generate a standardized time series. Three SDIs, including Standardized Precipitation Index (SPI) (McKee et al., 1993), Standardized Runoff Index (SRI) (Shukla and Wood, 2008), and Standardized Soil Moisture Index (SSI) (Hao and AghaKouchak, 2013), were computed by fitting parametric probability distributions, specifically the Gamma, log-normal, and normal distributions, respectively. The maximum likelihood estimation (MLE) method was utilized to estimate the parameters of the probability distributions, with the initial values determined based on the L-moments estimation (Stagge et al., 2015). Compared with other drought indices, the SDI is not only simple and efficient to calculate, but also applicable to diverse climates due to its consistent calculation procedure (Zarch et al., 2015; Adnan et al., 2018). Meanwhile, the SDI can be calculated using drought-related variables across multiple time scales, thus enhancing its effectiveness in analyzing drought propagation.

2.3 Response time analysis based on correlation coefficient

The response times among meteorological, hydrological, and agricultural droughts were analyzed using correlation analysis. Assuming that a high correlation coefficient indicates a strong relationship, the correlation analysis between drought indices with different accumulation periods can be conducted to determine the response time of different drought types (Zhang et al., 2022). For example, the response time from meteorological drought to agricultural drought is determined by the accumulation period of SPI that corresponds to the maximum correlation coefficient with the SSI at a 1-month accumulation period. A shorter accumulation period of SPI to 1-month SSI indicates a more rapid agricultural drought response to meteorological drought conditions. The correlation coefficient is calculated using Pearson's correlation coefficient formulated as follows (Pearson, 1896):

$$\begin{array}{}\text{(1)}& r_{\mathrm{P}}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n(x_i-\overline{x})(y_i-\overline{y})}{\sqrt{{\sum }_{i=\mathrm{1}}^n(x_i-\overline{x}{)}^{\mathrm{2}}}\sqrt{{\sum }_{i=\mathrm{1}}^n(y_i-\overline{y}{)}^{\mathrm{2}}}}}\end{array}$$

where r_P represents the Pearson's correlation coefficient between SPI-n (n is the accumulation period, n=1, 2, …, 24) and SSI-1; $\overline{x}$ and $\overline{y}$ represent the average value of SPI and SSI, respectively; x_i and y_i represents the SPI and SSI values in the time i, respectively. The Pearson's correlation coefficient is ranged from −1 to 1. We evaluated the response times from meteorological drought to runoff drought (RT_MR), from meteorological drought to agricultural drought (RT_MA), and from runoff drought to agricultural drought (RT_RA) by analyzing the correlations between SPI and SRI, between SPI and SSI, and between SRI and SSI, respectively.

2.4 Lag time analysis based on run theory

Run theory is a commonly used method for analyzing drought characteristics, which defines the initiation and termination of a drought event based on the drought index. In this study, the drought events were identified using a multi-threshold run theory, which has advantages in avoiding the unreasonable splitting of persistent droughts and filtering out minor drought episodes, thus providing more accurate identification of drought events (Fleig et al., 2006; Ma et al., 2021). Potential drought events were initially identified using an intermediate threshold (X₀=0). Subsequently, the adjacent drought events with an interval of one month and whose drought index values were below a high threshold (X₁=1) within that month were combined. Finally, the potential drought events with one month length and whose drought index value is greater than a low threshold ($X_{\mathrm{2}}=-\mathrm{1}$) were ruled out.

After using run theory to identify the initiation and termination of drought events, the propagation rate and lag time between the two types of droughts can be evaluated. Taking meteorological and agricultural droughts as an example, the propagation from meteorological drought to agricultural drought is defined as the occurrence of an agricultural drought event during the period in which a meteorological drought occurs. Thus, the propagation rate (PR_MA) and lag time (LT_MA) can be mathematically expressed as follows (Sattar et al., 2019):

$$\begin{array}{}\text{(2)}& {\displaystyle }{P}_{\mathrm{MA}}={\displaystyle \frac{n}{m}}\times \mathrm{100}\mathit{\%}\text{(3)}& {\displaystyle }{\text{LT}}_{\mathrm{MA}}={\displaystyle \frac{\sum _{i=\mathrm{1}}^n(T_{\mathrm{M},i}-T_{\mathrm{A},i})}{n}}\end{array}$$

where n is number of meteorological drought events that propagate to agricultural drought events; m is the total number of meteorological drought events during the study period; T_M,i is the starting time of meteorological drought event i, and T_A,i is the starting time of agricultural drought event i. To elucidate the propagation of drought across different types, the SPI, SRI, and SSI at a 1-month accumulation period were used to represent meteorological, runoff, and agricultural drought, respectively. Consistent with the analysis of drought response time, we analyzed the propagation rate and lag time between meteorological and runoff droughts (PR_MR and LT_MR), between meteorological and agricultural droughts (PR_MA and LT_MA), and between runoff and agricultural droughts (PR_RA and LT_RA).

2.5 Assessment of variation across multiple datasets

Two criteria were used to evaluate the variation in drought propagation across different datasets: the coefficient of variation (CV) and the mean absolute deviation (MAD). Specifically, CV and MAD are calculated as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle }\text{CV}={\displaystyle \frac{\sqrt{\frac{\sum _{i=\mathrm{1}}^{\mathrm{3}}(x_i-\overline{x}{)}^{\mathrm{2}}}{\mathrm{2}}}}{\overline{x}}}\text{(5)}& {\displaystyle }\text{MAD}={\displaystyle \frac{\sum _{i=\mathrm{1}}^{\mathrm{3}}\left|x_i-\overline{x}\right|}{\mathrm{3}}}\end{array}$$

where x_i represents the drought propagation characteristics (i.e., RT, PR and LT) for dataset i (i=1, 2, 3).

2.6 Trend analysis of drought propagation

Temporal evolution of drought propagation characteristics was assessed through a moving window approach. For each grid, the drought propagation characteristics (i.e., response time, propagation rate, and lag time) between different types of droughts were calculated using a 30-year moving window that advanced one year at a time, thereby generating an annual time series for the period from 1987 to 2024. The significance of the time series trend was assessed using the Mann-Kendall (MK) test, where a trend was considered statistically significant at the p value <0.05 (Mann, 1945; Kendall, 1975). Given that the series is derived from a moving window with overlapping data, we conducted the Trend-Free Pre-Whitening (TFPW) procedure prior to the MK test to address potential autocorrelation (Yue et al., 2002). The TFPW-MK test effectively removes the influence of serial correlation on trend significance, providing a more robust assessment. In addition, the magnitude of the trend was estimated using Sen's slope estimator (Sen, 1968).

2.7 SHAP-based attribution analysis

The characteristics of drought propagation are closely associated with regional climatic and geographical features. In this study, we quantified the factors that influence drought propagation characteristics across global land areas using the SHapley Additive exPlanations (SHAP) method. SHAP is an effective method for interpreting the outputs of machine learning models, which are often treated as black-box systems, based on principles derived from cooperative game theory (Nohara et al., 2022). The explanation model in the SHAP method can be represented as follows (Antwarg et al., 2021):

$$\begin{array}{}\text{(6)}& G\left(x\right)={\mathit{\phi }}_{\mathrm{0}}+\sum _n^{i-\mathrm{1}}{\mathit{\phi }}_i\end{array}$$

where G(x) is the simulated drought propagation characteristics; φ₀ is the average predicted drought propagation characteristics; and φ_i represents the SHAP value for factor i. For a machine learning model f and input instance x, the SHAP value φ_i for factor i can be formulated as:

$$\begin{array}{}\text{(7)}& {\mathit{\phi }}_i(f,x)=\sum _{S\subseteq N\backslash \mathit{\left\{}i\mathit{\right\}}}{\displaystyle \frac{\left|S\right|\mathrm{!}(\left|N\right|-\left|S\right|-\mathrm{1})\mathrm{!}}{\left|N\right|\mathrm{!}}}\left(f\right(S\cup \mathit{\left\{}i\mathit{\right\}})-f(S\left)\right)\end{array}$$

where N is the set of all factors, and S is a subset of factor excluding factor i. The SHAP value φ_i can quantify the magnitude and direction of factors influencing the model predictions.

According to previous studies, climatic conditions are among the most important factors influencing drought propagation characteristics (Aryal et al., 2024). To explore the relative importance of long-term climatic conditions for drought propagation, the average values (1958–2024) of eight climatic and physiographic variables, including precipitation, temperature, potential evapotranspiration, runoff, soil moisture, aridity index, elevation, and vegetation condition, were selected as model predictors. The corresponding drought propagation characteristics (i.e., response time, propagation rate, and lag time) were selected as target variables. The Extreme Gradient Boosting (XGBoost) model was employed to model the relationships between climatic predictors and drought propagation target variables. The XGBoost model is an efficient and robust gradient-boosted decision tree algorithm that is widely applied in classification and regression tasks within the field of water resources engineering (Chen and Guestrin, 2016; Niazkar et al., 2024). To account for spatial autocorrelation, spatial block cross-validation was employed on the training set to prevent overfitting. The global grid was partitioned into 43 spatially contiguous blocks according to the IPCC AR6 reference land regions (Iturbide et al., 2020). In each fold, ten blocks were held out for validation, and the XGBoost model was trained on the remaining blocks. Model performance was evaluated using the coefficient of determination (R²) and root mean square error (RMSE), averaged across all held-out blocks.

Figure 1 Spatial patterns of average response time from meteorological to runoff droughts (RT_MR), from meteorological to agricultural droughts (RT_MA), and from runoff to agricultural droughts (RT_RA), and the corresponding Pearson correlation coefficients derived from the ensemble of ERA5, GLDAS, and TerraClimate datasets. The blank grids indicate that the correlation between different drought indices is not statistically significant (p value <0.05). The inner plots show the histograms of response time and maximum correlation across global land areas.

3 Results

3.1 Response time and correlation of drought indices

Figure 1 illustrates the spatial patterns of response times among meteorological, runoff and agricultural droughts based on the ensemble of multiple datasets. The results indicate a more rapid propagation of meteorological drought to runoff drought when compared to agricultural drought. The average RT_MR, RT_MA and RT_RA across global land areas are 5.0 months, 8.7 months, and 5.8 months, respectively, with interquartile ranges (IQRs) of [2.7, 6.7] months, [5.0, 11.3] months, and [2.3, 7.3] months. The response time varies across different regions, with shorter propagation times in southeastern Asia, central Africa, South America, and North America. The robustness of the response time evaluation can be attributed to the relatively high correlation coefficients presented in Fig. 1. The IQRs of corresponding correlation coefficients of SPI–SRI, SPI–SSI, and SRI–SSI are [0.43, 0.80], [0.51, 0.68], and [0.52, 0.70], respectively (Fig. 1B, D, and F). The correlation coefficients of response times are relatively high in the mid- to low-latitude regions (30° S–30° N), suggesting strong robustness of the response time measurements in these regions. The response times among meteorological, runoff, and agricultural droughts also exhibit obvious seasonal variations (Figs. S1 and S2 in the Supplement). Shorter response times and higher correlation coefficients were observed during the summer season (June–August in the Northern Hemisphere, and December–February in the Southern Hemisphere).

Figure 2 Spatial patterns of CV and MAD across the ERA5, GLDAS, and TerraClimate datasets for the response time from meteorological to runoff droughts (RT_MR), from meteorological to agricultural droughts (RT_MA), and from runoff to agricultural droughts (RT_RA). Larger values of the CV and MAD signify a more substantial disparity among distinct datasets.

Moreover, the response time also varied across different datasets (Figs. S3 and S4). For example, consistently low RT_MR values were observed in the TerraClimate dataset, whereas high RT_MR values were observed in the ERA5 and GLDAS datasets. On average, the fluctuation ranges of RT_MR, RT_MA and RT_RA among different datasets were [1.96, 7.06] month, [7.87, 10.65] month, and [4.95, 8.00] month, respectively. To quantify the differences among various datasets, two metrics (i.e., CV and MAD) were calculated, and their spatial patterns are shown in Fig. 2. Larger values of CV and MAD indicate more substantial differences among different datasets. The lowest values of CV and MAD were observed in the RT_MA, followed by the RT_RA and RT_MR. This result indicates that the RT_MA showed relatively small variation across different datasets over global land areas, except in the high-latitude regions of the Northern Hemisphere. In contrast, substantial disparities existed in the RT_MR and RT_RA evaluations derived from different datasets, particularly in North America, the Sahara, central Asia, and central Australia. Specifically, consistently low RT_MR values were observed in the TerraClimate dataset, whereas high RT_MR values were predominantly found in the high latitudes of the Northern Hemisphere (Fig. S3). Regarding the maximum correlation of RT_MR, RT_MA, and RT_RA, the GLDAS dataset generally exhibited higher correlation coefficients, whereas the ERA5 and TerraClimate datasets displayed obviously spatial heterogeneity (Fig. S4). This indicates that the response time among different droughts is more reliably represented in the GLDAS dataset.

Figure 3 Spatial patterns of time series trends in RT_MR, RT_MA and RT_RA across global land areas. The blank grids signify that, within at least one time-window in the time series of response time obtained from the moving window, the correlation coefficient is not statistically significant. The black dots indicate the statistical significance of the time series trend, where the p value of the TFPW-MK test is less than 0.05. A significant increase (decrease) indicates that the Sen's slope is greater (less) than 0 and that the p value of the TFPW-MK test is less than 0.05. A nonsignificant increase (decrease) indicates that the Sen's slope is greater (less) than 0 and that the p value of the TFPW-MK test is greater than 0.05. A monotonic trend indicates that Sen's slope is equal to 0.

Figure 3 presents the spatial patterns of the time series trends for RT_MR, RT_MA, and RT_RA, calculated using 30-year moving windows based on the Sen's slope estimator and the M-K test. The temporal trends in response time across different drought types demonstrate notable spatial heterogeneity. The grids exhibiting a monotonic time series trend account for 53.4 %, 35.1 %, and 55.9 % of the total grid for RT_MR, RT_MA, and RT_RA, respectively. In contrast, the percentages of statistically significant increases range from 23.0 % to 26.9 %, whereas the percentages of statistically significant decreases fall within the range of 18.3 % to 31.8 %. For the RT_MR, regions exhibiting increasing trends are primarily located in Europe, northwestern Asia, central Africa, and North America, while decreasing trends are sporadically observed across Asia. The decreasing trends of RT_MA are primarily observed in the mid- to high-latitudes of the Northern Hemisphere, while increasing trends are predominantly evident in Central Asia and North America. In contrast, aside from the obvious increase observed in Central Asia, no clear spatial distribution pattern can be identified in the time series trend of RT_RA. This result indicates that the response time of various drought conditions can vary across different regions, exhibiting significant interannual variability.

Figure 4 Spatial patterns of propagation rate (PR_MR, PR_MA and PR_RA) and lag time (LT_MR, LT_MA and LT_RA) derived from the ensemble of ERA5, GLDAS, and TerraClimate datasets across global land areas. The inner plots show the histograms of propagation rate and lag time across global land areas. The value of LT_MR is lower those that of LT_MA and LT_RA, so it is assigned a different colour bar.

3.2 Propagation rate and lag time of drought events

Figure 4 illustrates the spatial patterns of propagation rate and lag time derived from the ensemble of ERA5, GLDAS, and TerraClimate datasets across global land areas. The propagation rate represents the linkages between two kinds of drought events. A high propagation rate indicates that subsequent drought events are more sensitive to prior drought events. Among three pairs of drought propagation, the propagation rate from meteorological to runoff droughts is highest, with a global average PR_MR value of 55.3 % and an IQR of [46.4, 63.2] %. Spatially, a higher PR_MR (larger than 60 %) tends to occur in mid- and low-latitude regions as well as in low-altitude areas, and is often associated with a shorter LT_MR (less than 1 month). The average LT_MR across global land areas is 1.23 months, with IQR of [0.68, 1.68] months. Consistent with the correlation analysis between SPI and SRI, the propagation from meteorological drought to runoff drought is significantly influenced by regional temperature conditions. Tropical and subtropical regions exhibit a higher propagation rate and a shorter lag time from meteorological drought to runoff drought. The average PR_MA and PR_RA over global land areas are 30.3 % and 35.0 %, respectively, with the IQRs being [19.3, 41.5] % and [23.0, 47.6] %; and average LT_MA and LT_RA are 2.60 and 2.49 months, with the IQRs being [1.71, 2.92] and [1.68, 2.51] months. The higher PR_MA and PR_RA (larger than 40 %) along with the shorter LT_MA and LT_RA (less than 2 months) trend to occur in humid regions, such as eastern North America, South America, central Africa, and southeastern Asia.

Figure 5 Spatial patterns of CV and MAD across the ERA5, GLDAS, and TerraClimate datasets for the propagation rate from meteorological to runoff droughts (PR_MR), from meteorological to agricultural droughts (PR_MA), and from runoff to agricultural droughts (PR_RA).

Figure 6 Spatial patterns of CV and MAD across the ERA5, GLDAS, and TerraClimate datasets for the lag time from meteorological to runoff droughts (LT_MR), from meteorological to agricultural droughts (LT_MA), and from runoff to agricultural droughts (LT_RA).

Figures 5 and 6 illustrate the spatial patterns of CV and MAD for propagation rate and lag time. The spatial patterns of propagation rate and lag time across the ERA5, GLDAS, and TerraClimate datasets are shown in Figs. S3 and S4. For different datasets, the average fluctuation ranges of PR_MR, PR_MA and PR_RA are [44.3, 72.8] %, [26.9, 33.6] %, and [24.9, 43.5] %, and those of LT_MR, LT_MA and LT_RA are [0.69, 1.49] month, [2.05, 3.03] month, and [1.87, 3.46] month, respectively. Consistent with the response time results, the PR_MA and LT_MA exhibited smallest differences across different datasets with low CV and MAD. In comparison, the relatively large differences were observed in the PR_MR, PR_RA, LT_MR, and LT_RA, especially in regions where in-situ observations are scarce. The highest PR_MR and lowest LT_MR values are found in the TerraClimate datasets, with PR_MR values in the low and middle latitudes approaching 90 % and LT_MR values approaching less than 1 month, respectively. In the ERA5 dataset, the relatively low PR_MR and high LT_MR values are found in the high latitudes of the Northern Hemisphere; while in the GLDAS dataset, the relatively low PR_MR and high LT_MR values are found in the hyper-arid regions, such as Sahara and the Arabian Peninsula. The highest PR_MR value is detected in the ERA5 dataset, with an average value of 33.6 % across global land areas, followed by GLDAS (30.2 %) and TerraClimate (26.9 %). The average PR_MA values are 2.05 months, 2.70 months, and 3.03 months for the ERA5, GLDAS, and TerraClimate datasets, respectively. This suggests that the ERA5 reanalysis reveal a significant global sensitivity of agricultural drought to meteorological drought. Similarly, agricultural drought exhibits greater sensitivity to runoff drought, as evidenced by the simulations in the GLDAS dataset. The average PR_RA values are 36.4 %, 43.5 %, and 24.9 % for ERA5, GLDAS, and TerraClimate datasets, respectively, while the corresponding average LT_RA values are 3.46 months, 1.87 months, and 2.08 months.

Figure 7 Spatial patterns of time series trends in the PR_MR, PR_MA and PR_RA across global land areas. The blank grids indicate that the data is missing in one of the ERA5, GLDAS, and TerraClimate datasets. The black dots indicate the statistical significance of the time series trend, where the p value of the TFPW-MK test is less than 0.05. A significant increase (decrease) indicates that the Sen's slope is greater (less) than 0 and that the p value of the TFPW-MK test is less than 0.05. A nonsignificant increase (decrease) indicates that the Sen's slope is greater (less) than 0 and that the p value of the TFPW-MK test is greater than 0.05. A monotonic trend indicates that Sen's slope is equal to 0.

Figure 8 Spatial patterns of time series trends in the LT_MR, LT_MA and LT_RA across global land areas. The blank grids indicate that the data is missing in one of the ERA5, GLDAS, and TerraClimate datasets. The black dots indicate the statistical significance of the time series trend, where the p value of the TFPW-MK test is less than 0.05. A significant increase (decrease) indicates that the Sen's slope is greater (less) than 0 and that the p value of the TFPW-MK test is less than 0.05. A nonsignificant increase (decrease) indicates that the Sen's slope is greater (less) than 0 and that the p value of the TFPW-MK test is greater than 0.05. A monotonic trend indicates that Sen's slope is equal to 0.

Figures 7 and 8 illustrate the spatial distributions of temporal trends in propagation rate and lag time across global land areas, based on a 30-year moving window analysis from year 1958 to 2024. The propagation rate and lag time derived from the runoff theory can be calculated from the grid using continuous meteorological and runoff data; therefore, blank grids indicate missing data in at least one of the ERA5, GLDAS, or TerraClimate datasets. The PR_MR shows a significant decline across 49.0 % of global land areas, while the corresponding LT_MH exhibits a notable decrease in 47.8 %. This indicates a reduced propagation from meteorological drought to runoff drought in nearly half of the global land area. In contrast, the P_MH (LT_MH) shows a significant increase across 28.1 % (29.6 %) of global areas, particularly in the North America and South America. Although more than 42.5 % (50.9 %) of the land area exhibits a significant decreasing trend in PR_MA (LT_MA), notable increasing trends are observed in the western North America, central South America, and central Africa. For the PR_RA, the land areas experiencing increasing trends (significant percentages of 46.7 %) are more extensive than those showing decreasing trends (significant percentages of 31.7 %). Notable increase trends are observed in the western North America, central South America, central Africa, Europe, and northern Asia. Similar spatial patterns are observed in the LT_RA, with significant increases and decreases accounting for 37.0 % and 38.3 %, respectively.

Figure 9 Scatter plots of SHAP values for factors influencing response time, propagation rate, and lag time across global land areas. (P: precipitation, T: temperature, PET: potential evapotranspiration, R: runoff, SM: soil moisture, AI: aridity index, DEM: elevation, NDVI: vegetation condition).

3.3 SHAP-based attribution of drought propagation characteristics

Using different drought propagation characteristics as the target variables, nine XGBoost models were established. The validation sets of these models yielded satisfactory evaluation results (Table S1 in the Supplement), which can substantiate the attribution results. Spatial analysis of drought propagation indicates that regional climatic conditions and physiographical characteristics significantly influence the patterns and dynamics of drought propagation. In this study, we conducted an in-depth analysis of the influence of climate, topography and vegetation conditions on the processes of drought propagation, using a SHAP-based attribution method. Figure 9 presents the SHAP values of eight factors that influence drought propagation characteristics across global land areas, including precipitation, temperature, potential evapotranspiration, runoff, soil moisture, aridity index, elevation, and vegetation condition. The SHAP value reflects both the direction and magnitude of the influence of factors on drought propagation characteristics. A positive SHAP value indicates that the factor has contributed to an increase in response time, propagation rate, and lag time. For example, considering the influence of temperature on the RT_MR, the SHAP value indicates that regions with higher average temperatures trend to have shorter response times from meteorological to runoff droughts (Fig. 9A). The scatter plots of SHAP values indicate that precipitation, temperature, potential evapotranspiration, runoff, aridity index, and NDVI are all positively correlated with the three drought propagation characteristics. Particularly, the elevation is negatively correlated with drought propagation characteristics, whereas soil moisture exhibits an inconsistent correlation with these characteristics.

Figure 10 Ranking of factors influencing drought propagation characteristics based on the absolute SHAP value.

Figure 10 presents the ranking of factors that influence drought propagation characteristics, based on the absolute SHAP values, which is consistent with the scatter plots of SHAP values. The greater the absolute value of SHAP, the more significant its influence on the characteristics of drought propagation. Temperature and potential evapotranspiration are the most significant factors that influence the propagation of meteorological drought into runoff drought. The average absolute SHAP values of temperature (potential evapotranspiration) for RT_MR, PR_MR, and LT_MR are 0.65 (0.37), 3.34 (2.14), and 0.13 (0.24), respectively, ranking first (fourth), first (second), and second (first) among all influencing factors. This result indicates that in warm regions characterized by higher average temperatures and potential evapotranspiration, runoff drought is more sensitive to meteorological drought, thus having a faster response time and a shorter lag time. In comparison, precipitation serves as the main influencing factor in the propagation from both meteorological and runoff droughts to agricultural drought. Among all factors considered for RT_MA, RT_RA, PR_MA, PR_RA, LT_MA, and LT_RA, precipitation exhibits the highest average absolute SHAP value. The differences in SHAP values among other factors are not obvious and vary depending on the distinct characteristics of drought propagation. In humid regions, agricultural drought trends to be more sensitive to meteorological and runoff drought conditions.

Figure 11 Box plots of drought propagation characteristics across global land areas classified by the percentiles of dominant factors. The dominant factor is temperature for RT_MR, PR_MR, PET for LT_MR, and precipitation for the other characteristics.

Figure 11 presents the boxplots of response time, propagation rate, and lag time across global land areas, grouped based on the percentiles of dominant factors identified through SHAP attribution analysis. The meteorological-to-runoff drought propagation characteristics are primarily influenced by regional temperature and PET, with the non-monotonic behaviour predominantly observed in the 30th to 70th percentiles of temperature and PET. In this percentile range, both RT_MR and LT_MR decrease as temperature and PET increase, while PR_MR increases as temperature and PET increase. This result indicates that as temperature and PET increase, runoff drought becomes more sensitive to meteorological drought, and this relationship remains robust within the intermediate ranges of temperature and PET. For example, the average RT_MR is 6.61 months in the 20th to 30th percentiles, decreasing to 3.43 months in the 60th to 70th percentiles; meanwhile, the average RT_MR is 6.88 months in the 0th to 10th percentiles and 4.60 months in the 90th to 100th percentiles. In comparison, the propagations from meteorological and runoff droughts to agricultural drought are mainly derived by the regional precipitation patterns. The average RT_MA, RT_RA, LT_MA, and LT_RA consistently decrease with increasing percentiles of precipitation, whereas the average PR_MA, PR_RA increase with increasing percentiles of precipitation. This indicates that agricultural drought driven by soil moisture deficits is more sensitive to meteorological and runoff droughts in humid regions, showing a strong linear correlation with regional precipitation patterns.

4 Discussion

4.1 Physical interpretation of drought propagation characteristics

In this study, two distinct methodological frameworks were employed to quantify drought propagation: (1) the response time derived from time-lag correlation analysis, and (2) the lag time based on event identification using the run theory. Response time is determined by identifying the accumulation period of a drought index (e.g., SPI) that maximizes its correlation with a target drought index (e.g., SSI at a 1-month accumulation timescale) (López-Moreno et al., 2013; Zhang et al., 2022). This approach reflects the overall synchronicity and statistical memory characteristics of various drought conditions. Thus, the response time values are strongly influenced by long-term variations in regional climatic and hydrological conditions, such as the seasonal cycle, multi-year climate oscillations, and water storage capacity. The response time refers to the system's long-term state that retains a memory of past drought conditions. The evaluation of response time is beneficial for seasonal drought predictability and long-term drought preparedness. The response time also functions as an indicator of the feasibility of using one type of drought index as a proxy for another. For example, due to the lack of comprehensive observational data, the SPI with varying accumulation periods can reflect hydrological, agricultural and groundwater drought conditions (Kumar et al., 2016).

In comparison, lag time is derived from discrete drought events identified using the multi-threshold run theory, which measures the time difference between the onset of one drought event and the onset of another drought event. By focusing on event-based dynamics, the lag time reflects the instantaneous triggering mechanism by which drought signals propagate from the atmosphere to the land surface. Numerous previous studies have analyzed the threshold of extreme stress that triggers drought propagation, using methods such as copula functions, hydrological models, and machine learning (Geng et al., 2024; Yang et al., 2025). The lag time captures the non-linear response mechanism between different drought conditions at a short time scale, which is crucial for real-time early warning and impact assessment.

Our results provide a globally consistent comparison of the response time and lag time for meteorological, runoff, and agricultural drought propagation. The response time of drought propagation (average RT_MR, RT_MA, and RT_RA of 5.0 [2.7, 6.7] months, 8.7 [5.0, 11.3] months, and 5.8 [2.3, 7.3] months) is generally longer than the lag time (average LT_MR, LT_MA, and LT_RA of 1.23 [0.68, 1.68] months, 2.60 [1.71, 2.92] months, and 2.49 [1.68, 2.51] months). This numerical gap arises from differences in the methodology, but both approaches indicate a consistent propagation pathway for meteorological, runoff, and agricultural droughts, with similar spatial patterns. In addition, the machine learning-based attribution method also identifies similar impact factors, which indicates the consistency of drought propagation mechanisms revealed by response time and lag time. This aligns with the conceptual framework of drought propagation, where precipitation deficits (meteorological drought) first influence runoff generation over the land surface (runoff drought), and subsequently affect soil moisture in the root zone (agricultural drought).

4.2 Main factors influencing drought propagation

Across the global land areas, the drought propagation characteristics exhibit notable spatial heterogeneity. Regions such as South America, eastern North America, central Africa, and southeastern Asia demonstrate shorter response times, higher propagation rates, and longer lag times across meteorological, runoff, and agricultural droughts (Figs. 1 and 4). In our analysis, XGBoost models incorporating the SHAP-based attribution method were employed to quantify the impacts of climatic, topographic, and vegetation-related features on drought propagation. Our findings demonstrate that temperature and PET are the key factors influencing the propagation from meteorological to runoff drought, while precipitation predominantly determines the propagation from meteorological/runoff to agricultural drought. This is consistent with previous studies highlighting the significant role of climate characteristics in drought propagation (Apurv et al., 2017). Specifically, the non-monotonic behaviour of meteorological-to-runoff drought propagation characteristics mainly occurred in the range of 20th to the 70th percentiles for temperatures and PET. In the subtropical regions, shorter RT_MR and LT_MR with low PR_MR trend to occur in regions characterized by higher average temperature and PET (Fig. 11). This is primarily attributed to the influence of temperature on the snow-related processes of the water cycle, resulting in a delayed response of runoff to changes in precipitation. During cold seasons, precipitation is stored in the form of snow and ice in glaciers, which subsequently melt and contribute to runoff during the following warm season. Therefore, the sensitivity of runoff drought to meteorological drought is significantly influenced by temperature variations. In comparison, the sensitivity of agricultural drought to meteorological and runoff droughts is closely associated with the regional average precipitation. This is because deep soil layers and aquifers in humid regions generally exhibit high saturation levels, where fluctuations in soil moisture demonstrate a more significant response to variations in precipitation and runoff.

4.3 Comparison of drought propagation in different datasets

A wide variety of meteorohydrological datasets are available for drought monitoring and evaluation. However, a consensus on the most suitable datasets for assessing drought propagation remains elusive across different applications and specific regions. Our results, derived from an ensemble of three different datasets (i.e., ERA5, GLDAS, and TerraClimate), reveal both robust global patterns and notable uncertainties in quantifying drought propagation characteristics. Consistent spatial patterns of drought propagation characteristics – such as shorter RT and LT in tropical and subtropical regions, and longer RT and LT in high-latitude and arid regions – across multiple datasets demonstrate the robustness of drought propagation mechanisms under climatic control. This agreement underscores the fundamental dynamics of drought propagation, which are independent of the methodology and forcing datasets. However, the magnitudes of drought propagation, especially in the meteorological to runoff pathway, also demonstrate significant inter-datasets variability (Figs. 2, 5, and 6). This divergence highlights the inherent uncertainty in drought propagation assessments and points to the distinct strengths and limitations of each dataset.

ERA5 is a high-resolution reanalysis dataset derived from the Integrated Forecasting System, which is forced by atmospheric observations. It generally exhibits higher values of RT_MR, RT_RA, LT_MR, and LT_RA, and lower values of PR_MR and PR_RA in high-latitude regions (Figs. S3–S6). This may more accurately represent the drought propagation in snow-dominated systems where runoff generation processes are complex and exhibit seasonal lags. GLDAS is an ensemble of multiple land surface models and exhibits intermediate drought propagation characteristics with relatively high spatial consistency in correlation coefficients (Fig. S4). This result indicates that the land surface model demonstrates a more robust pattern of drought propagation, and its process-consistent parameterizations may better represent the interrelationships among different drought types. TerraClimate, a statistically downscaled and bias-corrected dataset, consistently yielded the shortest RT_MR and LT_MR and the highest PR_MR, particularly in the mid- to low-latitudes. While this result aligns with the expectation of rapid response in rainfall-dominated regions, the empirical water balance model in TerraClimate may also lead to an overestimation of propagation speed and sensitivity. Our findings rely on the ensemble of multiple datasets, thus avoiding the bias of any single dataset and providing a more robust and consistent understanding of drought propagation.

4.4 Uncertainties and implications in drought propagation evaluation

Drought propagation evaluation relies heavily on drought indices for monitoring and characterizing various drought types. Considering the data availability and the continuity in both temporal and spatial dimensions at the global scale, we employed the SPI, SRI, and SSI to represent meteorological, runoff, and agricultural droughts. Our results demonstrated the propagation pathway of meteorological- runoff-agricultural droughts, which is consistent with previous studies that employed similar indices (Han et al., 2023; Mei et al., 2025). As a multifaceted phenomenon, hydrological drought is a broad term that is related not only to runoff but also to streamflow and the levels of groundwater, lakes, and reservoirs (Van Loon, 2015). Using the drought indices derived from streamflow, the propagation from agricultural to hydrological droughts has also been identified in many studies, particularly at the watershed scale (Odongo et al., 2023; Teutschbein et al., 2025). Runoff is the volume of water that originates from precipitation and flows over the land surface; it is not directly equal to the streamflow in stream channels. A deficit in runoff can directly affect the availability of soil moisture due to reduced recharge to the root zone, representing the propagation from hydrological drought to agricultural drought. In comparison, soil moisture retains precipitation that falls on the land surface and then delays the propagation time form precipitation to streamflow (McColl et al., 2017).

Due to the inherent variability of drought-related variables, significant uncertainties exist within hydrometeorological datasets (Bador et al., 2020). Our findings depend on an ensemble of three datasets (i.e., ERA5, GLDAS, and TerraClimate), which helps avoid biased and incomplete evaluations of drought propagation that could result from relying on a single dataset. We conducted a comparative analysis of drought propagation characteristics derived from multiple datasets, systematically evaluating their consistency and discrepancies (Figs. S3–S6). The results underscore the impact of input data uncertainties on the assessment of drought propagation, with notable discrepancies predominantly observed in the hyper-arid, high-latitude, and high-evaluation regions. Specifically, in hyper-arid regions – where precipitation is extremely low and highly erratic – the calculation of the SPI and SRI becomes statistically unstable; consequently, the correlation between SPI and SRI in such environments can indeed be uninterpretable. This is primarily attributed to the scarcity of in-situ stations capable of providing continuous spatial and temporal observations in these regions. The data assimilation systems and land surface models employed across different datasets to fill missing observations inevitably introduce uncertainties in both model parameters and structural configurations.

Generally, our study provides a comprehensive assessment of multiple drought propagation characteristics across global land areas, which has significant implications for the development and improvement of drought monitoring and early warning systems. In tropical and sub-tropical regions, real-time monitoring of meteorological drought can improve the forecasting of hydrological drought; whereas in humid regions, drought indices based on precipitation and runoff could provide more accurate predictions of agricultural drought risks. Future research could focus on improving the accuracy of predicting future drought changes by incorporating the mechanisms of drought propagation into predictive models. In addition, human activities – such as water abstraction, reservoir regulation, and land-use change – can profoundly modify natural drought propagation processes by altering catchment storage and flow pathways, thereby influencing drought propagation (Wens et al., 2019; Du et al., 2024). While disentangling the effects of human activities on drought propagation at a global scale remains challenging, existing studies have provided valuable evidence of such influences (Long et al., 2026). Future research could also focus on quantitatively disentangling the effects of human activities on drought propagation, particularly at the global scale.

5 Conclusions

In this study, we systematically assessed the propagation characteristics of multiple drought types from 1958 to 2024 across global land areas. Based on an ensemble of multiple datasets (i.e., ERA5, GLDAS, and TerraClimate), three standardized drought indices (SDIs) derived from precipitation, runoff, and soil moisture were employed to represent meteorological, runoff, and agricultural drought conditions, respectively. The lag time derived from correlation analysis, as well as the response time and propagation rate based on run theory, were examined to characterize the propagation of meteorological, runoff, and agricultural droughts. Furthermore, the XGBoost-SHAP model was utilized to quantify the crucial factors that influence drought propagation. Main findings are: (1) The average response time (RT_MR=5.0 months, RT_MA=8.7 months, RT_RA=5.8 months), propagation rate (PR_MR=55.3 %, PR_MA=30.3 %, PR_RA=35.0 %), and lag time (LT_MR=1.23 months, LT_MA=2.60 months, LT_RA=2.49 months) confirm the propagation pathway of meteorological- runoff-agricultural drought at global scale. (1) Over the past six decades, approximately 40 % of the landmass demonstrates temporal variations in the response time of drought propagation, while approximately 70 % of the landmass shows temporal changes in both propagation rate and lag time. (3) Among the eight climatic, topographic, and vegetation-related factors, temperature and potential evapotranspiration are the primary factors influencing the propagation from meteorological drought to runoff drought, primarily due to the lagging effects associated with snow melting processes. (4) In comparison, precipitation predominantly determines the propagation from meteorological or runoff drought to agricultural drought, due to the hydrological process effects of deep soil moisture and aquifer water. In summary, our study presents a multiple data-driven, global perspective on the propagation of meteorological, runoff, and agricultural drought conditions, offering significant implications for drought monitoring and early warning systems in the context of global warming.