Estimating propagation probability from meteorological to ecological droughts using a hybrid machine learning copula method

The impact of droughts on vegetation is essentially manifested as the transition of water shortage from the meteorological to ecological stages. Therefore, understanding the mechanism of drought propagation from meteorological to ecological drought is crucial for ecological conservation. This study proposes a method for calculating the probability of meteorological drought to trigger ecological drought at different magnitudes in northwestern China. In this approach, meteorological and ecological drought events during 1982–2020 are identified using the three-dimensional identification method; the propagated drought events are extracted according to a certain spatiotemporal overlap rule, and propagation probability is calculated by coupling the machine learning model and C-vine copula. The results indicate that (1) 46 drought events are successfully paired with 130 meteorological and 184 ecological drought events during 1982–2020, and ecological drought exhibits a longer duration but smaller affected area and severity than meteorological drought; (2) a quadratic discriminant analysis (QDA) classifier performs the best among the 11 commonly used machine learning models which are combined with four-dimensional C-vine copula to construct the drought propagation probability model; and (3) the hybrid method considers more drought characteristics and a more detailed propagation process which addresses the limited applicability of the traditional method to regions with large spatial extent.

1 Introduction

Drought is a multivariable and complex natural hazard with the characteristics of slow evolution, wide impact, and spatial extent (Feng et al., 2021; Wu et al., 2021; Zhang et al., 2021a, b). Conventionally, drought can be classified into meteorological drought, hydrological drought, agricultural drought, and socioeconomic drought. It is commonly accepted that all types of drought originate from meteorological drought (Mishra and Singh, 2010). Crausbay et al. (2017) argued that existing drought types are described through a “human-centric” lens to characterize a range of effects generated by meteorological drought. This implies that the response of the ecosystem to drought is generally ignored in policy development, which in turn elicits water use conflicts between humans and ecosystems (J. Zhang et al., 2021). The Ecological Drought working group of Science for Nature and People Partnership (SNAPP) proposed a framework of ecological drought from an “ecology-centric” lens, which incorporates ecological, meteorological, and hydrological information (Crausbay et al., 2017). Ecological drought was thus defined as an episodic deficit in water availability that drives ecosystems beyond thresholds of resilience into a vulnerable state, impacts ecosystem services, and triggers feedback in natural and/or human systems (Bradford et al., 2020; Crausbay et al., 2017; McEvoy et al., 2018; Munson et al., 2021; Raheem et al., 2019).

Vegetation is among the most important components in terrestrial ecosystems, and the distribution and growth of vegetation are largely influenced by meteorological factors (Wang et al., 2021; Zeng et al., 2020; Z. Zhang et al., 2021). Developments in remote sensing technology have facilitated the application of vegetation indices to reflect the response of vegetation to climate change (Lawal et al., 2021). For example, a simple linear relationship was found between the standardized precipitation evapotranspiration index (SPEI) and normalized difference vegetation index (NDVI) at a global scale (Vicente-Serrano et al., 2012). The correlation between SPEI and NDVI showed a positive relationship in most regions of northwestern China (NWC), with the exception of a few regions such as the western parts of the Tarim Basin, Qaidam Basin, and southeastern part of the area (Jiang et al., 2018). Actually, the impact of drought on vegetation is manifested as the transition of water shortage from the meteorological stage to the ecological stage. Therefore, this impact should be analyzed by quantifying the effect of decreasing precipitation on the variation of available ecological water, i.e., from the perspective of drought propagation.

Drought propagation refers to the transition of one drought type to another, and it is vital for drought monitoring and prediction (Fang et al., 2020; Warter et al., 2021). Accordingly, drought propagation has become a hot topic in meteorology and hydrology fields (Apurv et al., 2017; Guo et al., 2020). Approaches to drought propagation analysis are broadly divided into model simulations and statistical methods (Han et al., 2019). In the former approach, hydrological responses to meteorological drought are analyzed by using physical-based models that are considered to be effective in representing relevant hydrological processes. Nevertheless, this approach involves labor-intensive calibration processes and is not suitable at large spatial scales (Huang et al., 2017). In contrast, statistical methods with fewer assumptions are easier to use at different spatial scales (Huang et al., 2017). However, in such methods, the propagation process was analyzed using the time series of an average value of drought index in a region or subregion (explained in the Discussion section). In other words, the temporal connection between two drought types is only considered in the traditional statistical methods, but their spatial overlap is ignored, which may result in the miscalculation of drought propagation in regions with large spatial extent.

The probability information of one type of successive drought events is contained in another type of associated drought (Wu et al., 2021). Therefore, a number of studies have attempted to assess the propagation relationships between the two drought types based on the probabilistic method. A Bayesian network is a probabilistic model that acquires probabilistic inferences over interacting variables of interest based on a graphical structure. Therefore, this method has been proven to be suitable for quantifying the probability relationship between different drought types (Ayantobo et al., 2018; Chang et al., 2016; Das et al., 2020). For example, Guo et al. (2020) calculated the occurrence probability of hydrological drought based on different intervals of duration and severities of meteorological drought. Sattar et al. (2019) identified the occurrence probability of different classes and lag times of hydrological drought according to the intensity of meteorological drought. Xu et al. (2021) found that the probability of agricultural drought severity increased synchronously with meteorological drought in different regions of China. Jehanzaib et al. (2020) concluded that in the Korean Peninsula, the probability of meteorological drought propagating into hydrological drought increased significantly under climate change. In general, these studies primarily focused on the relationship between duration and severity between the two drought types but ignored the relationships among affected areas. Xu et al. (2015a) found that the probability of drought occurrence would be underestimated if drought-affected areas are not considered. Therefore, the traditional probabilistic model of drought propagation can be improved by introducing the three-dimensional clustering method, which would provide more drought information (Liu et al., 2019).

