Dynamic Evolution and Copula-Based Multivariable Frequency Analysis of Meteorological Drought Considering the Spatiotemporal Variability in Northwestern China

Abstract

Meteorological drought is a continuous spatiotemporal phenomenon that poses a serious threat to water resource security. Dynamic evolution and multivariable frequency analysis of meteorological drought are important for effective drought mitigation and risk management. Therefore, this study aims to analyze meteorological drought events in northwestern China between 1960 and 2018 based on the standardized precipitation evapotranspiration index (SPEI) through a three-dimensional identification method. This study investigates the meteorological drought dynamic evolution on different time and space scales and evaluates the frequency analysis considering the spatiotemporal variability based on Copula. The results show that SPEI presents an upward trend in Northwestern China. A trend towards increased humidity is observed in arid regions, contrasted by a trend towards aridification in semi-arid and semi-humid areas, indicating that the spatial distribution of drought in the study area tends towards homogenization. The possibility of high-intensity drought events occurring in the same area was relatively low, whereas low-intensity drought events were frequent. Additionally, this study analyzes the dynamic migration process of individual drought events from a three-dimensional perspective. Neglecting any one drought variable could significantly underestimate the occurring probability of severe drought events. Therefore, a multivariable frequency analysis considering the spatiotemporal variability plays a crucial role in the formulation of drought prevention and mitigation strategies, as well as drought forecasting.

1. Introduction

