Nonstationary Spatiotemporal Projection of Drought Across Seven Climate Regions of China in the 21st Century Based on a Novel Drought Index

1. Introduction

Drought is a frequent natural variability phenomenon that develops slowly, persists for long periods, and affects wide areas. It damages water supply, agricultural production, ecosystems, and socioeconomic development to varying degrees; nearly every country has experienced drought[1,2]. For example, the once-in-50-year extreme drought in Russia in 2010 seriously threatened local food security[3]. The severe East African drought of 2011—the worst in 60 years—caused food shortages and affected 12.4 million people[4]. Situated in East Asia, China is prone to natural hazards; meteorological disasters account for about 70% of all natural disasters, with droughts comprising over half of meteorological events. After 2000, extreme droughts have become more frequent. The 2019 drought in the middle and lower Yangtze River basin affected seven provinces and severely impacted agriculture[5]. In the summer of 2022, China was hit by its most severe heatwave in six decades, exacerbating a drought that has impacted food and factory production, power supplies, and transport in a vast area of the country[6]. The IPCC Sixth Assessment Report (AR6) identifies greenhouse gas emissions as the main driver of warming. The global mean temperature for 2001–2020 was 0.99 °C (0.84–1.10 °C) higher than 1850–1900. Global warming has become undeniable, weakening climate-system stability and altering the water cycle. Human activities—excessive water use, land-surface changes, and other disturbances—further disrupted hydrological systems, increasing drought intensity and frequency, exacerbating regional water supply–demand conflicts, and causing water-quality deterioration, crop losses, and ecological degradation. Extreme droughts directly threaten national food security and socioeconomic stability, inflicting substantial social and economic losses[7,8,9,10].

In addition, climate change enhances uncertainty in drought projections. Many studies used climate model outputs as inputs to offline drought indices or hydrological-impact models and found amplified future drought frequency and intensity[11,12]. However, these projected increases conflict with observations showing rising global vegetation cover and little change in runoff over recent decades, and with projections indicating slight future increases in vegetation cover and runoff[13,14].

This contradiction arises partly from differences among drought metrics. Another factor is the method used to compute atmospheric water demand. Sheffield, et al. [15] showed that temperature-based potential evapotranspiration (PET) algorithms (e.g., Thornthwaite) tend to overestimate PET and thus overstated drought. Some studies also indicated that ignoring the effect of rising CO₂ on vegetation physiology during evapotranspiration calculation contributed to overly dry projections[14]. Therefore, future drought projections should evaluate how different drought indices affect results to select those appropriate for climate-change scenarios, and should fully account for CO₂-induced physiological effects on vegetation, PET estimates, and drought assessment to produce more realistic forecasts and better inform early warning and drought resilience measures.

Drought characteristics, such as severity, duration, and areas, are key indicators for drought evaluation. Current studies on future drought projections often break the spatiotemporal continuity of drought characteristics, limiting analysis to temporal or spatial aspects separately and reducing high-dimensional drought events to low-dimensional problems. This simplification causes substantial information loss and prevents accurate representation of drought structure in space and time. Moreover, frequency analyses of multiple drought characteristic variables typically assume these variables are stationary, but under intensifying climate change, this assumption no longer holds. Therefore, a nonstationary framework for multidimensional frequency analysis and risk assessment of drought characteristics is required[16].

This study aims to address the aforementioned issues by developing a CO₂-aware Standardized Moisture Anomaly Index (SZI[CO₂]). It employs a three-dimensional drought identification method to extract key drought characteristics—duration, area, and severity—and analyzes their spatiotemporal dynamics. Additionally, the study utilizes Copula-based approaches with time-varying parameters to estimate joint drought probabilities and project the spatiotemporal evolution of droughts under SSP2-4.5 and SSP5-8.5 scenarios.

2. Materials and Methods

2.1. Study Regions

China is vast with complex terrain—higher in the west and lower in the east—where mountains, plateaus, and hills cover two-thirds of the land and plains and basins less than one-third. Climates vary from south to north, including tropical monsoon, subtropical monsoon, temperate monsoon, temperate continental, and alpine plateau climates. Regional climate change affects drought differently, so climatic zoning is needed to explore future drought characteristics by region. Based on existing climate classifications[13], China is divided into seven climate regions (Figure 1).