Taking a typically ecologically fragile region, northwestern China (NWC), as an example, the motivation of this study is to identify meteorological drought and ecological drought during 1982–2020 in NWC from a three-dimensional perspective, and propose a novel method to investigate the response probability of ecological drought to meteorological drought. The remainder of the current paper is organized as follows: Sect. 2 briefly overviews the geographic information of NWC and describes the datasets used in this paper, and the procedure for estimating propagation probability from meteorological to ecological drought. The results and the comprehensive analysis of the proposed approach are presented in Sects. 3 and 4, respectively. Finally, the conclusions are given in Sect. 5.

Figure 1 Elevation and four divisions of northwestern China.

2 Materials and methods

2.1 Study area

Northwestern China (NWC; 31${}^{\circ }{\mathrm{35}}^{\prime }$–$\mathrm{49}{}^{\circ }{\mathrm{15}}^{\prime }$ N, 73${}^{\circ }{\mathrm{25}}^{\prime }$–$\mathrm{111}{}^{\circ }{\mathrm{15}}^{\prime }$ E) includes the provinces of Shaanxi, Gansu, Qinghai, and the autonomous regions of Xinjiang Uyghur and Ningxia Hui, covering a total surface area of 3.1 million km² (Fig. 1) (Zheng et al., 2021). The terrain of NWC constitutes mountains, basins, and the Gobi. The altitude ranges from −156 to 6647 m, showing the characteristics of “west high and east low”. Four climatic divisions, including humid, semi-humid, semi-arid, and arid areas were demarcated, based on the dryness index (Z. Zhang et al., 2021). As NWC is located at upstream of the Yangtze, Yellow, and other large rivers, it is significant to study the impact of drought on its ecosystem (Liu et al., 2021).

2.2 Datasets