A meteorological drought is a water shortage phenomenon caused by an imbalance between the water supply and demand. It is prolonged by the decrease in precipitation below the climate-normal precipitation. A meteorological drought is one of the most serious natural disasters in the world. In recent years, drought has become a critical environmental problem, threatening human survival, and has aroused widespread concern in the scientific community. The IPCC report (https://www.ipcc.ch/ accessed on 24 October 2023) pointed out that the intensified global warming since the end of the 20th century has led to enhanced surface evaporation and significant regional differences in precipitation, further exacerbating the drought trend. It is therefore urgent to carry out an analysis on the spatiotemporal evolution characteristics and risk of drought to provide a scientific basis for regional drought prevention, drought resistance, and water resource management.

Drought indices can well reflect the wet and dry conditions. The research mechanisms of various drought indices are different, and their applicability in each region is also different. Therefore, an objective and reasonable selection of drought indices is of great practical significance to drought research. Scholars have proposed various drought indices, such as the Palmer Drought Severity Index (PDSI), the Standardized Precipitation Index (SPI), and the Standardized Precipitation Evapotranspiration Index (SPEI). Among them, PDSI is based on the principle of surface water balance to consider the water conditions in the early drought stage. The application of the PDSI is limited by its complex calculations, regional parameters, and single time scale. The SPI has the advantages of simple calculations and multiple time scales, but it ignores the role of temperature in drought formation. In the context of climate change, it is necessary to comprehensively consider the joint effects of precipitation and temperature in drought assessments. In view of this, the SPEI was proposed by Vicente-Serrano et al. [1], which combines the advantages of the PDSI and SPI and performs well in regional drought research.

Understanding the occurrence and spatiotemporal progression of drought is essential for formulating effective strategies to mitigate its impacts. Numerous studies have delved into the spatiotemporal features of drought, aiming to comprehend its characteristics on a regional scale. For example, Tong et al. [2] utilized the Standardized Precipitation Evapotranspiration Index (SPEI) to discern the fluctuations and trends of drought occurrences on the Mongolian Plateau, spanning the period from 1980 to 2014. Sun et al. [3] investigated the spatiotemporal evolution of the 2009–2011 meteorological drought in Southwest China based on the SPEI. Previous research on drought predominantly focused on spatial pattern analysis utilizing data accessible within the specified area [4]. For instance, statistical techniques such as correlation analysis and empirical orthogonal functions were employed to assess the regional attributes of drought occurrences [5,6,7]. Nevertheless, these approaches often simplify the three-dimensional (3D) structure of drought events by reducing the temporal and spatial information into one or two dimensions. For example, this includes examining the temporal changes of fixed-area drought indices or spatial distribution of drought conditions within a specific time frame [8]. In practice, drought events frequently exhibit evident heterogeneity and continuity across various time and space scales. As a result, drought characterization analysis in low-dimensional subspaces cannot well reflect the actual and continuous spatiotemporal evolution characteristics of drought events [9,10].

Copula functions have been widely applied in constructing multivariate joint probability distributions for drought analysis. Over the past few decades, there has been a significant increase in the use of Copula functions for multivariate drought frequency analysis [11,12]. For example, Lee et al. [13] employed Frank and Gumbel Copula functions for bivariate drought frequency analysis. Huang et al. [14] used Bayesian Copula as the best joint distribution function for probabilistic characterization analysis of dry-wet conditions. Ma et al. [15] constructed the trivariate joint distribution of drought duration, severity, and peak based on elliptical Copulas and symmetrical and asymmetrical Archimedean Copulas. Ma et al. [16] proposed the SDF relationships of streamflow drought in the source area of the Yellow River (SAYR) using a Copula-based approach. However, most of these studies focused on temporal features of drought, such as the duration and severity and failed to consider the drought-affected area in the frequency analysis. The spatial characteristics (e.g., drought-affected area) play an important role in characterizing the spatiotemporal structures of drought [17].

Located in the hinterland of the Eurasian continent, Northwestern China is one of the most arid regions at the same latitude and is also an ecologically fragile and climate-sensitive area. Its ecological environment, water resources, and socio-economic development are greatly impacted by regional climate change, and the regional climate responds significantly to global climate warming. Many studies in the 20th century believed that the climate in Northwestern China is generally warm and dry [18,19]. This is confirmed by the rapid retreat of glaciers, the significant degradation of vegetation and oases, the continuous disappearance of wetlands, and a significant reduction in lake areas [20]. In contrast to previous studies, some scholars in the early 21st century discovered that the northwest region has experienced a continuous increase in runoff since the mid-to-late 1980s [21,22]. Based on these findings, researchers proposed that the climate in Northwestern China may be transitioning from warm-dry to warm-wet, with the transition well underway in the western region. Although there has been a certain level of precipitation increase in the western region since the mid-to-late 1980s, the eastern region has experienced a decrease in precipitation.

Therefore, this paper uses the SPEI to study the spatial-temporal patterns and multivariate frequency characteristics of drought in northwest China over the past 60 years. By analyzing the warm-wet and warm-dry patterns in this region, this paper aims to provide a scientific decision-making basis for regional drought prevention, drought resistance, and disaster reduction. The results have important practical significance for understanding the dynamic evolution of drought in northwest China under climate change.

2. Materials and Methods

2.1. Study Area

The study area encompasses Northwestern China, including the three provinces of Shanxi, Gansu, and Qinghai, as well as the Ningxia Hui Autonomous Region and the western portion of the Inner Mongolia Autonomous Region (Figure 1). This region is situated between longitudes of 89°25′~111°27′ E and latitudes of 31°33′~42°48′ N. It covers an approximate land area of 1.87 million km², accounting for roughly 19.46% of the total area in China. The elevation ranges from 184 to 6672 m, exhibiting significant spatial variation with a decrease from west to east. The study area has a diversified topography of mountains, plateaus, and basins, presenting as the Qinghai-Tibet Plateau, Gobi desert, desert steppe, and Loess Plateau from west to east [23].

The study area is mainly affected by the westerlies, Asian monsoons, and plateau monsoons, resulting in dry and cold winters and wet and hot summers [24]. The annual precipitation varies from 50 to 900 mm. The majority of rainfall occurs between June and September, contributing to around 72% of the total annual precipitation. The annual potential evapotranspiration ranges from approximately 800 to 1600 mm, displaying notable spatial diversity and a descending trend from north to south. The mean annual temperature is around 9 °C [25].

2.2. Dataset

In recent years, satellite remote sensing data has been widely applied in Earth science research due to their long time-series, extensive coverage, and high applicability [25]. In this study, monthly precipitation (P) and potential evapotranspiration (PET) data from the CRU TS v.4.03 dataset (https://crudata.uea.ac.uk/cru/data/hrg/ accessed on 18 August 2023) are employed to compute the Standardized Precipitation Evapotranspiration Index (SPEI). The time series covers the years from 1960 to 2018, with a spatial resolution of 0.5° × 0.5°. Additionally, monthly meteorological data of 70 observation stations in the study area obtained from the National Meteorological Information Center are utilized to assess the accuracy of the remote sensing dataset. It is observed that over 85% of the observation stations in the study area exhibited the determination coefficient (R²) values concentrated in the range of 0.7 to 0.9 (with a significance level of p = 0.01). This demonstrates a good consistency between the meteorological data from the observation stations and the CRU dataset [26,27].

2.3. Methods

2.3.1. Standardized Precipitation Evapotranspiration Index

The Standardized Precipitation Evapotranspiration Index (SPEI), developed by Vicente-Serrano [1], is effective for assessing drought conditions under climate change [28,29,30]. The SPEI reflects the drought situation at different time scales through the difference between precipitation and potential evapotranspiration (PET) [31]. The PET is calculated using the Penman-Monteith method, incorporating various climatic variables (e.g., temperature, rainfall, wind speed, and insolation), as shown in formula (1).

$$PET=\frac{0.408\times \Delta \times (R_n-G)+\gamma \times \frac{900}{T+273}\times U_2\times (e_s-e_a)}{\Delta +\gamma \times (1+0.34\times U_2)}$$

(1)

where ∆ is the slope of the saturation vapor pressure temperature relation curve, kPa/°C; R_n is the net amount of radiation reaching the ground, MJ/m²d; G is the soil heat flux density, MJ/m²d; γ is the hygrometer constant; U₂ is the wind speed 2 m above ground level, m/s; e_s is the saturated water vapor pressure of air, kPa; and e_a is the actual vapor pressure, kPa.

Then, the drought severity values of SPEI can be obtained via the following steps:

First, the difference series between precipitation and potential evaporation is given by

$$D_n^k={\displaystyle \sum _{i=0}^{k-1}(P_{n-i}-PE{T}_{n-i})} n\ge k$$

(2)

where P is the precipitation in a particular time scale; PET is the potential evapotranspiration on the same period; i is the month; n is the number of calculations; and k is the time scale for the calculations.

Second, the 3-parameter log-logistic probability distribution function is used to fit the D_n^k series. The probability distribution function F(x) is

$$F(x)={\left[1+{\left(\frac{\alpha }{x-\gamma }\right)}^{\beta }\right]}^{-1}$$

(3)

where α is the scale parameter; β is the shape parameter; and γ is the positional parameter.

Third, the F(x) value is transformed to a standard normal distribution, and then, the SPEI can be calculated using

$$SPEI=W-\frac{C_0+C_{1}W+C_{2}W}{1+d_{1}W+d_2{W}^2+d_3{W}^3}$$

(4)

where C₀ = 2.5155, C₁ = 0.8028, C₂ = 0.0103, d₁ = 1.4327, d₂ = 0.1892, and d₃ = 0.0013.

Table 1. Drought level classification based on SPEI.

Drought Level	SPEI	Drought Severity
I	−0.5 < SPEI	No drought
II	−1.0 < SPEI ≤ −0.5	Mild drought
III	−1.5 < SPEI ≤ −1.0	Moderate drought
IV	−2.0 < SPEI ≤ −1.5	Severe drought
V	SPEI ≤ −2	Extreme drought

presents the drought levels classified according to the SPEI. Negative SPEI values indicate that the climate condition is dry (drought), while positive SPEI values represent a wet climate condition, and SPEI values near zero refer to normal climate conditions. And 1–12 month time scale SPEIs have been calculated for analysis. Droughts are defined when the SPEI value is less than −0.5 in this study.

2.3.2. Modified Mann-Kendall Test (MMK)

The conventional Mann-Kendall test (MK) assumes that the temporal series are random and independent. However, in practice, the time series of meteorological elements often exhibit autocorrelation, which can potentially compromise the reliability of the test results [32]. To address this problem, the Modified Mann-Kendall (MMK) test has been developed, a more rational non-parametric method for assessing the variation in time series. This method effectively addresses autocorrelation within the time series, leading to more reliable trend test results [33]. In this study, the MMK method is applied to analyze the temporal variation trend of drought as well as its spatial distribution characteristics at the grid scale. For a detailed calculation process of the MMK method, please refer to the research of Wang et al. [34].

In terms of the trend statistic value (Z), a positive value indicates an upward trend in the SPEI time series, whereas a negative value indicates a downward trend. Additionally, if the absolute value of Z exceeds 1.64, 1.96, and 2.58, the change trend of SPEI time series is considered significant at the levels of p = 0.1, p = 0.05, and p = 0.01, respectively.

2.3.3. R/S Analysis

The R/S diagnostic method, a non-parametric statistical method, is commonly used to investigate the fractal features and long-term memory dynamics within the time series. Under the context of the R/S analysis, the Hurst index is utilized to predict future trends and reflect the prolonged persistence characteristics and degree of memory in the time series. The detailed steps for the R/S analysis are available in the study of Araujo and Celeste [35].

According to Zhang et al. [36], the Hurst index H ranges from 0 to 1. When 0 < H < 0.5, it signifies that the future trend of the SPEI time series is opposite to its past state, with stronger anti-persistence as H approaches 0. In the case of H = 0.5, it indicates that the SPEI time series is an independent random process with no correlation between the past and future changes. When 0.5 < H < 1, it indicates that the future trend of the SPEI time series maintains consistency with its past state, with stronger persistence as H approaches 1.

2.3.4. Identification of Drought Events by a 3D Clustering Method

Drought identification is the precondition for evolution and frequency analysis of drought. Generally, drought events are characterized by spatial and temporal continuity (longitude, latitude, and time). It is therefore important to consider the spatiotemporal drought structure in the drought analysis [37,38]. In view of this, a 3D clustering method is used to identify the drought events and extract drought variables at different space-time scales. This method consists of two main steps: spatial identification of drought clusters and temporal connection of drought clusters (Figure 2). The detailed calculation process can be found in our previous research [12,25].

Based on the identification results, several spatiotemporal drought variables are defined and calculated. The definitions of these variables are as follows. Drought duration (D) refers to the time interval from the initiation to the termination of one event. Drought-affected area (A) represents the cumulative expanse impacted by a drought event throughout its entire duration. Drought severity (S) indicates the degree of the total water deficit during a drought episode, which is the volume of the three-dimensional drought structure. Drought intensity (I) represents the magnitude of drought in the study area, which can be reflected by the SPEI; the greater the absolute value of SPEI, the more severe the drought.

2.3.5. Copula-Based Multivariable Frequency Analysis of Drought Events

Marginal Distribution of Drought Variables

The accurate fitting of marginal distribution function of drought variables is the basis of the frequency analysis. Inaccuracies in fitting can lead to amplified errors in the calculation of joint distributions. Therefore, it is essential to precisely fit marginal distributions to drought variables. In this study, seven three-parameter distributions are considered [14,39,40]. Seven probability distribution functions, including Gamma, Log-logistic, Lognormal, Weibull, Pearson-III, Generalized Extreme Value, and Generalized Pareto distribution, are used to fit different drought variables.

Table 2. Cumulative probability distribution function and parameters.

Distribution Types	Cumulative Probability Distribution	Parameters
Gam	$F(x)=\frac{{\Gamma }_{x/\beta }(\alpha )}{\Gamma (\alpha )}$	α: shape parameter β: scale parameter
LogL	$F(x)={\left(1+{\left(\frac{\beta }{\alpha }\right)}^{\alpha }\right)}^{-1}$	α: shape parameter (α > 0) β: scale parameter (β > 0)
LogN	$F(x)=\Phi \left(\frac{\ln x-\mu }{\sigma }\right)$	μ: location parameter σ: scale parameter
Wb	$F(x)=1-\exp \left({\left(\frac{x}{\beta }\right)}^{\alpha }\right)$	α: shape parameter β: scale parameter
P-III	$F(x)=\frac{{\displaystyle {\int }_0^{\frac{x-\mu }{\beta }}{t}^{\alpha -1}\exp \left(-t\right)dt}}{\Gamma \left(\alpha \right)}$	α: shape parameter β: scale parameter μ: location parameter
GEV	$F(x)=\{\begin{array}{l}\exp \left(-\left(1+k{\left(\frac{x-\mu }{\sigma }\right)}^{-\frac{1}{k}}\right)\right), k\ne 0\\ \exp \left(-\exp \left(-\left(\frac{x-\mu }{\sigma }\right)\right)\right), k=0\end{array}$	k: shape parameter σ: scale parameter (σ>0) μ: location parameter
GP	$F(x)=\{\begin{array}{l}1-{\left(1+k\frac{\left(x-\mu \right)}{\sigma }\right)}^{-\frac{1}{k}}, k\ne 0\\ 1-\exp \left(-\frac{\left(x-\mu \right)}{\sigma }\right), k=0\end{array}$	k: shape parameter σ: scale parameter (σ > 0) μ: location parameter

shows the probability distribution function and related parameters.

In addition, the Kolmogorov-Smirnov (K-S) test is used to determine whether the assumed distribution significantly differs from the underlying distribution of the actual data at a significance level of p = 0.05. When p > 0.05 and the K-S statistic is below the critical value, it indicates that the assumed distribution agrees well with the underlying distribution of the drought variable. A smaller K-S statistic corresponds to a better fit of the distribution.

Determination of the Optimal Copula Function

The determination of the Copula function also greatly impacts the frequency analysis results of multivariate drought variables. Initially, the K-S method is used to assess whether alternative Copula functions pass the goodness-of-fit test. Subsequently, the relevance between alternative Copula functions and empirical Copula functions is measured using Root-Mean-Square Error (RMSE), Akaike Information Criterion (AIC), and Bayesian Information Criterion (BIC) [41]. The optimal Copula function is then determined based on the goodness-of-fit test results, which is used to establish the joint distribution among drought variables.

For simplicity and representativeness [42], six Copula functions are employed as candidate Copula functions. These include theoretical Archimedean Copulas (Frank, Clayton, Gumbel, and Joe) and elliptical Copulas (Normal and t) [43]. The estimation of Copula parameters was carried out through the inversion of Kendall’s (tau) method. Detailed information about the Copula functions is given in

Table 3. Function expressions of six Copula functions.

Copula	Function Expression	Parameters
Frank	$C_F\left(u_1,u_2,\cdots ,u_d;\theta \right)=-\frac{1}{\theta }\ln \left[1+\frac{{\displaystyle \prod _{j=1}^d{e}^{-\theta u_j-1}}}{{\left(e^{-\theta }-1\right)}^{d-1}}\right]$	$\theta \in R\backslash 0$
Clayton	$C_C\left(u_1,u_2,\cdots ,u_d;\theta \right)=\left[\left({\displaystyle \sum _{j=1}^d{u}_j^{-\theta }}\right)-d+1\right]$	$\theta \in \left[-1, \infty \right)\backslash 0$
Gumbel	$C_G\left(u_1,u_2,\cdots ,u_d;\theta \right)=\exp \left\{-{\left[{\displaystyle \sum _{j=1}^d{\left(-\ln u_j\right)}^{\theta }}\right]}^{\frac{1}{\theta }}\right\}$	$\theta \in \left[1, \infty \right)$
Joe	$C_J\left(u_1,u_2,\cdots ,u_d;\theta \right)=1-{\left[{\displaystyle \sum _{j=1}^d{\left(1-u_j\right)}^{\theta }}-{\displaystyle \prod _{j=1}^d{\left(1-u_j\right)}^{\theta }}\right]}^{\frac{1}{\theta }}$	$\theta \in \left[-1, \infty \right)$
Normal	$C_N\left(u_1,u_2,\cdots ,u_d;\sum \right)=\Phi \left({\Phi }^{-1}\left(u_1\right),\cdots ,{\Phi }^{-1}\left(u_d\right)\right)$	$\sum =\left[\begin{array}{cc}1& \cdots \\ ⋮& \ddots \\ {\rho }_{1d}& \cdots \end{array}\begin{array}{c}{\rho }_{1d}\\ ⋮\\ 1\end{array}\right]$
Student t	$C_t\left(u_1,u_2,\cdots ,u_d;\sum ,\upsilon \right)=T_{\sum ,\upsilon }\left(T_{\upsilon }{}^{-1}\left(u_1\right),\cdots ,T_{\upsilon }{}^{-1}\left(u_d\right)\right)$	$\sum =\left[\begin{array}{cc}1& \cdots \\ ⋮& \ddots \\ {\rho }_{1d}& \cdots \end{array}\begin{array}{c}{\rho }_{1d}\\ ⋮\\ 1\end{array}\right]$

.

Copula-Based Multivariable Probability Calculation

It is assumed that D, S, and A are random variables, and the corresponding marginal cumulative distribution functions are denoted as FD(d), FS(s), and FA(a), respectively. The joint cumulative distribution function of drought variables based on optimal bivariate and trivariate Copula functions are as follows:

$$\begin{array}{l}F(d,s)=P\left(D\le d,S\le s\right)=C\left(F_D\left(d\right),F_S\left(s\right)\right)\\ F(d,a)=P\left(D\le d,A\le a\right)=C\left(F_D\left(d\right),F_A\left(a\right)\right)\\ F(s,a)=P\left(S\le s,A\le a\right)=C\left(F_S\left(s\right),F_A\left(a\right)\right)\\ F(d,s,a)=P\left(D\le d,S\le s,A\le a\right)=C\left(F_D\left(d\right),F_S\left(s\right),F_A\left(a\right)\right)\end{array}$$

(5)

The bivariate and trivariate joint occurrence probabilities are estimated under two distinct scenarios: (1) D and S (D and A, S and A, D and S and A) both exceed a specific threshold, denoted as $P_{DSA}^{and}$; (2) either D or S (D or A, S or A, D or S or A) exceeds a given threshold, denoted as $P_{DSA}^{or}$. The calculation formulas are as follows:

$$P_{DS}^{and}=P\left(D>d\cap S>s\right)=1-F_D\left(d\right)-F_S\left(s\right)+C\left(F_D\left(d\right),F_S\left(s\right)\right)$$

(6)

$$P_{DS}^{or}=P\left(D>d\cup S>s\right)=1-C\left(F_D\left(d\right),F_S\left(s\right)\right)$$

(7)

$$\begin{array}{c}{P}_{DSA}^{and}=P\left(D>d\cap S>s\cap A>a\right)=1-F_D\left(d\right)-F_S\left(s\right)-F_A\left(a\right)+C\left(F_D\left(d\right),F_S\left(s\right)\right)\\ +C\left(F_D\left(d\right),F_A\left(a\right)\right)+C\left(F_S\left(s\right),F_A\left(a\right)\right)-C\left(F_D\left(d\right),F_S\left(s\right),F_A\left(a\right)\right)\end{array}$$

(8)

$$P_{DSA}^{or}=P\left(D>d\cup S>s\cup A>a\right)=1-C\left(F_D\left(d\right),F_S\left(s\right),F_A\left(a\right)\right)$$

(9)

The conditional drought probability plays an important role in water resource planning and management [40]. The detailed calculation process of conditional drought probability is as follows. For example, the expression for the conditional probability distribution of S > s given D ≥ d is

$$F\left(D\ge d|S\ge s\right)=\frac{F\left(D\ge d,S\ge s\right)}{F\left(S\ge s\right)}=\frac{1-F_D\left(d\right)-F_S\left(s\right)+C\left(F_D\left(d\right),F_S\left(s\right)\right)}{1-F_S\left(s\right)}$$

(10)

Similarly, the conditional distribution of D > d given S ≥ s and A ≥ a can be expressed as

$$\begin{array}{cc}\hfill F\left(D\ge d|S\ge s,A\ge a\right)& =\frac{F\left(D\ge d,S\ge s,A\ge a\right)}{F\left(S\ge s,A\ge a\right)}\hfill \\ \hfill & =\frac{1-F_D\left(d\right)-F_S\left(s\right)-F_A\left(a\right)+C\left(F_D\left(d\right),F_S\left(s\right)\right)}{1-F_S\left(s\right)-F_A\left(a\right)+C\left(F_S\left(s\right),F_A\left(a\right)\right)}+\hfill \\ \hfill & \frac{C\left(F_D\left(d\right),F_A\left(a\right)\right)+C\left(F_S\left(s\right),F_A\left(a\right)\right)-C\left(F_D\left(d\right),F_S\left(s\right),F_A\left(a\right)\right)}{1-F_S\left(s\right)-F_A\left(a\right)+C\left(F_S\left(s\right),F_A\left(a\right)\right)}\hfill \end{array}$$

(11)

3. Results

3.1. Spatial-Temporal Variation of Drought at Multiple Time Scales

3.1.1. Temporal Evolution Characteristics of Drought at Different Time Scales

Figure 3 presents the temporal variation characteristics of the SPEI at the 1–12 month timescales in northwest China and its three climatic zones from 1960 to 2018. It can be observed that the meteorological drought was a gradual process in temporal evolution. With the increasing time scale, the drought variability declined, the inter-annual variation trend became more prominent, and the drought and wet cycles increased significantly. This implies that the frequency of drought decreased significantly, while the duration and intensity of drought increased. At different time scales, the SPEI in the plateau climate zone, westerly climate zone, and northwest China showed an ascending trend, indicating that droughts were abating (Figure 3a–c). However, the SPEI in the southeast climate zone presented a downward trend, indicating that droughts were aggravating (Figure 3d). The temporal variation characteristics of the SPEI in northwest China and the plateau climate zone were similar: a relatively severe meteorological drought occurred in around 1965; frequent droughts and floods alternated during the 1970–1990s; and drought events reduced significantly after 2005, with wet events dominating.

In the westerly climate zone, the 1970s was characterized by fluctuations in drought and flood. Drought and flood events alternated during this period. This region was dominated by drought events in the 1980s, with a severe and protracted drought from 1985 to 1988. In the 1990s, droughts significantly decreased with a decline in intensity. From the 2000s onward, drought and flood events fluctuated again, and droughts showed a deceasing trend. The drought trend in the southeast climate zone differed from that in other regions, with drought and flood events alternating until the 1990s, followed by a prolonged drought in the 1990s.

3.1.2. Spatial-Temporal Characteristics of Seasonal and Annual Drought Variation Trends

The variations in the SPEI at different time scales reflect its sensitivity to precipitation and temperature changes, as well as its cumulative effect over time. With the increasing time scale, the SPEI becomes more responsive to changes in precipitation and temperature, and short-term changes are less likely to affect its overall trend. Therefore, the SPEI is a reasonable indicator of long-term water variation.

Temporal Characteristics of Drought Variation Trend

Figure 4 illustrates the annual variation in the SPEI in the northwest region from 1960 to 2018. It can be observed that the SPEI in the northwest region, plateau climate zone, and westerly climate zone all showed a fluctuating upward trend, indicating a trend towards increased moisture, with trend rates of 0.086/10a, 0.188/10a, and 0.046/10a, respectively. In contrast, the SPEI in the southeast climate zone showed a trend towards aridity, with a rate of −0.038/10a. Additionally, the variation trend of drought was particularly distinct in the plateau climate zone (R² = 0.21).

Figure 5 shows the seasonal variation in the SPEI in the northwest region from 1960 to 2018. Obviously, the SPEI in all four seasons showed a fluctuating upward trend during this period. According to the long-term trend analysis, the seasonal drought in the northwest region was decreasing, reflected as a trend towards increased moisture in all four seasons. The trend rates for spring, summer, autumn, and winter were 0.062/10a, 0.07/10a, 0.035/10a, and 0.087/10a, respectively. The interdecadal variation in the SPEI varied significantly among different seasons. In spring, the SPEI showed an upward trend in the 1960s, followed by a sharp decline in the 1970s, another ascending trend in the 1980s, and fluctuations thereafter, reaching its minimum in 1995 (−1.46). In summer, the SPEI showed a general upward trend in the 1960s and 1970s and a slight downward trend after 1980. Since then, it has been alternating between droughts and floods. In autumn, the SPEI was roughly on the decrease in the 1960s–1980s, indicating increased aridity, but shifted to an upward trend after 1990, indicating increased moisture. No drought events occurred after 2000. In winter, the SPEI showed an upward trend before 1975, followed by a downward trend in 1975–1985. After a brief upward trend from 1986 to 1990, it turned into a rapid downward trend in the 1990s. After 2000, the SPEI fluctuated up and down, but no drought events occurred. Based on the Z-index, Zhao et al. [44] analyzed the temporal and spatial characteristics of drought in northwest China. Their results showed that the Z-index of four seasons in northwest China showed an upward trend, with alternating droughts and floods and no obvious consistency among the seasons. This was consistent with the results of this study. In conclusion, the difference in climate conditions among the four seasons led to different SPEI variation characteristics in each season. Furthermore, the seasonal variation trends of the SPEI differed significantly from its annual variation trend.

Temporal Characteristics of Drought Variation Trend

The historical and future changes in seasonal SPEI in northwest China are reported in

Table 4. Past and future variation trends of seasonal SPEI series in northwestern regions.

SPEI Series		Trend Statistics	Hurst Index	Past Trend	Persistence and Future Trend
Northwestern regions	Spring	1.73	0.54	Significant increase	Persistence: increase
	Summer	1.75	0.52	Significant increase	Persistence: increase
	Autumn	0.90	0.62	Increase	Persistence: increase
	Winter	2.09	0.59	Significant increase	Persistence: increase
Plateau climate zone	Spring	2.99	0.63	Significant increase	Persistence: increase
	Summer	3.06	0.57	Significant increase	Persistence: increase
	Autumn	1.49	0.60	Increase	Persistence: increase
	Winter	2.19	0.64	Significant increase	Persistence: increase
Westerlies climate zone	Spring	0.18	0.48	Increase	Unsustainability: decrease
	Summer	0.56	0.53	Increase	Persistence: increase
	Autumn	0.90	0.56	Increase	Persistence: increase
	Winter	0.34	0.49	Increase	Unsustainability: decrease
Southeast climate zone	Spring	0.10	0.47	Increase	Unsustainability: decrease
	Summer	0.12	0.55	Increase	Persistence: increase
	Autumn	−0.86	0.65	Decrease	Persistence: decrease
	Winter	1.77	0.53	Significant increase	Persistence: increase

. According to the table, only the autumn SPEI in the southeast climate zone showed a decreasing trend, but the trend was not statistically significant. In the future, this index is likely to continue to decline, with a potentially high intensity of drought. The seasonal SPEI in the other regions showed a rising trend. Specifically, the spring and summer SPEIs in the northwest region and the winter SPEI in the southeast climate zone exhibited a statistically significant p-value of 0.1 for their ascending trends. In addition, the winter SPEI in the northwest region and the plateau climate zone showed a significant upward trend with a p-value of 0.05. Moreover, the spring and summer SPEIs in the plateau climate zone also showed an upward trend, with a statistically significant p-value of 0.01. These ascending trends appear to continue in the future, which indicates continuous moisture. Nevertheless, the upward trends of other SPEIs were not statistically significant. The spring and winter SPEI in the westerly climate zone, as well as the spring SPEI in the southeast climate zone, are likely to decline in the future, indicating the possibility of a low-intensity drought.

The characteristic value Z for the SPEI trend of each month is shown in Figure 6. The Z-values in the northwest region from January to December were 0.81, 0.28, −0.20, −0.11, 1.19, 1.8, 0.59, −0.25, 0.52, 0.14, −0.19, and 0.73. Among them, the SPEI for March, April, August, and November presented a decreasing trend, indicating a rising trend in meteorological drought. However, the drought trend was insignificant. The SPEI for the remaining months showed an upward trend, indicating a slowing trend in drought. Furthermore, the SPEI upward trend in June passed the significance test with a p-value of 0.1. The Z-values for January, May, June, and December in the three climate zones were all greater than 0, suggesting a slowing trend in meteorological drought in these months. Additionally, the moistening trend in the plateau and westerly climate zones in June was significant, with a p-value of 0.05.

Spatial Characteristics of Drought Variation Trend

The spatial evolution characteristics of drought trends in different seasons in northwest China are shown in Figure 7. In the spring, 73.7% of the northwest region showed an upward trend in the SPEI, with a significant increase in southern Qinghai. The areas with a significant downward trend in the SPEI were mainly concentrated in the northwest of Qinghai Province, accounting for only 0.9% of the total area. In the summer, the areas with a significant upward trend in the SPEI accounted for 28.6% of the total area, mainly concentrated in Qinghai Province. However, the drought in eastern Inner Mongolia, northern Shaanxi, southern Gansu, and part of southern Qinghai showed an increasing trend, especially in southern Gansu. In the autumn, the areas with a humidification trend accounted for 71.1% of the total area, mainly concentrated in Qinghai, Inner Mongolia, and some parts of Gansu. The areas with a significant descending trend in the SPEI were concentrated in southern Shaanxi, accounting for 3.7% of the total. In the winter, the areas with an upward trend in the SPEI expanded, accounting for 83.2% of the total, and the areas with a significant humidification trend also increased. The area with a significant trend of drought accounted for only 0.8%.

In summary, the trend of droughts in different seasons in northwest China was mainly downward in the west and upward in the east, indicating a humidification trend in arid areas, while there was an aridification trend in semi-arid and semi-humid areas. From spring to winter, the humidification trend was expanding towards the east, and the areas with a significant humidification trend gradually shifted towards the east. These results indicate that the spatial distribution of droughts in northwest China tends to be uniform in the long run.

3.1.3. Spatial Variation Characteristics of Drought Intensity and Frequency

Spatial Characteristics of Drought Intensity

Figure 8 illustrates the spatial distribution of drought intensity in different seasons in northwest China. The high drought intensity areas in spring, covering 7.5% of the total area, were mostly found in southeastern Gansu, northwestern Qinghai, and some areas in Shaanxi, while the low drought intensity areas included central Inner Mongolia and western Gansu, accounting for 6.6% of the total area. The high drought intensity areas in summer accounted for approximately 8.3% of the total area. They were mainly distributed in northern Shaanxi, central and western Gansu, central Inner Mongolia, and northern Qinghai. Additionally, the maximum drought intensity was 1.49. This implied a more severe level of drought as compared to the other seasons, as well as a larger area of high-intensity drought. Nevertheless, the areas with low drought intensity accounted for only 0.9%. As for autumn, the high drought intensity areas, accounting for approximately 4.1% of the total area, were mainly concentrated in southwestern and northeastern Qinghai and southeastern Gansu. A considerable portion of Shaanxi, northwestern Gansu, and northwestern Qinghai had a low level of drought intensity, which accounted for 6.6% of the total area. The minimum drought intensity recorded was 0.91. In the case of winter, the proportions of the areas with high and low drought intensity were similar and small. The high-value areas included southeastern and northwestern Gansu and northern Qinghai, accounting for approximately 4.4% of the total area. The low-value areas consisted of southwestern Qinghai, northern Inner Mongolia, and southern Shaanxi, accounting for about 4.6% of the total area.

Spatial Characteristics of Drought Frequency

Figure 9 displays the spatial distribution characteristics of drought frequency in four seasons in Northwest China. As shown in the figure, the high-frequency areas of drought during spring were mainly in central Inner Mongolia, northwestern Gansu, and southwestern and northeastern Qinghai, with the drought frequency maintained at 25% or above. In contrast, southeastern Gansu, southeastern Qinghai, and Shaanxi remained at a low level of drought frequency, accounting for 15.9% of the total area in northwest China. During summer, the maximum drought frequency was 28.4%. The high-frequency areas mainly included eastern Inner Mongolia, southern Gansu, and southern Qinghai, covering 2.4% of the total area. The low-frequency areas were in northwestern Qinghai, central and western Gansu, and northern Shaanxi, with the minimum of 12.5%. The drought frequency during autumn drought had the maximum of 31.8% among the four seasons. The high-frequency areas mainly concentrated in western Gansu, northwestern and southeastern Qinghai, and the junction of Ningxia, Inner Mongolia, and Shaanxi. In contrast, the low-frequency areas of autumn drought were mainly distributed in central and northeastern Qinghai and southeastern Gansu. During winter, the high and low-frequency areas accounted for the smallest proportion of the total area. The areas with a high frequency of winter drought included southern Qinghai and southern Shaanxi, accounting for only 4.9% of the total area in northwest China. By contrast, the low-frequency areas were mainly concentrated in southeastern Gansu and central Shaanxi, accounting for only 13.7% of the total area.

The analysis shows that the spatial distribution of drought intensity and frequency in the study area exhibited opposite characteristics in the four seasons. For example, the southeastern part of the study area experienced high drought intensity but low drought frequency in spring, while the northeastern part experienced low drought intensity but high drought frequency in summer. Similarly, the northwestern part was characterized by low drought intensity but high drought frequency in autumn, whereas the southwestern part had low drought intensity but high drought frequency in winter. This indicates that the likelihood of high-intensity meteorological drought events occurring in the same area was relatively low, and low-intensity drought events were frequent.

3.2. Dynamic Evolution of Typical Drought Event

Temporal Evolution Characteristics of Drought at Different Time Scales

Figure 10 and Figure 11 illustrate the spatiotemporal evolution and migration path of a recent severe drought event. This drought event lasted from December 2017 to May 2018, covering a duration of 6 months. It had a severity of 2.73 × 10⁶ month·km². The drought event initially struck an approximate area of 0.31 × 10⁶ km², accounting for 18.10% of the total study area. The drought center was situated in the northwest part of Haixi Mongolian-Tibetan Autonomous Prefecture, Qinghai. In January 2018, the drought area rapidly expanded to 0.63 × 10⁶ km², accounting for 35.90% of the total area. The drought center shifted southeastwards to the central part of Haixi Mongolian-Tibetan Autonomous Prefecture with a migration distance of 191.73 km. In February of the same year, the drought area and severity declined to 0.34 × 10⁶ km² and 0.54 × 10⁶ month·km², respectively. The drought center migrated southwestwards to the northern part of the Yushu Tibetan Autonomous Prefecture at an average velocity of 255.36 km·month⁻¹. In March, the drought magnitude continued to weaken, and the drought area reached its minimum of 0.12 × 10⁶ km², accounting for only 6.90% of the total study area. It was primarily concentrated in the western part of Jiuquan City, Gansu. Subsequently, the drought magnitude showed an upward trend, with an area of 0.26 × 10⁶ km² and an intensity of 1.75, signifying the most severe level. It expanded further into central Gansu and the western part of Alxa League. In May, the drought area reduced to 0.15 × 10⁶ km², indicating a weakened drought situation and the termination of the drought. In summary, this drought event was mainly concentrated in the western part of the study area. The migration path of the drought center was characterized by a north-south oscillation pattern, and the drought had generally experienced five processes: occurrence, aggravation, mitigation, re-aggravation, and termination.

3.3. Multivariable Frequency Analysis of Drought

3.3.1. Correlation Analysis of Drought Variables

The correlation analysis of the various combination features is the foundation for selecting Copula functions and estimating parameters. In this paper, drought events with a duration of two months or longer (a total of 160 events) were selected for multi-feature frequency analysis. Figure 12 and

Table 5. Pearson, Kendall, and Spearman correlation coefficients of drought variables.

Drought Variables	Pearson	Kendall	Spearman
Duration-Severity	0.83 **	0.67 **	0.82 **
Duration-Area	0.77 **	0.57 **	0.72 **
Severity-Area	0.93 **	0.84 **	0.97 **

Note: “**” indicates significant at the 0.01 level.

, respectively, show the scatter diagram and significance test results of the correlations among drought duration, intensity, and area in the study area from 1960 to 2018. It can be observed that there was a strong linear relationship (R² > 0.6) among the drought feature variables. Moreover, each set of drought variables had high mutual dependence and a strong correlation, with a statistically significant p-value of 0.01. Therefore, it is reasonable to use the Copula function for joint distribution modeling for the frequency analysis.

As shown in Figure 12, the scatter plots of the correlation between drought duration-intensity and drought duration-area exhibit noticeable bands distributed along the vertical axis for each corresponding scale of the abscissa. This is attributed to the monthly SPEI used for drought identification, resulting in the drought duration being integer months and further a much smaller actual sample size. For example, the drought durations in the figure are only present in 2–11 months, which totals to only 10 sample values. This can affect the fitting goodness test results and parameter estimation. Furthermore, according to the Sklar theorem, the marginal distributions used for the Copula frequency analysis must be continuous functions. To achieve this, the randomization scheme proposed by De Michele et al. was adopted to transform the discrete drought duration sequence into a continuous one. This randomization method does not alter the probability distribution characteristics of the original sequence.

3.3.2. Selection of Marginal Distributions for Drought Variables

According to the K-S and A-D goodness-of-fit test results for seven candidate distribution functions, the optimal marginal distribution type for each drought feature variable was determined, and the maximum likelihood estimation was used for parameter estimation. At a significance level of p = 0.05, if the result passed the K-S test, the distribution with the smallest A-D statistic was chosen as the optimal distribution of the drought feature variable. Otherwise, the distribution was rejected directly.

Table 6. The optimal marginal distributions and its parameters of drought variables.

Drought Variables	K-S Test							A-D Statistics							Optimal Distribution	Parameters
	Gam	LogL	LogN	Wb	P-III	GEV	GP	Gam	LogL	LogN	Wb	P-III	GEV	GP
Duration	√	√	√	×	√	√	√	2.63	2.13	1.82	5.69	0.42	1.25	0.35	GP	k = −0.077 σ = 2.438 μ = 1.581
Severity	√	√	√	×	√	×	√	2.39	1.08	0.67	2.92	5.24	2.19	1.10	LogN	σ = 1.134 μ = −0.114
Area	√	√	√	√	√	√	√	0.64	1.50	1.12	0.81	0.31	1.13	0.16	GP	k = −0.199 σ = 0.561 μ = 0.034

Note: “√” indicates that it passes the K-S test, and “×” indicates that it fails the K-S test.

lists the goodness-of-fit test results and optimal distribution parameters for each distribution function of the drought feature variables, with the optimal distribution in bold black font. It can be seen that the GP distribution is the optimal marginal distribution function for the drought duration and area, while the LogN distribution is optimal for the drought intensity. Figure 13 shows the fitting effect of the empirical frequency and theoretical frequency of the optimal marginal distribution function for drought feature variables. The data points are concentrated around the 1:1 line, indicating that the theoretical and empirical probabilities of the best-fit distribution function for drought feature variables have a good correspondence.

3.3.3. Selection of Optimal Copula Functions

According to the goodness-of-fit test results, one optimal Copula function was determined for each combination of drought feature variables from six candidate Copula functions (Gumbel, Clayton, Frank, Joe, Normal, and t) to construct the joint distribution.

Table 7. The GOF tests and parameters of Copula functions.

Copula Function		Gumbel	Clayton	Frank	Joe	Normal	t	Optimal Copula	Parameters
D-S	AIC	−1195.90	−1094.77	−1235.56	−1149.28	−1181.84	−1139.16	Frank	8.01
	BIC	−1192.82	−1091.70	−1232.49	−1146.20	−1178.76	−1136.09
	RMSE	0.023	0.032	0.021	0.027	0.025	0.028
D-A	AIC	−1226.21	−1050.26	−1205.33	−1215.57	−1174.88	−1150.98	Gumbel	2.04
	BIC	−1223.14	−1047.19	−1202.26	−1212.49	−1171.81	−1147.90
	RMSE	0.021	0.037	0.023	0.022	0.025	0.027
S-A	AIC	−1231.89	−1300.85	−1270.51	−1061.98	−1266.06	−1259.59	Clayton	10.63
	BIC	−1228.82	−1297.77	−1267.43	−1058.90	−1262.99	−1256.51
	RMSE	0.021	0.017	0.019	0.036	0.019	0.02
D-S-A	AIC	−1204.64	−1056.04	−1202.06	−983.64	−1147.71	−1140.11	Gumbel	3.62
	BIC	−1201.57	−1052.97	−1198.99	−980.56	−1144.63	−1137.04
	RMSE	0.022	0.037	0.023	0.046	0.028	0.029

lists the goodness-of-fit test results and the parameters of the optimal Copula functions. Smaller AIC, BIC, and RMSE values indicate better fitting results. The optimal Copula function is in bold black font. As shown in the table, the Frank Copula function has the smallest AIC, BIC, and RMSE values for the duration-intensity combination, making it the optimal joint model for this combination. The Gumbel Copula function performed the best in constructing the joint probability distribution of the duration-area combination, while the optimal model for the intensity-area combination was the Clayton Copula function. The optimal joint probability distribution model for the duration-intensity-area combination was constructed using the Gumbel Copula function.

Figure 14 shows the PP (probability-probability) of the theoretical cumulative probability of the optimal Copula function and the empirical Copula cumulative probability for each combination of drought feature variables. The data points of the empirical and theoretical cumulative probabilities of the joint distribution of each group of drought feature variables are evenly distributed on both sides of the diagonal line, indicating that the results of the two are roughly equal. The optimal joint distribution function can better describe the joint probability distribution between drought feature variables.

3.3.4. Joint Occurrence Probability of Drought

The joint occurrence probabilities between drought variables can be divided into two situations. One is when a variable, such as drought duration, severity, or area, is greater than or equal to a specific value (referred to as the “or” situation), such as P (D ≥ 5 ∪ S ≥ 2) and P (D ≥ 5 ∪ S ≥ 2 ∪ A ≥ 1). The other is when all the three drought variables simultaneously exceed a specific value (referred to as the “and” situation), such as P (D ≥ 5 ∩ S ≥ 2) and P (D ≥ 5 ∩ S ≥ 2 ∩ A ≥ 1). The bivariate and trivariate joint occurrence probabilities of drought variables are shown in Figure 15. It can be seen that the range of high joint occurrence probability in the “or” situation is much larger than that in the “and” situation. The joint occurrence probability of the same combination of drought variables is higher in the “or” situation. For example, when the drought duration and severity are greater than or equal to 4.09 months and 11.02 × 10⁶ month·km², respectively, the joint occurrence probability in the “or” situation is 49.4%, while it is 30.2% in the “and” situation. In both situations, the joint occurrence probability increases as the drought variables decrease. For example, when the drought duration and severity exceed 6.65 months and 2.73 × 10⁶ month·km², respectively, in the “or” situation, the joint occurrence probability is 20.1%; when they exceed 3.30 months and 1.97 × 10⁶ month·km², respectively, the probability is 49.8%.

3.3.5. Conditional Probability of Drought

The conditional probabilities of meteorological droughts can be easily obtained based on the joint probabilities of droughts calculated using the Copula function. Figure 16 presents the two-dimensional conditional probability among drought duration, severity, and area. Figure 16a shows the conditional probability of drought severity when the drought duration exceeds a certain threshold. Figure 16b shows the conditional probability of drought duration when the drought intensity is larger than a certain threshold. The results indicate that with increasing conditional factors, both the conditional probabilities of drought severity and duration show a rising trend. Specific probabilities can be determined for a given conditional factor from the figure. For instance, given a drought duration of more than 5 months and drought severity of more than 1 × 10⁶ month·km², the probability of drought is 93.4%. It means that there is a high likelihood of meteorological drought under such conditions. Furthermore, if the drought severity is larger than 1, 3, 5, 7, or 9 × 10⁶ month·km² and the drought duration is greater than 5 months, the probability of meteorological drought will be 46.2%, 73.9%, 79.8%, 81.7%, and 82.6%, respectively. The distribution characteristics of probability show a dense distribution at both ends and a sparse distribution in the middle, indicating the difference in probability results at different conditional factor levels. For example, when the drought severity is high (Figure 16a), the probability results for different conditional factors have little difference. Additionally, an increase in the conditional factor (drought duration) will significantly raise the probability of meteorological drought events with moderate intensity, similar to the results of Yusof et al. [45].

Figure 17 presents the conditional probability of meteorological drought for three drought variable combinations. As shown in Figure 16 and Figure 17, the conditional probability of meteorological drought under three drought variables is greater than that under two. For example, in the case of a drought severity above 5 × 10⁶ month·km² and a drought duration longer than 5 months, the conditional probability of meteorological drought is 79.8%. Meanwhile, given a drought area greater than 1 × 10⁶ km², this figure increases to 99.8%. It can be concluded that ignoring any drought variable may significantly underestimate the occurrence probability of severe meteorological droughts, as well as the severity of drought events [17].

4. Discussion

This study found that the temporal variation in the drought index represented an ascending trend in Northwestern China due to a continuous increase in precipitation and temperature [46].Especially, this phenomenon was significant in arid regions, which indicates that runoff and glacier ablation continued to increase, and the water level of inland lakes increased significantly, the rise and flood disasters increased, while the vegetation improved and the number of sandstorm days decreased [47,48]. Moreover, the warm-wet tendency in northwest China exhibits an obvious response to the “global warming stagnation” [49].

The present climate in northwest China has been changing from warm-dry to warm-wet [50], and the 1980s was the turning point of drought. Taking the annual scale drought index series as an example, the moving T test is used to detect the change point of drought in the study area. The step is set to 6, for a fixed significance level of p = 0.05 in the present paper, and the corresponding critical value is 2.31 (Figure 18). Obviously, the statistical value exceeded the reliable lines of the significance level (p = 0.05) between 1981 and 1985. Therefore, the change points of the annual drought index were concentrated between 1981 and 1985, and the property (positive or negative) of the maximum change point (in 1993) could be used to estimate the main trend of the drought index series.

Similar studies have pointed out the late 1980s were the abrupt wetting point of drought. For example, Wang et al. [51] found that there were two evident contrasting periods of drought between 1961 and 2010. Sustaining drought conditions existed between 1961 and 1986, whereas wet conditions dominated after 1986. Li et al. [52] reported a weakening drought in Northwestern China based on the increasing trends of the SPI series. Additionally, as the variation of annual precipitation in northwest China showed a significant upward trend in the central and west of the study region, and a downward trend in the eastern regions [53,54], the SPEI exhibited a fluctuating upward trend in the plateau climate zone and westerly climate zone, while there was a trend towards aridity in the southeast climate zone.

In summary, joint occurrence probabilities can quantitatively capture the occurrence probability of meteorological drought under different drought variables, providing important information for drought mitigation and water resource management. For example, the occurrence probability of meteorological drought with a certain duration and intensity can serve as a critical condition for identifying specific water supply systems and triggering drought emergency plans [55]. Overestimation of drought risk may occur if only the “or” situation is considered, while underestimation may occur if only the “and” situation is considered. Therefore, a comprehensive analysis of both situations is necessary to ensure an accurate assessment of the occurrence probability of meteorological drought. Additionally, the conditional occurrence probability is useful for assessing water resource systems under extreme drought conditions and determining if the systems can meet the water demand under given drought conditions. This information allows water resource managers to quickly determine the amount of additional water needed to alleviate drought.

5. Conclusions

In this study, the SPEI calculated based on the gridded data from 1960 to 2018 was used to estimate meteorological drought. Drought events were identified through a 3D identification method, showing the dynamic evolution of drought events. The optimal marginal distribution of drought variables was selected from seven alternative distributions. Moreover, the best-fitted Copula was selected from six candidate Copulas using the RMSE, AIC, and BIC methods. These selected Copulas were employed to estimate drought risk through bivariate and trivariate joint probability and conditional probability analyses, respectively. The major conclusions are as follows:

(1)

Overall, the SPEI showed an upward trend in the plateau climate zone and the westerly climate zone, with rates of 0.188/10a and 0.046/10a, indicating a mitigation of drought. Conversely, the SPEI demonstrated a descending trend (−0.038/10a) in the southeast climate zone, suggesting an intensified drought situation.

(2)

The variation trend of drought in different seasons was mainly downward in the west part and upward in the east part of Northwestern China. From spring to winter, the humidification trend expanded towards the east, and the areas with a significant humidification trend gradually shifted towards the east.

(3)

The spatial distribution of drought intensity and frequency in different seasons exhibited opposite characteristics. For example, the southeastern part of the study area experienced a high drought intensity but a low drought frequency in spring. This indicated that the likelihood of high-intensity meteorological drought events occurring in the same area was relatively low, whereas low-intensity drought events were frequent.

(4)

The most severe drought event occurred from January 1961 to October 1962 and experienced five processes: occurrence, aggravation, mitigation, re-aggravation, and termination. The migration path was characterized by north-south oscillation.

(5)

The joint occurrence probabilities were consistently higher in the “or” situation than in the “and” situation for the same combination of drought variables. Furthermore, the conditional probability of drought variables, given specific conditional factors, declined as the values of these factors increased. Notably, as the conditional factors increased, there was a noticeable reduction in the drought occurrence probability for drought variables with lower values.