2.2. Historical Data

Observed monthly precipitation and daily maximum and minimum temperature data used in this chapter were obtained from the Climate Change Research Center (https://ccrc.iap.ac.cn/resource/detail?id=228). This dataset is based on the latest compilation of observations from 2,472 ground stations by the National Meteorological Information Center, interpolated using spline methods and possessing high observational accuracy[17,18,19]. Monthly wind speed, pressure, relative humidity, and radiation data were obtained from the ERA5 reanalysis dataset (https://cds.climate.copernicus.eu/datasets/reanalysis-era5-single-levels-monthly-means?tab=download), covering 1940 to present at 0.25°×0.25° resolution. The above data were resampled to the spatial resolution of the climate models using bilinear interpolation for comparison with CMIP6 model outputs.

2.3. Global Climate Models Data

Monthly precipitation and temperature datasets from 10 GCMs in CMIP6 were selected (Table 1). The CMIP6 historical period is 1850–2014, and the future-scenario period is 2015–2100.

2.4. Downscaling GCMs Data

Climate model output data typically have coarse spatial resolution and significant simulation biases, necessitating bias correction and downscaling procedures. The Bias Correction and Spatial Disaggregation (BCSD) method, introduced by Wood, et al. [20], is a widely used statistical downscaling technique that has been extensively applied in GCM data downscaling studies[21], achieving good results. Therefore, the BCSD method was used to downscale GCMs’ data in this study.

2.5. Drought Index Calculation

To calculate SZI, we first need to estimate PET. Yang, Roderick, Zhang, McVicar and Donohue [14] used meteorological and hydrological data and CO₂ concentration outputs from GCMs in CMIP5 to derive a PET calculation formula that accounts for CO₂ concentration (PM[CO₂]) based on the relationship between canopy resistance parameter (r_s) and CO₂ concentration:

$$PET={\displaystyle \frac{0.408\mathsf{\Delta }\left(R_n-G\right)+\frac{900}{T+273}{u}_2\left(e_s-e_a\right)}{\mathsf{\Delta }+\mathsf{\gamma }\{1+u_2\left[0.34+2.4\times {10}^{-4}\left(\left[C{O}_2\right]-300\right)\right]\}}}$$

(4)

where R_n is the net radiation (MJ m^-2 day^-1), G is the soil heat flux (MJ m^-2 day^-1), e_s and e_a represent saturation and actual vapor pressure of the air (kPa), T is the mean daily air temperature at 2m height (°C), u₂ is the wind speed at 2m height (m s^-1), $\mathsf{\Delta }$ is the slope vapor pressure curve (kPa °C^-1), $\mathsf{\gamma }$ is the psychrometric constant (kPa °C^-1), and [CO₂] is the CO₂ concentration (ppm).

Then, the monthly climatically appropriate precipitation was calculated:

$\widehat{P}={\mathsf{\alpha }}_{j}PET+{\mathsf{\beta }}_{j}PR+{\mathsf{\gamma }}_{j}PRO-{\mathsf{\delta }}_{j}PL,\hspace{1em}j=\mathrm{1,2},3,\dots ,12$ (5)

The right-hand side of the equation represents, respectively, the climatically appropriate evapotranspiration, soil moisture recharge, runoff, and soil moisture loss — meaning that, under climate-appropriate conditions in a given month, precipitation must supply water for evapotranspiration, for replenishing soil moisture, and for generating runoff. In addition, soil moisture will appropriately lose some water to evapotranspiration and runoff, which reduces the demand on precipitation and therefore must be subtracted as the corresponding climatically appropriate soil moisture loss. α_j, β_j, γ_j, and δ_j denote the coefficients for the water-balance components corresponding to different months:

$${\mathsf{\alpha }}_j={\displaystyle \frac{{\overline{ET}}_j}{PE{T}_j}},\hspace{1em}{\mathsf{\beta }}_j={\displaystyle \frac{{\overline{R}}_j}{P{R}_j}},\hspace{1em}{\mathsf{\gamma }}_j={\displaystyle \frac{{\overline{RO}}_j}{PR{O}_j}},\hspace{1em}{\mathsf{\delta }}_j={\displaystyle \frac{{\overline{L}}_j}{P{L}_j}}$$

(6)

where ${\overline{ET}}_j$, ${\overline{R}}_j$, ${\overline{RO}}_j$, and ${\overline{L}}_j$ represent the multi-year mean values of the actual evapotranspiration, actual soil moisture recharge, actual runoff, and actual soil moisture loss for month j, respectively. PET_j, PR_j, PRO_j, and PL_j are the potential evapotranspiration, potential soil moisture recharge, potential runoff, and potential soil moisture loss, respectively, calculated using a two-layer soil model. After obtaining $\widehat{P}$, the water deficit d can be calculated as:

$$d=P-\widehat{P}$$

(7)

where P represents monthly precipitation. The empirical distribution was used to estimate the cumulative probability for d:

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

(8)

where n denotes the number of samples, i denotes the ranked position of the variable, and p denotes the non-exceedance probability of variable x in month t. Then, the empirical probability values are converted to standardized values, that is SZI[CO₂], using the inverse Gaussian distribution:

$SZI\left[{CO}_2\right]={\varphi }^{-1}\left(p\right)$ (9)

2.6. Three-Dimensional Drought Characteristics Identification

Drought is essentially a spatiotemporally continuous phenomenon. The run theory–based drought identification process cannot simultaneously account for both temporal and spatial dimensions of drought, so it is necessary to use spatiotemporal identification methods to extract spatiotemporally continuous drought events[22,23,24]. In this study, a 3-D array of SZI[CO₂] (longitude–latitude–time) was computed from gridded meteorological data; each grid cell is represented as SZI[CO₂](i, j, k), where i denotes longitude, j denotes latitude, and k denotes the time index. The following steps were used here to identify three-dimensional drought events:

(1) Delineation of drought patches. First, for each monthly 2-D drought-index grid, we identify cells below the drought threshold (SZI[CO₂] < −1). Then, using a 3×3 neighborhood to group spatially adjacent drought cells: cells that are drought-affected in adjacent positions are grouped and assigned a common identifier, merged into a single drought patch (Figure 2). If a drought cell has no adjacent drought neighbors, assign a new identifier for a new patch; repeat until all cells for that month are processed. This yields multiple drought patches in different areas. Apply a given area threshold (A₀) to screen patches: patches larger than A₀ are defined as a drought event (Figure 2 A₁, A₂, and A₄), while patches smaller than A₀ are discarded (Figure 2 A₃ and A₅). A₀ is also used to determine temporal continuity between patches, preventing the merging of unrelated or weakly related drought events across adjacent months[25].

(2) Temporal connection of drought patches. After identifying monthly drought patches, we need to determine whether patches in adjacent months are connected and can form a single drought event. Let a patch at time t have area A_t, the corresponding patch at t+1 have area A_t+1, and their overlapping area be A*. If A* > A₀, A_tand A_t+1 are considered temporally continuous and belong to the same drought event (A₃ and A₄); otherwise, they are treated as two independent drought events (A₁ and A₂) (Figure 3). Following this rule, evaluate the overlap A* between patches at successive times; when A* < A₀, the drought event is deemed to have ended, and the patches identified as the same event are assigned the same identifier. Repeat this process to link all patches in time, producing 3-D continuous drought bodies and yielding multiple 3-D drought events.

In 3-D drought identification, the drought-area threshold A₀ is a key parameter. Its selection strongly affects identification results: if A₀ is too small, unrelated droughts may be linked into a single event; if A₀ is too large, some moderate droughts may be excluded, causing misclassification. Studies show that setting A₀ to about 1.6% of the study area is appropriate [23].

Extraction of 3-D drought event characteristics can more comprehensively reflect the spatiotemporal continuity of drought. This study extracted the drought duration, area, severity, intensity, and centroid of 3-D drought events. For more details, one can refer to Feng et al. [23].

2.7. Non-Stationary Frequency Analysis

The fitting of cumulative distribution functions to drought characteristic variables is the foundation of drought characteristic frequency analysis. However, there’s significant uncertainty regarding the types of distribution functions that different drought characteristic variables follow. Therefore, it’s necessary to screen the distribution types to determine the optimal marginal distribution function for each drought characteristic variable. We selected seven marginal probability distribution functions: Normal (NO), Log-Normal (LNO), Gamma (GAM), Weibull (WB), Exponential (EXP), Log-Logistic (LOG), and Generalized Extreme Value (GEV) distributions.

As future climate change intensifies, the assumption of stationarity for drought characteristic variables may no longer hold. Therefore, it’s necessary to use non-stationary methods to perform frequency analysis of the marginal distributions. In recent years, Generalized Additive Models for Location, Scale, and Shape (GAMLSS)[26] have been widely used in non-stationary frequency analysis research. The GAMLSS model is a semi-parametric regression model that can flexibly simulate the changes in distribution parameters with covariates. In this paper, the time variable was selected as the covariate, and the linear relationship between the location parameter μ and the covariate can be expressed as ${\mu }_t=b_0+b_{1}t$. The relationship between the scale parameter σ and the covariate is consistent with the expression of the location parameter. This paper considers one stationary model (M0) and three non-stationary models: a time-varying scale parameter model (M1), a time-varying location parameter model (M2), and a time-varying location-scale parameter model (M3).

This paper selected the Copula model to analyze the joint characteristics of multivariate drought. The Copula function can connect marginal distribution functions following different distributions to perform multivariate dependency measurement and joint probability calculation. In recent years, the Copula function has been widely used in multivariate hydrological frequency calculation. By combining non-stationary marginal distribution fitting and Copula, the optimal marginal distribution function of different drought characteristic variables and the optimal Copula combination of different drought characteristic variables can be obtained. Parameter estimation was performed using the maximum likelihood method. After determining the optimal parameters, the goodness of fit is tested using the K-S method. Then, the AIC and BIC criteria were used to optimize the selection of the Copula function, in order to determine the optimal Copula function between different characteristic variables. Furthermore, the joint probability and conditional probability of different variable combinations can be calculated with Bayes’ theorem and the law of total probability.

3. Results

3.1. Temporal Trends of SZI[CO₂] in Different Regions of China

From 1985 to 2100, under the SSP2-4.5 and SSP5-8.5 scenarios, Central China (Figure 4g) and South China (Figure 4h) show a pronounced decline in SZI[CO₂], with values decreasing from about 0.3 in the 1980s to −0.3 in the 2090s, at an average rate of −0.005 per year, indicating a clear drying trend in the southern humid regions under future climates. Other climate zones (Northwest, Inner Mongolia, the Tibetan Plateau, Northeast, and North China) exhibit varying degrees of wetting, with SZI[CO₂] increasing in future scenarios. Northwest (Figure 4b) and Inner Mongolia (Figure 4c) have the highest increase rates, about 0.010 per year. The Tibetan Plateau (Figure 4d) shows a higher increase under SSP5-8.5 (0.008 per year) than under SSP2-4.5 (0.005 per year). Northeast (Figure 4e) and North China (Figure 4f) have lower increase rates, around 0.005 per year. Because most of the country trends toward wetter conditions, the national-scale SZI[CO₂] also increases, at a rate of 0.004 per year under SSP2-4.5 and 0.005 per year under SSP5-8.5.

The trends in relative drought area (RDA) across climate regions are similar to those of SZI[CO₂]. Northwest (Figure 5b) shows the fastest RDA decline: 0.0019/yr under SSP2-4.5 and 0.0015/yr under SSP5-8.5. Next is Inner Mongolia (Figure 5c), with RDA change rates of −0.0014/yr (SSP2-4.5) and −0.0020/yr (SSP5-8.5). On the Tibetan Plateau (Figure 5d), RDA declines at 0.0008/yr (SSP2-4.5) and 0.0012/yr (SSP5-8.5). North China (Figure 5f) shows declines of 0.0017/yr (SSP2-4.5) and 0.0010/yr (SSP5-8.5). Northeast (Figure 5e) has the slowest decline, about 0.0010/yr. Central and South China (Figure 5g,h) show increasing RDA; under SSP5-8.5, the increase rate is 0.0023/yr, and under SSP2-4.5, the rates are 0.0018/yr and 0.0015/yr, respectively. Overall, drought area increases in southern humid regions and decreases in northern regions. At the national scale (Figure 5a), RDA declines at 0.005/yr, indicating a marked reduction in future drought area.

3.2. Spatiotemporal Dynamics of a Representative Future Drought Event

Drought events are inherently spatiotemporally continuous structures, so it is necessary to analyze their 3-D (longitude–latitude–time) evolution from a spatiotemporal continuity perspective. This section used 3-month SZI[CO₂] and the 3-D drought identification method to detect droughts in seven climate regions under two climate scenarios, in order to analyze future spatiotemporal dynamics and development patterns of droughts in different regions.

Using the projected 2094–2097 drought in Central China under SSP5-8.5 as an example, we analyzed its spatiotemporal evolution and time series of drought metrics. The event is projected to start in August 2094 and end in May 2097, lasting 34 months, with a severity of 5.2×10⁷ km²·month (Figure 6). The drought initiates in the southwest of Central China with an area of 43,000 km² (2.2% of the region). By October 2094, it rapidly spreads to the central–western area (Figure 7), reaching 780,000 km² (40% of the region) and a severity of 1.3×10^6 km²·month, peaking in November. It then weakens, re-expands in July 2095, and reaches a second peak in August 2095 (covering 75% of the region, severity 3.1×10⁶ km²·month). The drought would undergo further cycles of decline–expansion–peak, reaching its overall maximum in Jan–Mar 2096; in February 2096 drought area would peak at 93% and the severity reaches 4.4×10⁶ km²·month. After this peak, the drought gradually weakens; after April 2096 affected area falls below 50%, and by May 2097, the drought would end. As shown in Figure 7, the event starts in southwest Central China, spreads to the central area, develops into a region-wide drought, then shifts northeast before disappearing. The drought centroid is mainly in the central–western area, where severity is markedly higher than elsewhere; five months have intensity > 2.0. The event lasts 34 months and undergoes six peaks, with February 2096 being the most severe (intensity 2.4). Two minor peaks occur in November 2096 and March 2097, but each affects less than 50% of the region and has low intensity.

3.3. Frequency Analysis of Multiple Drought Characteristic Variables

Using 28 stationary/nonstationary marginal distributions formed from 7 marginal distribution families and 4 parameter-model combinations, we fitted marginal distributions to drought duration, area, and severity for the seven regions under the two climate scenarios, and selected the best models using AIC and BIC criteria. The LON distribution was selected most often, followed by GEV, GAM, and WB.

Among the four parameter models, the M2 model — which only varies the location parameter — was selected most frequently. The M2 model represents a distribution whose location parameter changes over time while the scale parameter remains fixed; the location parameter describes the series mean. Under future climate scenarios, the trend term is the main deterministic component of drought characteristics, so a nonstationary distribution with a time-varying mean best represents this type of nonstationary time series.

Identifying the marginal distribution types of drought variables provides the basis for building a Copula-based probability model. The AIC- and BIC-selected marginal distributions were matched to the drought variables; parameters were estimated by maximum likelihood, and fits were tested by the K–S test. All selected models passed at the p = 0.05 significance level.

The joint occurrence probabilities among drought characteristic variables include two types. One is the “and” case, representing the situation where drought duration, area, or severity simultaneously exceed certain thresholds, e.g., P(D > 5 ∩ A > 1×10⁵ km²). The other is the “or” case, indicating that any one of the variables exceeds a threshold, e.g., P(D > 5 ∪ A > 1×10⁵ km²).

Figure 8 shows the joint distribution of drought duration and area under both “and” and “or” cases across the seven climatic regions in China under the SSP2-4.5 and SSP5-8.5 scenarios. It can be observed that the range with higher joint occurrence probability under the “or” case is much larger than under the “and” case. For example, in Northwest China under SSP2-4.5, the joint probability of drought duration exceeding 4 months and area exceeding 4×10⁵ km² is 21% under the “or” case, and 13% under the “and” case. In both cases, the joint probability increases as the thresholds for drought characteristics decrease. Continuing with Northwest China under SSP2-4.5 as an example, the joint probability of drought duration longer than 6 months and area larger than 5×10⁵ km² is 10% under “or,” while the joint probability of duration longer than 4 months and area larger than 3×10⁵ km² reaches 31%. Similar patterns are observed in other regions: the joint probability for the same variable combinations under the “or” case is higher than under the “and” case, and the joint probability increases as the drought characteristic thresholds decrease.

The joint occurrence probabilities of drought duration–severity and drought area–severity show similar patterns with respect to the variation of drought features as those observed for drought duration–area. The three-dimensional joint probabilities of drought duration, area, and severity are shown in Figure 9. When the drought characteristic variables are held at the same fixed values, the joint probability under the “or” case is higher than under the “and” case. Taking Northwest China under SSP2-4.5 as an example, the joint probability for drought duration exceeding 3.2 months, drought area exceeding 3×10⁵ km², and drought severity exceeding 1.2×10⁶ km²·month is 47% under the “and” case, and 21% under the “or” case (Figure 9).

4. Discussion

In statistics, stationarity is a special type of stochastic process where statistical features such as mean and variance do not change over time. Conversely, non-stationarity refers to processes where these statistical properties vary with time[27]. The assumption of stationarity is foundational for fitting time series distributions; however, with increasing climate change and human activities, the stationarity assumption for drought characteristic variables is often violated. Therefore, when performing frequency analysis, the influence of time variables on parameters is considered. Using time or large-scale circulation factors as covariates for parameters are two common approaches. Choosing time as a covariate does not require evaluating or screening large-scale circulation indices (such as El Niño–Southern Oscillation, Indian Ocean Dipole, etc.). Moreover, obtaining large-scale circulation factors from climate model data is challenging and uncertain. As a result, in non-stationary frequency analysis of future drought characteristics, it is difficult to apply large-scale circulation factors as covariates broadly.

The Archimedean Copula is currently the most widely used type of Copula in hydrology. It features a simple form, good flexibility, and symmetry. Copulas can effectively model the joint distribution of two variables for bivariate analysis, but extending to high dimensions presents challenges, particularly in parameter estimation. In this paper, we used the three-dimensional Copula formula and employed maximum likelihood estimation to directly estimate parameters, thereby obtaining the high-dimensional joint distribution. This method is suitable for single-parameter Copulas; however, for multi-parameter forms, directly estimating parameters becomes more difficult. Hence, methods such as nested Copulas or mixture Copulas can be used for parameter estimation. Each step of nested or mixture Copulas involves estimating parameters for bivariate Copulas, which avoids the difficulties associated with directly estimating parameters for high-dimensional multivariate Copulas. Additionally, because various types of bivariate Copulas can be used within this framework, with numerous options available, these methods can achieve good fitting results.

5. Conclusions

This study used downscaled climate model data to calculate the SZI[CO₂] drought index under China’s SSP2-4.5 and SSP5-8.5 scenarios. A three-dimensional spatiotemporal drought identification method was employed to extract drought events and examine their dynamic spatial-temporal evolution. A non-stationary drought characteristic frequency analysis framework was constructed, utilizing Copula functions for multivariate frequency analysis of drought features. The main conclusions are as follows:

(1) In Northwest China, Inner Mongolia, the Tibetan Plateau, Northeast China, and North China, SZI[CO₂] shows a significant increasing trend, while the affected drought area exhibits a decreasing trend, indicating a future wetting trend in these regions. Conversely, Central China and South China show signs of becoming drier, with increases in drought frequency, duration, and severity.

(2) Drought characteristics (duration, area, severity) identified through the three-dimensional method display obvious trend components. Comparative analysis of seven stationary and non-stationary marginal distributions reveals that non-stationary LON and GEV distributions are suitable for modeling the frequency distributions of drought features in most regions.

(3) When the same values of drought features are considered, the joint occurrence probability of drought under SSP5-8.5 is higher than under SSP2-4.5. Regions with notable differences include the Tibetan Plateau, Central China, and South China. The conditional probability of drought occurrence considering three features is significantly higher than with only two features, indicating that ignoring any one drought characteristic is likely to underestimate the probability of severe drought events.