Monthly meteorological data, including surface reflectance, temperature, relative humidity, atmospheric pressure, downward shortwave radiation, wind speed, and longwave radiation, were obtained from the ERA5-Land reanalysis dataset (https://cds.climate.copernicus.eu, last access: 18 November 2022) issued by the European Centre for Medium-Range Weather Forecasts (ECWMF), which has a spatial resolution of $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ and covers the period of 1981–2021. Root soil moisture data were obtained from the hydrological dataset, simulated by the Noah model of the Global Land Data Assimilation System (GLDAS, 0.25${}^{\circ }\times \mathrm{0.25}{}^{\circ }$; https://ldas.gsfc.nasa.gov/gldas, last access: 18 November 2022), covering the period of 1948–2021. NDVI data covering the period 1981–2021 were obtained from the National Centers for Environmental Information (NCEI) (https://www.ncei.noaa.gov/, last access: 18 November 2022), with a spatial resolution of 0.05${}^{\circ }\times \mathrm{0.05}{}^{\circ }$. Land use type data (LUTD) with a spatial resolution 1 km was downloaded from China's multi-period land use/cover change monitoring dataset (http://www.resdc.cn, last access: 18 November 2022); it includes the years of 1980, 1990, 1995, 2000, 2005, 2010, 2015, 2018, and 2020. In order to uniform the spatial resolution of root soil moisture, all spatial datasets were resampled to 0.25${}^{\circ }\times \mathrm{0.25}{}^{\circ }$ using the bilinear interpolation method. The temporal range of all datasets was extracted from January 1982 to December 2020.

2.3 Meteorological and ecological drought index

Previous studies have found that the standardized precipitation evaporation index (SPEI) overestimated the meteorological drought in NWC where actual atmospheric water demand is determined by precipitation variation (Ayantobo and Wei, 2019; B. Zhang et al., 2019; J. Zhang et al., 2019). Additionally, precipitation is the main water resource for vegetation growth in most regions of NWC due to the great depth to groundwater (Cao et al., 2021). A standardized precipitation index (SPI) is thus used in the current study to represent meteorological drought. SPI at different time scales was calculated by aggregating n month moving sums, allowing the identification of various drought types (McKee et al., 1993). At short time scales, drought events are characterized by high frequency and short duration, while at long time scales, they have longer duration and lower frequency. SPI-3 has been reported to be highly representative of the impacts of meteorological conditions on vegetation, as vegetation variation is sensitive to precipitation accumulated over 3 months (McKee et al., 1993; Vicente-Serrano et al., 2012, 2010a). Therefore, SPI-3 is used to characterize meteorological drought in this study. Further details on the SPI calculation are available in McKee et al. (1993).

Figure 2 A schematic diagram illustrating the procedure of the drought propagation identification method.

Commonly used drought indices indirectly reflect the influence of drought on ecosystems, and they do not comprehensively reflect the homeostasis between ecological water consumption and requirement in drought evolution (Jiang et al., 2021). Additionally, decreases in vegetation coverage are not only caused by a persistent deficit in available water for ecosystems but also other aspects, such as wildfire, hail, flood, and human activities (Bento et al., 2020). This limited the ability of vegetation indices to reflect drought conditions. Therefore, a new drought index, the standardized ecological water deficit index (SEWDI), was constructed to monitor terrestrial ecological drought in our previous study (Jiang et al., 2021). SEWDI follows a similar procedure as SPI, which includes the calculation of ecological water deficit (EWD), the selection of an optimal distribution for fitting monthly EWD series, and the inverse normal transformation of the cumulative density distribution of EWD. EWD is the difference between ecological water requirement (EWR) and ecological water consumption (EWC) (Chi et al., 2018; Jiang et al., 2021). Among them, EWR was calculated using the single crop coefficient method recommended by the Food and Agriculture Organization (FAO). EWC equals the actual evapotranspiration, which is derived from latent heat fluxes calculated by the surface energy balance system (SEBS) algorithm. Therefore, SEWDI can reflect the dynamics of energy and water balance under human activities and climate change. Additionally, the standardization method facilitates the same threshold and evaluation criteria in monitoring two drought types (Peng et al., 2019; Zang et al., 2020), which reduces the influence of other algorithms on final results and guarantees spatiotemporal comparability (Liu et al., 2017). The procedure for calculating the SEWDI calculation is detailed in Jiang et al. (2021).

2.4 Drought propagation probability method

Since a reliable understanding of the drought propagation process is beneficial for drought forecasting, research interest in the probability of drought propagation from meteorological droughts to other types of droughts has been increasing (Zhou et al., 2021). The current study thus proposes a novel method coupling spatial and temporal connection methods of two drought types with a machine learning model and C-vine copula to investigate the relationship between meteorological and ecological drought. A flow diagram of the method is depicted in Fig. 2.

Investigating the relationship between the characteristics of the two drought types is key to constructing a probability model. The approach is summarized in three steps as follows:

• Step 1. Meteorological and ecological drought events were identified from a three-dimensional perspective (Sect. 2.4.1).

• Step 2. The two drought types with a genetic relationship were paired on the basis of a certain spatiotemporal matching rule. Their drought characteristics, including drought-affected area, drought severity, and drought duration, were calculated according to the method described in Sect. 2.4.2.

• Step 3. Taking the characteristics of meteorological drought extracted in step 2 as inputs, and propagation results as outputs, the optimal model was selected from 11 machine learning classification models to calculate the propagated probability of meteorological drought (Sect. 2.4.3). Then, a conditional probability model of the paired meteorological and ecological drought events was constructed based on the C-vine copula (Sect. 2.4.4). According to the severities of all identified ecological drought events, cumulative probabilities of 0.5, 0.75, and 0.9 were selected to demarcate moderate, severe, and extreme drought, respectively (Guo et al., 2020). The probabilities of ecological drought at different magnitudes triggered by meteorological drought were obtained by multiplying their conditional probability with the propagated probability of meteorological drought.

2.4.1 Drought identification based on the three-dimensional clustering method

According to Andreadis et al. (2005), the evolution of a drought event should be viewed as a spatiotemporal continuum (longitude, latitude, and time). Different from the traditional one- or two-dimensional drought identification method, the three-dimensional array of SPI-3 and SEWDI-3 were extracted to characterize the degree of meteorological and ecological droughts. The extraction procedure involves two steps (Fig. 2) (Andreadis et al., 2005; Xu et al., 2015a, b): firstly, the clustering method was used to identify drought patches in each month; secondly, the drought continuum was constructed by selecting the overlapping areas of drought patches between two adjacent months, which was greater than 1.6 % of the total area (see Sect. 3.1 for the reason).

For each drought event, three drought characteristics were extracted as follows: (1) affected area was calculated by cumulating the area affected by drought in each month during the entire drought period. (2) Duration denotes the time that a drought event persisted. (3) Severity is a cumulative value of SEWDI-3 or SPI-3 for the entire drought duration and areal extent and equals the volume of the three-dimensional continuum.

Figure 3 Sensitivity test of overlapping areas of drought patches between two adjacent months.

2.4.2 Spatiotemporal connection of two drought types

Liu et al. (2019) developed a new method for identifying the propagation between two related drought types based on the drought identification method from a three-dimensional perspective. The current study employs this method to identify the propagation from meteorological to ecological drought. The key to this method is the determination of the temporal and spatial connection between two drought types. The specific steps are as follows.

Firstly, the identified meteorological and ecological drought events are sorted in chronological order. Secondly, whether the two drought types overlap in time is judged according to Eqs. (1)–(2):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\mathrm{Overlap}}_t=\\ & \left\{\begin{array}{l}\mathrm{1}\text{ if }\left\{\begin{array}{l}{\mathrm{MBT}}_i\le {\mathrm{EBT}}_j\text{ and min }\left({\mathrm{MET}}_i,{\mathrm{EET}}_j\right)\\ -\max \left({\mathrm{MBT}}_i,{\mathrm{EBT}}_j\right)\ge \mathrm{2}\\ {\mathrm{MBT}}_i>{\mathrm{EBT}}_j\text{ and min }\left({\mathrm{MET}}_i,{\mathrm{EET}}_j\right)\\ -\max \left({\mathrm{MBT}}_i,{\mathrm{EBT}}_j\right)\ge \mathit{\alpha }\end{array}\right.\\ \mathrm{0}\text{ if }{\mathrm{MBT}}_i\ge {\mathrm{EBT}}_j\text{ and min }\left({\mathrm{MET}}_i,{\mathrm{EET}}_j\right)\\ -\max \left({\mathrm{MBT}}_i,{\mathrm{EBT}}_j\right)<\mathit{\alpha },\end{array}\right.\end{array}\end{array}$$

$$\begin{array}{}\text{(2)}& \mathit{\alpha }=min\left({\displaystyle \frac{{\mathrm{MDD}}_i}{\mathrm{3}}},{\displaystyle \frac{{\mathrm{EDD}}_j}{\mathrm{3}}}\right),\end{array}$$

where 1 and 0 denote the existence and absence of time overlap between two drought types, respectively; MBT_i and EBT_j represent the beginning time of the ith meteorological and jth ecological drought events, respectively. Similarly, MET_i and EET_j represent the end time of the ith meteorological and jth ecological drought events, respectively; MDD_i and EDD_j indicate the duration of the ith meteorological and jth ecological drought events, respectively.

Thirdly, whether the meteorological and ecological drought patches connecting at a spatial scale is judged according to Eqs. (3) and (4)

$$\begin{array}{}\text{(3)}& {\displaystyle }{\mathrm{Overlap}}_s=\left\{\begin{array}{l}\mathrm{1}\text{ if }{\mathrm{MDA}}_i\cap {\mathrm{EDA}}_i\ge \mathit{\beta }\\ \mathrm{0}\text{ if }{\mathrm{MDA}}_i\cap {\mathrm{EDA}}_i<\mathit{\beta },\end{array}\right.\text{(4)}& {\displaystyle }\mathit{\beta }=\max \left(\mathrm{1.6}\mathit{\%}\cdot A_{\mathrm{NWC}},\min \left({\mathrm{MDA}}_i,{\mathrm{EDA}}_i\right)\cdot b\right),\end{array}$$

where 1 and 0 denote the existence and absence of spatial overlap between two drought types, A_NWC represents the total area of the NWC, and MDA_i and EDA_j represent the projected area of the ith meteorological and jth ecological drought events, respectively. b is set as 15 % in the current study (see Sect. 3.1 for the reason).

Fourthly, successfully matched drought events are encoded following chronological order. Cells in Fig. 2 represent the relationship between preliminarily identified events of the two drought types. The propagation type from meteorological to ecological drought can be classified into four categories: one ecological drought event induced by one meteorological drought event (one-to-one), multiple ecological drought events induced by one meteorological drought event (one-to-many), one ecological drought event induced by multiple meteorological drought events (many-to-one), and multiple ecological drought events induced by multiple meteorological drought events (many-to-many). The codes of cells are identical if the propagation type belongs to one-to-many, many-to-one, and many-to-many.

Finally, the characteristics of meteorological and ecological drought events that belong to the same paired drought event are integrated. Among them, total duration is the difference between the latest-ending and earliest-starting drought events, total affected area is the projected area of all individual drought events, and total severity is the sum of severities of individual drought events.

Figure 4 Identification results of paired meteorological and ecological drought events.

2.4.3 Drought propagation identification based on the machine learning model

The purpose of this part is to identify whether a meteorological drought event has the potential to trigger ecological drought. Eleven commonly used machine learning classification models, including the K-neighbors (KN) classifier (Parzen, 1962), support vector machine (SVM) classifier (Ben-Hur et al., 2000), Gaussian process (GP) classifier (Chen et al., 2020), decision tree (DT) classifier (Quinlan, 1986), multi-layer perceptron (MP) classifier (Cybenko, 1989), AdaBoost (AB) classifier (Freund and Schapire, 1997), Gaussian Naïve Bayes (GNB) (Chan et al., 1982), quadratic discriminant analysis (QDA) (Cover, 1965), gradient boosting (GB) classifier (Friedman, 2001), XGBoost (XGB) classifier (Chen and Guestrin, 2016), and random forest (RF) classifier (Pal, 2005), were employed for propagation judgment. Drought duration, severity, and affected area of meteorological drought were set as the model inputs (Fig. 2). One and zero were set as model target which represents propagation occurrence and non-occurrence, respectively. In this study, each binary classifier was constructed using a Python package called PyCaret, which wraps several machine learning libraries, including scikit-learn, XGBoost, LightGBM, CatBoost, spaCy, Optuna, and Hyperopt (Ali, 2020). The tune_model() function in the PyCaret package offers simple selection of optimal hyperparameters of each model. A five-fold cross-validation was used to train and validate the classifiers in each model by setting “fold=5” in the create_model() function. In using the compare_models() function, the classifier with the highest summation of accuracy, precision, recall, F1 score, and Matthew's correlation coefficient was selected as the optimal model. To avoid overfitting and maintain high calculation efficiency, the L2 regularization method was selected for each model by setting the parameter “penalty=`l2' ”.

$$\begin{array}{}\text{(5)}& {\displaystyle }\mathrm{accuracy}={\displaystyle \frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}}},\text{(6)}& {\displaystyle }\mathrm{precision}={\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}},\text{(7)}& {\displaystyle }\mathrm{recall}={\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}},\text{(8)}& {\displaystyle }{F}_{\mathrm{1}}\mathrm{score}={\displaystyle \frac{\mathrm{2}\cdot \mathrm{TP}}{\mathrm{2}\cdot \mathrm{TP}+\mathrm{FP}+\mathrm{FN}}},\end{array}$$

$$\begin{array}{}\text{(9)}& \mathrm{MCC}={\displaystyle \frac{\mathrm{TP}\cdot \mathrm{TN}-\mathrm{FP}\cdot \mathrm{FN}}{\sqrt{\left(\mathrm{TP}+\mathrm{FP}\right)\cdot \left(\mathrm{TP}+\mathrm{FN}\right)\cdot \left(\mathrm{TN}+\mathrm{FP}\right)\cdot \left(\mathrm{TN}+\mathrm{FN}\right)},}}\end{array}$$

where TP and FN represent actual positives that are correctly and wrongly predicted, respectively; TN and FP represent actual negatives that are correctly and wrongly predicted, respectively. MCC represents Matthew's correlation coefficient, which is used to evaluate the accuracy of binary classification tasks.

2.4.4 Drought propagation probability model based on C-vine copula

Five univariate distributions, including Johnson S_B (Soukissian, 2013), gamma (Thom, 1958), exponential (Marshall and Olkin, 1967), Pearson III (Wallis and Wood, 1985), and Weibull distributions (Thoman et al., 1969), were used to fit affected area, duration, and severity of meteorological drought and severity of ecological drought. The optimal distribution was selected according to the goodness of fit (GOF), which was estimated with the Kolmogorov–Smirnov (KS) test (Marsaglia et al., 2003) and root mean square error (RMSE).

Commonly used copulas, including an elliptical copula (Gaussian) and four Archimedean copulas (Clayton, Gumbel, Frank, and Joe), were used to join two marginal distributions (Chang et al., 2016). The GOF of these copulas was estimated with RMSE and the Cramer-von Mises (CvM) test (Genest et al., 2009).

The vine copula function is an effective tool for integrating different bivariate distributions and calculating the conditional probability of multiple variables (Ni et al., 2020). In a vine copula, an n-dimensional multivariate density is decomposed into $n(n-\mathrm{1})/\mathrm{2}$ bivariate copula densities and arranged into n−1 trees. Among numerous vine copula structures, the C-vine copula has a relatively simple structure and good robustness for constructing multivariate distributions (Wu et al., 2021). Therefore, it is of primary significance to this study. The GOF of C-vine copulas was estimated with RMSE and the CvM test. The joint density function of an n-dimensional C-vine copula is expressed as Eq. (10):

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& f\left(x_{\mathrm{1}},\mathrm{\dots },x_n\right)=\prod _{i=\mathrm{1}}^n{f}_i\left(x_i\right)\times {\prod }_{i=\mathrm{1}}^{n-\mathrm{1}}\prod _{j=\mathrm{1}}^{n-i}{c}_{i,i+j\mathrm{|}\mathrm{1}:\left(i-\mathrm{1}\right)}\\ & \left\{F\left(x_i\mathrm{|}{x}_{\mathrm{1}},\mathrm{\dots },x_{i-\mathrm{1}}\right),F\left(x_{i+j}\mathrm{|}{x}_{\mathrm{1}},\mathrm{\dots },x_{i-\mathrm{1}}\right)\right\},\end{array}\end{array}$$

where $f(x_{\mathrm{1}},\mathrm{\dots },x_n$) represents the joint density function. c represents bivariate copula densities, which include Gumbel, Gaussian, Frank, and Clayton copula functions; F represents the cumulative distribution function of marginal distribution. i and j represent root nodes. More detailed information about the n-dimensional C-vine copula can be referred to Wu et al. (2021). By this means, the conditional probabilities of ecological drought at different magnitudes under impacts of meteorological drought are calculated using Eq. (11):

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& F\left(X>x\mathrm{|}D>d,A>a,S>s\right)\\ & ={\displaystyle \frac{F\left(D>d,S>s,A>a,X>x\right)}{F\left(S>s,A>a,D>d\right)}}\\ & =\left(\mathrm{1}-F\left(d\right)-F\left(s\right)-F\left(a\right)-F\left(x\right)\right.\\ & +C\left(F_{\mathrm{D}}\left(d\right),F_{\mathrm{S}}\left(s\right)\right)+C\left(F_{\mathrm{D}}\left(d\right),F_{\mathrm{A}}\left(a\right)\right)\\ & +C\left(F_{\mathrm{D}}\left(d\right),F_X\left(x\right)\right)+C\left(F_{\mathrm{A}}\left(a\right),F_{\mathrm{S}}\left(s\right)\right)\\ & +C\left(F_{\mathrm{A}}\left(a\right),F_X\left(x\right)\right)+C\left(F_X\left(x\right),F_{\mathrm{S}}\left(s\right)\right)\\ & -C\left(F_{\mathrm{D}}\left(d\right),F_{\mathrm{S}}\left(s\right),F_{\mathrm{A}}\left(a\right)\right)\\ & -C\left(F_{\mathrm{D}}\left(d\right),F_{\mathrm{S}}\left(s\right),F_X\left(x\right)\right)\\ & -C\left(F_{\mathrm{D}}\left(d\right),F_{\mathrm{A}}\left(a\right),F_X\left(x\right)\right)\\ & -C\left(F_{\mathrm{S}}\left(s\right),F_{\mathrm{A}}\left(a\right),F_X\left(x\right)\right)\\ & +C\left(F_{\mathrm{D}}\left(d\right),F_{\mathrm{S}}\left(s\right),F_{\mathrm{A}}\left(a\right),F_X\left(\left.x\right)\right)\right)/\left(\mathrm{1}-F\left(d\right)\right.\\ & -F\left(s\right)-F\left(a\right)+C\left(F_{\mathrm{D}}\left(d\right),F_{\mathrm{S}}\left(s\right)\right)\\ & +C\left(F_{\mathrm{D}}\left(d\right),F_{\mathrm{A}}\left(a\right)\right)+C\left(F_{\mathrm{S}}\left(s\right),F_{\mathrm{A}}\left(a\right)\right)\\ & -C\left(F_{\mathrm{D}}\left(d\right),F_{\mathrm{S}}\left(s\right),F_{\mathrm{A}}\left(\left.a\right)\right)\right),\end{array}\end{array}$$

where D, A, and S represent duration, area, and severity of propagated meteorological drought, respectively; X represents ecological drought at moderate, severe, and extreme magnitudes, which equals the cumulative probability of 0.5, 0.7, and 0.9, respectively. C represents the cumulative distribution function of the joint distribution.

Figure 5 A box plot showing the intensity, duration, and affected area of paired meteorological–ecological drought among different types.

Table 1 Sensitivity test of parameter b.

Table 2 Top 10 meteorological drought events according to severity.

Figure 6 Spatiotemporal continuums of (a) meteorological drought event No. 87 and (b) ecological drought event No. 127.

3 Results

3.1 Threshold selection

Determining overlapping areas of drought patches between two adjacent months is critical in the identification of drought events from a three-dimensional perspective. Sheffield et al. (2009) used 500 000 km² as the area threshold in global scales. For mainland China, 150 000 km² was used as the area threshold in some studies (Wang et al., 2011; Xu et al., 2015b). Liu et al. (2019) took 1.5 % of the total area as the threshold in the Loess Plateau. To determine an optimal area threshold, the number of meteorological and ecological drought events, as well as the ratio of minor drought events, were calculated under different area thresholds. Here, a minor drought is defined as a drought event with 2 months' duration and average SPI/SEWDI larger than −1. As shown in Fig. 3, 1328 and 2305 meteorological and ecological drought events were identified with an area threshold of 0.48 % of the total area of NWC, and the proportions of minor drought events were 44 % and 32 %, respectively. The number of drought events and the proportion of minor drought events decreased with increasing area threshold. However, this trend gradually stabilized when the area threshold was set to be larger than 1.6 % of the total area of NWC, indicating that most minor drought events with the relatively small area were excluded. Therefore, 1.6 % of the total area of NWC was used as the area threshold in this study.

Similarly, the sensitivity of b in Eq. (4) for matching two drought types was tested. The binding mode of absolute and relative thresholds was employed to extract spatial intersection. b is set as 10 %, 15 %, 30 %, 50 %, 70 %, and 90 % to match two drought types. Although some of the successful matching drought events may be merged into one drought event under larger b, the number of successful matching drought events showed little difference under different b (Table 1). In the current study, b=15 % is set, because the most paired drought events could be identified for fitting machine learning models and C-vine copula.

Table 3 Top 10 ecological drought events according to the severity.

3.2 Top 10 meteorological and ecological drought events according to drought severity

A total of 130 meteorological drought events were identified based on SPI-3 from a three-dimensional perspective. The first 10 meteorological drought events in terms of severity in NWC during 1982–2020 are shown in Table 2. Meteorological drought events with longer duration exhibited a relatively larger affected area and were mainly concentrated between 1982 and 2000. Zou et al. (2005) estimated meteorological droughts with the Palmer drought severity index (PDSI) from 1951 to 2003 in China and found that most parts of NWC experienced severe droughts during 1997–2003, which is similar to the results of this study. As shown in Table 2, 2 of 10 meteorological drought events occurred during this period. Moreover, according to the historical record, Xinjiang and Gansu experienced severe meteorological drought during 1985–1986 (J. Zhang et al., 2019). The three-dimensional identification method could sensitively capture these events. Event No. 9 started from southern Gansu in August 1985 and ended in May 1987, and ranked first.

Figure 7 Temporal evolution of severity and area of (a) meteorological drought event No. 87 and (b) ecological drought event No. 127.

Figure 8 Cumulative SPI/SEWDI and migration trajectory of (a) meteorological drought event No. 87 and (b) ecological drought event No. 127.

A total of 184 ecological drought events during 1982–2020 were identified using the three-dimensional identification method. Table 3 lists the top 10 ecological drought events in terms of severity. The most severe ecological drought event started in April 1982 and originated from central Gansu, which was induced by the persistent meteorological drought Nos. 0 and 3. Compared with the characteristics of meteorological droughts, ecological droughts showed a longer duration and a smaller affected area. This reveals that a longer recovery time is required for the mitigation of ecological droughts.

3.3 Identifying propagation from meteorological to ecological drought

A total of 46 paired drought events were successfully matched based on the spatiotemporal connection criterion. As shown in Fig. 4, points representing paired drought events were mainly distributed along a diagonal line, illustrating a relatively high consistency between the two types of droughts on the temporal scale. The number of one-to-one, many-to-one, one-to-many, and many-to-many were 8, 8, 4, and 26, accounting for 17.4 %, 17.4 %, 8.7 %, and 56.5 % of the total number of paired drought events, respectively. Meteorological drought of type one-to-many showed a longer duration, a larger affected area, and a greater severity than ecological drought. However, this is contrary to type many-to-one. Simultaneously, ecological drought of type many-to-one showed a longer duration, a larger affected area, and a greater severity than those of type one-to-many (Fig. 5).

Table 4 Estimations of 11 machine learning models in identifying the potential of meteorological drought to trigger ecological drought. Bold formatting represents the machine learning model with the highest accuracy.

Paired drought event No. 36, comprising meteorological drought event No. 87 and ecological drought event No. 127, was taken as an example to show their spatiotemporal continuums (Fig. 6). The affected area of meteorological and ecological drought in each month was extracted to show their temporal variation. Meteorological drought No. 87 (Fig. 6a) started two months ahead of ecological drought (Fig. 6b), and its effects lasted for two months. It is noteworthy that the most severe meteorological and ecological droughts mainly occurred in central Xinjiang. The affected area and severity of meteorological drought event No. 87 and ecological drought event No. 127 maintained a similar trend of increase–decrease (Fig. 7). Among them, the peaks of the meteorological drought event appeared 2 months ahead (December 2007) of that of the ecological drought (February 2008). In terms of drought trajectory (Fig. 8), they all originated from the Yili basins and showed a counterclockwise shift.

Figure 9 Three-dimensional diagram showing characteristics of meteorological drought events. Larger circles indicate greater severity.

Table 5 Goodness of fit of the marginal distribution.

Table 6 Goodness of fit of the bivariate distribution.

3.4 Propagation probability from meteorological to ecological drought

To estimate the propagated potential of meteorological drought, 11 commonly used machine learning models were trained based on the characteristics of 81 integrated meteorological drought events. As can be seen in Fig. 9, propagated meteorological droughts have greater severity, larger affected area, and longer duration than non-propagated droughts. Table 4 lists the evaluation results of five-fold cross-validations, including accuracy, precision, recall, F1 score, and MCC metrics. The closer these values are to 1, the higher the precision of the model. Therefore, the five metrics were summed to compare the performances of the 11 models. Most models showed a good performance except for Gaussian process and multi-layer perceptron. The QDA classifier with maximum total value was chosen as the best model to identify the propagation potential of meteorological drought.

Table 7 Goodness of fit of the multivariate distribution.

The reliability of the copula function is highly dependent on the dependence between two variables, which was measured by Kendall's τ and Spearman's ρ (Chang et al., 2016; Feng et al., 2021). The τ and ρ between affected area (M_Area), duration (M_Duration), severity (M_Severity) of meteorological drought, and severity of ecological drought (E_Severity) both reached significance at level 0.01 (Tables S1–S2 in the Supplement). The optimal marginal distributions of M_Area, M_Duration, M_Severity, and E_Severity are listed in Table 5. All the distributions passed the KS test and their RMSE was small. Similarly, the parameters of bivariate distribution were estimated using the τ method. The copula estimation can be eased by the itau method, which inverts Kendall's τ method (Demarta and McNeil, 2005). The CvM test and RMSE were used to evaluate their goodness of fit (Table 6). The selected bivariate copulas also demonstrated a good applicability. In the end, the C-vine copula was constructed centered on E_Severity. The CvM test, RMSE (Table 7), and P–P plots (Fig. S1 in the Supplement) indicated that the distribution can be used in probability analysis. The copula structure of M_Area–M_Duration–M_Severity–E_Severity was shown in Table 7.

Figure 10 Conditional probability of ecological drought at (a) extreme, (b) severe, and (c) moderate levels, given that characteristics of meteorological drought exceed a certain value.

Table 8 E_DS with polynomial functions based on meteorological drought characteristics.

Conditional probability is helpful in providing valuable information for the effective allocation of water resources under a certain drought level (Guo et al., 2020). In the current study, the occurrence probabilities of ecological drought at different levels were determined according to the characteristics of meteorological drought (Fig. 10). For example, the occurrence probabilities of moderate, severe, and extreme ecological drought events were 80 %, 63 %, 14.7 %, respectively, when M_DA >17.6×10⁵ km²∩ M_DD >11.8 months ∩ M_DS >7.5×10⁶ months km². Furthermore, the occurrence probability was found to increase more rapidly with increasing M_DS and M_DD compared with M_DA, indicating that the duration and severity of meteorological drought had stronger effects on ecological drought than affected area. Additionally, meteorological drought events with a duration of 2 months but great severity have a high potential to trigger ecological drought. This may be attributable to water shortage induced by meteorological droughts with extremely high intensity (intensity is the drought severity divided by the product of drought duration and affected area).

For comparison, ternary linear and ternary quadratic models were constructed based on 46 pairs of meteorological–ecological drought events (Table 8). The comparisons were made in terms of three independent variables, M_DS, M_DD, and M_DA, and one dependent variable, E_DS. As shown in Table 8, the R² of the ternary quadratic model was evidently higher than that of the ternary linear model, whereas the RMSE, AIC, and BIC were lower. This illustrates that M_DS, M_DD, M_DA, and E_DS follow a nonlinear relationship and that the ternary quadratic model is more suitable for simulating their relationship. According to the ternary quadratic model, E_DS equals 1.4×10⁶ months km⁻² when M_DA >17.6×10⁵ km${}^{-\mathrm{2}}\cap$ M_ DD >11.8 months ∩ M_DS >7.5×10⁶ months km⁻². These values correspond to the thresholds of moderate (1.7×10⁶ months km⁻²), severe (2.4×10⁶ months km⁻²), and extreme (4.6×10⁶ months km⁻²) ecological drought.

Figure 11 Conceptual graph depicting (a) traditional and (b) spatial and temporal connectivity rule of two drought types.

4 Discussion

4.1 Advantages of the proposed approach

Many studies have linked meteorological drought to hydrological drought at different time scales (Ding et al., 2021; Fang et al., 2020; Feng and Su, 2020; Han et al., 2019; Huang et al., 2017; Ma et al., 2019). In these studies, propagated drought events were identified on the basis of the time series between two drought types, and they focused on their lagging, attenuation, lengthening, and pooling (Fig. 11a). The spatial and temporal drought propagation identification method used in the current study not only preserved the characteristics identified by the low dimensional method but also considered the spatial overlap of two drought types (Fig. 11b). Using this method, two types of drought events without spatial connection would be excluded (only 103 out of 184 ecological drought events were induced by 81 out of 108 meteorological drought events), and more drought characteristics, such as affected area, and migration path could be extracted. This addresses the limited applicability of the traditional method to regions with large spatial extent and provides more reliable information for quantifying the relationship between characteristics of meteorological drought and ecological drought. Additionally, we improved the method for calculating the affected area and duration of paired drought events developed by Liu et al. (2019), represented by a simple sum of characteristics of multiple drought events. However, this method overestimates the duration and affected area of some paired drought events, which is inconsistent with the real situation. In this study, the enhanced method could reflect the characteristics of paired drought during the propagation process more accurately.

The conditional probability model was constructed based on paired meteorological and ecological drought events; it is not suitable for calculating the probability of ecological drought at different levels according to meteorological drought events without propagation potential. For example, the probability of moderate ecological drought was 63.3 % if the characteristics of meteorological drought event No. 122 ($\mathrm{M}\mathrm{\_}\mathrm{Area}=\mathrm{5.1}\times {\mathrm{10}}^{\mathrm{5}}$ km⁻², M_Duration=6 months, $\mathrm{M}\mathrm{\_}\mathrm{Severity}=\mathrm{1.89}\times {\mathrm{10}}^{\mathrm{6}}$ months km⁻²) were directly input to the conditional probability model. In reality, this meteorological drought event did not trigger ecological drought. The QDA model added before the C-vine copula was used to address this issue, which could estimate the propagation potential of the corresponding meteorological drought. After this modification, the probability of the propagation of meteorological drought event No. 122 to moderate ecological drought changed to 24.8 %.

4.2 Uncertainty of the model and its improvement measures

The QDA model could simulate the propagation potential of most meteorological drought events well (Table 4). However, some errors occurred in humid southern Shaanxi. For example, meteorological drought event No. 22 showed the potential to trigger ecological droughts, which were incorrectly classified as propagation occurrence. This could be attributed to the compensation of rich water resources for short-term ecological water deficit. Additionally, this paper provides a method for estimating the occurrence probability of ecological drought under the condition of a certain precipitation deficit. The effects of human activities and climate change on ecological drought were not distinguished in the current study. The proposed method may not be accurate for regions with complex water supply systems and strong anthropogenic impacts on vegetation growth.

To improve the accuracy of the method, future studies should consider the non-consistence of ecological drought to quantify the impacts of human activities on drought propagation. Moreover, SPI can be replaced by PDSI or scPDSI to represent meteorological drought through which multiple water balance processes are considered to analyze their relationship with ecological drought (Altunkaynak and Jalilzadnezamabad, 2021). However, such modification may lead to new problems associated with spatiotemporal incomparability. Nevertheless, this approach is worth applying in ecological drought warning. For example, when a meteorological drought event occurs, its characteristics can be applied as input to the tuned model to estimate propagation probability from meteorological to ecological drought in different degrees.

5 Conclusions

This study proposes a method in identifying the propagation probability of meteorological drought events to trigger ecological drought in different magnitudes. Taking NWC as an example, 130 meteorological drought and 185 ecological drought events during 1982–2020 were extracted using the three-dimensional identification method. Compared with meteorological drought, ecological drought events exhibited longer duration but smaller affected area and severity, suggesting that a longer recovery time is required for mitigating ecological drought.

A total of 46 drought events were successfully matched according to a certain spatiotemporal connection principle. The paired drought events were divided into four categories, including one-to-one, many-to-one, one-to-many, and many-to-many. The four categories accounted for 17.4 %, 17.4 %, 8.7 %, and 56.5 % of the total number of paired drought events, respectively. Then, a drought propagation probability model was constructed by coupling QDA and C-vine copula. Compared with the traditional propagation probability model, the proposed model intuitively provides more objective probabilities of ecological drought at different magnitudes.

The current study certainly provides a more robust method for estimating propagation probability from meteorological to ecological drought in similar ecologically fragile regions.