Copula-Based Assessment and Regionalization of Drought Risk in China

Droughts are widespread in China and have brought considerable losses to the economy and society. Droughts are intricate, stochastic processes with multi-attributes (e.g., duration, severity, intensity, and return period). However, most drought assessments tend to focus on univariate drought characteristics, which are inadequate to describe the intrinsic characteristics of droughts due to the existence of correlations between drought attributes. In this study, we employed the standardized precipitation index to identify drought events using China’s monthly gridded precipitation dataset from 1961 to 2020. Univariate and copula-based bivariate methods were then used to examine drought duration and severity on 3-, 6-, and 12-month time scales. Finally, we used the hierarchical cluster method to identify drought-prone regions in mainland China at various return periods. Results revealed that time scale played an essential role in the spatial heterogeneity of drought behaviors, such as average characteristics, joint probability, and risk regionalization. The main findings were as follows: (1) 3- and 6-month time scales yielded comparable regional drought features, but not 12-month time scales; (2) higher drought severity was associated with longer drought duration; (3) drought risk was higher in the northern Xinjiang, western Qinghai, southern Tibet, southwest China, and the middle and lower reaches of the Yangtze River, and lower in the southeastern coastal areas of China, the Changbai Mountains, and the Greater Khingan Mountains; (4) mainland China was divided into six subregions according to joint probabilities of drought duration and severity. Our study is expected to contribute to better drought risk assessment in mainland China.


Introduction
Drought is one of the most severe climate-related threats to human civilization, with randomness, creep, and complexity [1,2]. It refers to a temporary lack of precipitation or water shortage in a specific period, irrespective of the region's typical aridity. Consequently, droughts can take place in both wet and dry regions [3,4]. Drought can generally be defined from various aspects such as meteorology, hydrology, agriculture, and socioeconomics. Currently, droughts are commonly monitored and analyzed through various drought indices obtained from gauge-based or remote-sensing data. The widely-used climate-related drought indices include the Palmer Drought Severity Index (PDSI) [5], Standardized Precipitation Index (SPI) [6], Standardized Precipitation Evaporation Index (SPEI) [7], and so on. Our research is limited to meteorological drought as evaluated by SPI, which has been proven to be an effective index for assessing drought severity and duration from regional to global scales [8][9][10][11][12]. The World Meteorological Organization and a few international meteorological drought centers have also suggested the SPI as one reference drought index to track drought features, due to its easy computation, multi-time scale properties, and spatiotemporal comparability [13,14].

Data Source and Collection
In this study, the utilized monthly gridded precipitation dataset (version 2.0) has a spatial resolution of 0.5 • × 0.5 • and encompasses a time span ranging from 1961 to 2020. The dataset was acquired from the China Meteorological Data Service Centre [35]. Thin-plate smooth spline interpolation was employed to construct the gridded dataset based on the precipitation observations from 2472 meteorological stations across China [36]. Excluding Taiwan and South China Sea Islands, the dataset comprises 4189 grids and was verified by average root-mean-square and cross-validated by the Tropical Rainfall Measurement Mission product [36,37]. As a result, the dataset has an average error of 0.49 mm/month, and the average correlation coefficient between the grids and the observations is 0.93 (arriving at a significance level of α = 0.01) [38]. The dataset has also been extensively applied for research on climate change in China [1,30,39]. The multi-year average precipitation was calculated using this dataset and is shown in Figure 1. The spatial pattern of annual precipitation indicates a gradual decrease from the southeastern coast to the northwest inland. In addition, the seasonal distribution of precipitation is also uneven, which mainly falls in May-October in the form of high-intensity rainforms [40,41]. The uneven spatial distribution of precipitation and the considerable disparities in the time domain have led to frequent droughts and floods in China, which have greatly threatened agricultural production, as well as people's lives and properties. (3) to perform traditional univariate and copula-based bivariate frequency analyses on drought duration and severity; and (4) to conduct drought risk regionalization.

Data Source and Collection
In this study, the utilized monthly gridded precipitation dataset (version 2.0) has a spatial resolution of 0.5° × 0.5° and encompasses a time span ranging from 1961 to 2020. The dataset was acquired from the China Meteorological Data Service Centre [35]. Thin-plate smooth spline interpolation was employed to construct the gridded dataset based on the precipitation observations from 2472 meteorological stations across China [36]. Excluding Taiwan and South China Sea Islands, the dataset comprises 4189 grids and was verified by average root-mean-square and cross-validated by the Tropical Rainfall Measurement Mission product [36,37]. As a result, the dataset has an average error of 0.49 mm/month, and the average correlation coefficient between the grids and the observations is 0.93 (arriving at a significance level of α = 0.01) [38]. The dataset has also been extensively applied for research on climate change in China [1,30,39]. The multi-year average precipitation was calculated using this dataset and is shown in Figure 1. The spatial pattern of annual precipitation indicates a gradual decrease from the southeastern coast to the northwest inland. In addition, the seasonal distribution of precipitation is also uneven, which mainly falls in May-October in the form of high-intensity rainforms [40,41]. The uneven spatial distribution of precipitation and the considerable disparities in the time domain have led to frequent droughts and floods in China, which have greatly threatened agricultural production, as well as people's lives and properties.

Standardized Precipitation Index
The SPI index proposed by Mckee et al. [6] can characterize droughts on various time scales, such as 1-, 3-, 6-, and 12-month periods. The detailed computation formulas can be found in the previous studies [10,42,43], and the key steps are as follows [44,45]: (1) Calculating the cumulative precipitation on a specific time scale; (2) Matching the probability distribution of accumulative precipitation using a gamma distribution function; (3) Determining the non-exceedance probabilities for cumulative precipitation; (4) Converting the probabilities to a standard normally distributed variable. SPI is classified into five groups based on the "Meteorological drought grade (GB/T20481-2017)" [46], as shown in Table 1.

Standardized Precipitation Index
The SPI index proposed by Mckee et al. [6] can characterize droughts on various time scales, such as 1-, 3-, 6-, and 12-month periods. The detailed computation formulas can be found in the previous studies [10,42,43], and the key steps are as follows [44,45]: (1) Calculating the cumulative precipitation on a specific time scale; (2) Matching the probability distribution of accumulative precipitation using a gamma distribution function; (3) Determining the non-exceedance probabilities for cumulative precipitation; (4) Converting the probabilities to a standard normally distributed variable. SPI is classified into five groups based on the "Meteorological drought grade (GB/T20481-2017)" [46], as shown in Table 1.  [27]. In this study, we employed SPI on the time scales of 3-month (SPI-3), 6-month (SPI-6), and 12-month (SPI-12) to explore the short-, medium-, and long-term drought characteristics in mainland China (excluding Taiwan and South China Sea Islands). The SPI value of a specific accumulation period was assigned to the final month of that period, so the SPI-6 value for June 2020 referred to the SPI for the accumulative precipitation from January to June 2020.

Run Theory
The run theory introduced by Yevjevich [15] makes it easy to identify drought features including start and end times, as well as the duration and severity of drought. For this study, a drought threshold of −0.5 was set based on Table 1. There were two aspects to be dealt with in the process of drought identification: (1) Minor drought events. When SPI was less than −0.5, we preliminarily judged that the month was one drought-month. It can be seen in Figure 2 that there were six distinct droughts labelled E 1 , E 2 , E 3 , E 4 , E 5 , and E 6 . A single-month drought event was considered to have occurred if its associated SPI value was less than −1.5 (such as E 6 ). Otherwise, it was considered a minor drought event and should be ignored (such as E 1 ); (2) A confluence of droughts. In the case of two consecutive droughts separated by a month, if the SPI value was between −0.5 and 0.5 in the interval month, the two adjacent droughts were considered subordinate droughts. Then, the two droughts should be combined into one drought event (such as E 4 and E 5 ). The total drought duration equaled E 4 (d) + E 5 (d) + 1, while the drought severity equaled E 4 (s) + E 5 (s). Otherwise, it was two independent drought events. As a result, four drought events were finally identified, namely E 2 , E 3 , E 4 + E 5 and E 6 . Wherein drought duration (d) refers to the number of months elapsed between its onset and its cessation, while drought interval (l) is the amount of time between its onset and that of the next adjacent drought. The drought severity (s) is the cumulative difference that the SPI falls below −0.5 for each drought episode. To facilitate computation, we doubled the drought severity by −1 to get a positive value [47,48], namely:

Copula Function
Sklar [50] pointed out that copula functions could show how drought variables are related even if they don't have the same marginal distribution. They can depict the effects of drought on a specific region with more precision than conventional statistical methods.
The Archimedean copula function family, which includes Clayton, Gumbel, and Frank copulas, has been widely used in hydrology and meteorology. Selecting the appropriate copula function directly impacts the analysis and computation outcomes [51]. Thus, it is crucial to identify the copula function that can effectively represent the variables in the analysis. Here, the best-fit copula was selected based on the Akaike Information Criterion (AIC) [52]. The selection process followed a criterion where a lower value indicated better performance. The formula for AIC was provided in Equation (5). This paper discussed the drought duration, severity, and maximum (the peak value of drought severity) [49] as defining features of drought episodes. To calculate these attributes for each grid, the following formulas were used: where n represents the total number of droughts.

Copula Function
Sklar [50] pointed out that copula functions could show how drought variables are related even if they don't have the same marginal distribution. They can depict the effects of drought on a specific region with more precision than conventional statistical methods.
The Archimedean copula function family, which includes Clayton, Gumbel, and Frank copulas, has been widely used in hydrology and meteorology. Selecting the appropriate copula function directly impacts the analysis and computation outcomes [51]. Thus, it is crucial to identify the copula function that can effectively represent the variables in the analysis. Here, the best-fit copula was selected based on the Akaike Information Criterion (AIC) [52]. The selection process followed a criterion where a lower value indicated better performance. The formula for AIC was provided in Equation (5). After selecting the appropriate copula function, its parameters were evaluated by the marginal inference function, which consisted of two separate estimation steps: first estimating each univariate marginal distribution function based on the maximum likelihood method; and then estimating the copula dependence parameter [53]. It should be also noted that the presence of a correlation between variables is essential in the creation of a bivariate distribution using a copula function. Without this correlation, the use of a copula function is not possible [54]. To determine the correlation between drought severity and duration in this study, Kendall's correlation coefficient [55] was employed. Prior research [24,56] revealed that the cumulative distribution functions of drought duration (F D (d)) and severity (F S (s)) follow the exponential and gamma distributions, as defined in Equations (6) and (7), respectively.
where the variable L max represents the maximum log-likelihood that can be attained by the model for a given dataset; k refers to the number of parameters that can be freely estimated for the model. λ denotes the parameter to describe exponential distribution, α and β determine the shape and scale of gamma distribution. The copula-based joint cumulative distribution function of drought duration and severity may be represented as follows: where θ is a measure for quantifying the strength of the link between F D (d) and F S (s).

Return Period and Joint Probability
The return period of a drought refers to an estimated average time between two drought events with a specified value or higher [20]. Given that droughts may persist longer than one year, Shiau and Shen [57] developed the following formula to calculate the return period of drought duration (R D ) and drought severity (R S ): where E(l) denotes the expected value of drought interval time as described in Section 2.2.2.
Drought risk assessment and drought management rely heavily on information on the joint probability of drought features. In this study, we assessed the joint probability that both drought duration and severity concurrently surpass the specified threshold (i.e., D ≥ d and S ≥ s). The formula for the probability was given by Shiau [23]:

Drought Risk Regionalization
Drought risk regionalization is a vital task in the field of drought risk management, as it is the process of identifying areas that are likely to be affected by drought hazards. This is achieved through the aggregation of regions that have similar climatic characteristics and spatial continuity. The process of drought risk regionalization is complex, as it involves the use of sophisticated statistical and geographical techniques, such as clustering and spatial analysis. By implementing these methods, researchers can effectively identify areas that are susceptible to drought and subsequently develop strategies to mitigate its impact. In this study, we employed the HiClimR package [58] in R software to accomplish the drought risk regionalization in mainland China. The scheme for drought risk regionalization can be summarized as follows: (1) Extract the longitude and latitude of each grid point and generate the corresponding raster layer. (2) Select the dataset of input data for the implementation of hierarchical clustering.
The dataset includes latitude, longitude, and grid-based joint probability of drought duration and severity calculated by Equation (11). (3) Preprocess the input data by removing the effect of magnitude through standard deviation normalization. (4) Identify adjacent grids by using the latitude and longitude variables. (5) Calculate the distance between the grids using the Pearson correlation coefficient and the distance between classes using the sum of squares method. (6) Determine the adjacency and homogeneity of the regions by minimizing the interregional correlation coefficient and maximizing the intra-regional correlation coefficient. (7) Obtain the final group numbers based on the silhouette width criterion. (8) Group the adjacent grids based on their similarity in terms of joint probability of drought duration and severity. (9) Evaluate the results of the drought risk regionalization and identify the areas that are susceptible to drought hazards.

Regional Average Drought Attributes
Drought attributes, such as average duration and severity, maximum severity, and the number of drought episodes, were depicted spatially in Figure 3 at 3-, 6-, and 12-month time scales. On the 3-month time scale, most of Xinjiang, northern and southern Tibet, western Qinghai, and eastern Heilongjiang experienced relatively longer drought durations, while the average drought durations were relatively short in eastern Inner Mongolia, southwestern Heilongjiang, eastern Jilin, and most of Shaanxi province (Figure 3a). The distribution patterns of average drought severity showed that western and northern Xinjiang, northern Tibet, southwestern and eastern Qinghai, central and eastern Yunnan, southern Yangtze River Basin, central Anhui, and eastern Heilongjiang experienced relatively more severe droughts. In contrast, the values of drought severity were relatively low in southeastern Xinjiang, western Tibet, western Inner Mongolia, and eastern Jilin (Figure 3b). Northern Xinjiang, northeastern Tibet, southern Qinghai, and central-northern Yunnan have experienced more intense droughts (Figure 3c). Compared with southern China, more droughts occurred in northern China (Figure 3d). and northern Xinjiang, northern Tibet, southwestern and eastern Qinghai, central and eastern Yunnan, southern Yangtze River Basin, central Anhui, and eastern Heilongjiang experienced relatively more severe droughts. In contrast, the values of drought severity were relatively low in southeastern Xinjiang, western Tibet, western Inner Mongolia, and eastern Jilin (Figure 3b). Northern Xinjiang, northeastern Tibet, southern Qinghai, and central-northern Yunnan have experienced more intense droughts (Figure 3c). Compared with southern China, more droughts occurred in northern China (Figure 3d). The average duration and severity followed similar spatial patterns on the 6-month time scale as they did on the 3-month time scale (Figure 3e,f). However, the most severe droughts were also observed along the southeastern coast (Figure 3g), and the drought frequency pattern was more pronounced, with more droughts occurring north of the Yangtze River than south of it (Figure 3h).
On the 12-month time scale, the average drought durations in Northeast China, Northwest China, North China, and Qinghai-Tibet Plateau were longer than that in the southeast coastal provinces (Figure 3i). Drought severity followed a similar pattern as drought duration (Figure 3j), which meant that longer drought duration tended to result in more severe drought. The most severe droughts occurred in southwestern China and at the junction of Inner Mongolia, Jilin, and Liaoning provinces (Figure 3k). The pattern of drought frequency was quite different from that on the 3-and 6-month time scales. Although the drought durations in central and eastern China were shorter, the drought frequencies were relatively high (Figure 3l). Similar to the drought severity, the areas with higher numbers of drought episodes decreased as time scales increased. Compared to the spatial patterns of short-and medium-term drought attributes, long-term drought attributes had quite distinct regional distributions. So, when examining drought-related The average duration and severity followed similar spatial patterns on the 6-month time scale as they did on the 3-month time scale (Figure 3e,f). However, the most severe droughts were also observed along the southeastern coast (Figure 3g), and the drought frequency pattern was more pronounced, with more droughts occurring north of the Yangtze River than south of it (Figure 3h).
On the 12-month time scale, the average drought durations in Northeast China, Northwest China, North China, and Qinghai-Tibet Plateau were longer than that in the southeast coastal provinces (Figure 3i). Drought severity followed a similar pattern as drought duration (Figure 3j), which meant that longer drought duration tended to result in more severe drought. The most severe droughts occurred in southwestern China and at the junction of Inner Mongolia, Jilin, and Liaoning provinces (Figure 3k). The pattern of drought frequency was quite different from that on the 3-and 6-month time scales. Although the drought durations in central and eastern China were shorter, the drought frequencies were relatively high (Figure 3l). Similar to the drought severity, the areas with higher numbers of drought episodes decreased as time scales increased. Compared to the spatial patterns of short-and medium-term drought attributes, long-term drought attributes had quite distinct regional distributions. So, when examining drought-related issues, an appropriate time scale should be selected depending on the purpose of the study. Figure 4 illustrated the drought severity for 3-, 6-, and 12-month time scales at different return periods, including 5-, 10-, 20-, 50-, and 100-year. The study findings suggested that drought severity increased with longer drought return periods across various time scales [27]. On the 3-month time scale, severe droughts were observed in most regions of China except the Changbai Mountains, Greater Khingan Mountains, Pearl River Basin, western Tibet, and eastern Xinjiang province. Moreover, similar spatial patterns of drought severity were observed at 3-and 6-month time scales. It should be noted that severe droughts have also been found along the southeastern coast of China. On the 12-month time scale, the Qinghai-Tibet Plateau was identified to be more susceptible to droughts. ous time scales [27]. On the 3-month time scale, severe droughts were observed in most regions of China except the Changbai Mountains, Greater Khingan Mountains, Pearl River Basin, western Tibet, and eastern Xinjiang province. Moreover, similar spatial patterns of drought severity were observed at 3-and 6-month time scales. It should be noted that severe droughts have also been found along the southeastern coast of China. On the 12-month time scale, the Qinghai-Tibet Plateau was identified to be more susceptible to droughts.  Drought duration for 3-, 6-, and 12-month time scales with return periods of 5-, 10-, 20-, 50-, and 100-year were depicted in Figure 5. In comparison to Figure 4, we discovered that areas where droughts lasted longer were likewise associated with more severe droughts. This means that longer drought durations tended to indicate more severe droughts and vice versa. From Figures 4 and 5, it can be observed that more severe and longer drought events were detected in western and northern Xinjiang, southwestern China, and the southern regions of the mid-lower reaches of the Yangtze River. In contrast, less intense and shorter drought events were monitored mainly in the Greater Khingan Mountains, Changbai Mountains, and Pearl River Basin. To confirm the positive relationship between drought duration and severity, we analyzed the spatial distribution of Kendall correlation coefficients for these two variables at different time scales. Specifically, we plotted the correlation coefficients for the 3-, 6-, and 12-month time scales in Figure 6. The findings indicate that there is indeed a positive correlation between drought duration and severity, and this relationship reaches a significance level of α = 0.05. Moreover, the Kendall correlation coefficient shows an increase with longer time scales, suggesting that the longer the drought lasts, the more severe it becomes. Nonetheless, univariate analysis is insufficient to completely characterize drought episodes due to the correlation between drought attributes.

Spatial Patterns of Drought Severity and Duration with Various Return Periods
Khingan Mountains, Changbai Mountains, and Pearl River Basin. To confirm the positive relationship between drought duration and severity, we analyzed the spatial distribution of Kendall correlation coefficients for these two variables at different time scales. Specifically, we plotted the correlation coefficients for the 3-, 6-, and 12-month time scales in Figure 6. The findings indicate that there is indeed a positive correlation between drought duration and severity, and this relationship reaches a significance level of α = 0.05. Moreover, the Kendall correlation coefficient shows an increase with longer time scales, suggesting that the longer the drought lasts, the more severe it becomes. Nonetheless, univariate analysis is insufficient to completely characterize drought episodes due to the correlation between drought attributes.

Joint Probability of Drought Duration and Severity
To gain a better understanding of the relationship between drought severity and duration, it is essential to evaluate the drought risks in mainland China using bivariate analysis. Firstly, the best-fit copula function was determined grid-by-grid based on the AIC criteria, and the Frank copula was found to have the best performance for each grid at various time scales. Next, we employed the joint probability, as described in Equation (11), with the thresholds s and d derived from the preceding univariate study of drought severity and duration for 5-, 10-, 20-, 50-, and 100-year return periods, respectively. The results displayed in Figure 7 showed a clear spatial distribution of the joint probability P

Joint Probability of Drought Duration and Severity
To gain a better understanding of the relationship between drought severity and duration, it is essential to evaluate the drought risks in mainland China using bivariate analysis. Firstly, the best-fit copula function was determined grid-by-grid based on the AIC criteria, and the Frank copula was found to have the best performance for each grid at various time scales. Next, we employed the joint probability, as described in Equation (11), with the thresholds s and d derived from the preceding univariate study of drought severity and duration for 5-, 10-, 20-, 50-, and 100-year return periods, respectively. The results displayed in Figure 7 showed a clear spatial distribution of the joint probability P at various time scales and return periods. On the 3-and 6-month time scales, relatively large P values were mainly observed in northern Xinjiang, western Qinghai, southern Tibet, central Yunnan, and eastern Heilongjiang, indicating a higher likelihood of experiencing droughts in these regions. Conversely, the P values in the southeastern coastal areas of China, Changbai Mountains, and Greater Khingan Mountains were mostly smaller, implying a lower drought risk in these regions. On the 12-month time scale, the drought risk on the Tibetan Plateau was the highest, which was consistent with the results obtained from the univariate analysis.

Drought Risk Regionalization
In order to identify drought-prone regions in mainland China, we mapped the drought risk regionalization at various return periods, as shown in Figure 8. With the aid of the silhouette width criterion, a preliminary estimate of the optimal number of clusters was made. According to the findings, the optimal number of drought risk zones ranged from four to six categories over various return periods. For convenience of analysis, we divided mainland China into six sub-regions considering the topography, drought risk, and the continuity and independence of regions. The map visualization

Drought Risk Regionalization
In order to identify drought-prone regions in mainland China, we mapped the drought risk regionalization at various return periods, as shown in Figure 8. With the aid of the silhouette width criterion, a preliminary estimate of the optimal number of clusters was made. According to the findings, the optimal number of drought risk zones ranged from four to six categories over various return periods. For convenience of analysis, we divided mainland China into six sub-regions considering the topography, drought risk, and the continuity and independence of regions. The map visualization was accomplished using Arcmap 10.4. The mapping process involved the use of the focal statistics function for raster smoothing, Chinese border extraction, raster to polygon conversion, and modifying isolated classification grid points by changing their attributes to that of the adjacent major class. It was clear from Figure 8 that the spatial pattern of drought risk zones exhibited a high degree of similarity on 3-and 6-month time scales with return periods ranging from 5 to 100 years. The six sub-regions were identified as Xinjiang (I), Qinghai-Tibet Plateau (II), central and western Inner Mongolia (III), South China (IV), North China (V), and Northeast China (VI). However, for the 12-month time scale, the zoning results under different return periods mostly separated Southwest China, suggesting that this region was particularly susceptible to long-term droughts. In addition, the zone (central and western Inner Mongolia) was absorbed by neighboring zones. China, suggesting that this region was particularly susceptible to long-term droughts. In addition, the zone (central and western Inner Mongolia) was absorbed by neighboring zones.

Drought Characteristics and Risk
The study observed that regions with less precipitation and high rainfall variabilities, such as Northwest China and the Qinghai-Tibet alpine region, experienced longer

Drought Characteristics and Risk
The study observed that regions with less precipitation and high rainfall variabilities, such as Northwest China and the Qinghai-Tibet alpine region, experienced longer and more severe droughts. Several studies have highlighted that the water resource supply is decreasing, which makes these regions more susceptible to drought [59][60][61][62]. However, other studies have shown an increase in precipitation in the Qinghai-Tibet Plateau and Northwest China, indicating a possible improvement in drought conditions in these regions in the future [63,64]. Additionally, severe drought episodes were observed in the middle and lower reaches of the Yangtze River and the Yunan-Guizhou Plateau, which might have been due to a significant reduction in precipitation [12,16,17,65]. Another study conducted by He et al. [66] identified a higher drought risk in the northern Xinjiang, Loess Plateau, Northeast Plain, southern areas of the Yangtze Plain, and Yunnan-Guizhou Plateau. Our research, as depicted in Figure 7, also reveals similar findings, although we identified high drought risk in eastern Heilongjiang instead of most of the Northeast Plain. The reason for the discrepancy might be that our primary focus was meteorological drought rather than agricultural drought, which was the main concern in their study.

Regionalization
In China, previous studies on regionalization focused on employing various indices to determine climate zones. For example, Yang and Li [67] applied the rotated empirical orthogonal function (REOF) method to divide China into 11 subregions based on the precipitation anomaly percentage. However, due to the empirical nature of the REOF method, there is a certain subjectivity and fuzziness in determining the boundaries of the subregions. Li et al. [68] identified eight homogeneous regions in China by applying Ward's and k-means clustering techniques based on SPI-3, but still encountered subjectivity in the determination of the subregion boundaries attributed to the utilization of station-based data. Drought risk regionalization has only been explored in the research conducted by Wu et al. [69], who took into account the impact of the long-term evolution of drought variables on the results of drought risk regionalization, but only considered the drought severity. In contrast, our study examined the joint probabilities of drought duration and severity at different drought levels, providing more detailed and helpful information for drought risk regionalization. As a result, we divided mainland China into six subregions based on the joint probability. However, it should be noted that the spatial patterns of long-term drought features differed from those on short-and medium-term time scales (seen from Figures 3-6) because the time scale played an essential role in the spatial heterogeneity of drought behaviors. Therefore, when discussing drought-related scientific issues, the time scale must be considered according to the research purpose, as pointed out by McKee et al. [6], due to the complexity of drought phenomena.

Limitations and Future Work
This study focused only on a bivariate analysis of drought, specifically duration and severity. However, given that droughts are complex and multifaceted phenomena, relying solely on a bivariate analysis may lead to less accurate estimates of drought probabilities compared to more complex three-or four-variate analyses [70]. Hence, researchers have developed multivariate extensions of bivariate copula that can provide a more comprehensive assessment of droughts [70][71][72][73]. Nonetheless, the parameter estimate and calculation [70][71][72]74] will get more challenging as the number of drought-related variables increases [75]. Additionally, under the backdrop of climate change, the uncertainties surrounding the incidence and progression of drought present an imminent threat to both China's ecological and food security. Thus, predicting drought in the future becomes crucial. To this end, our forthcoming research will center on employing SPI and highdimensional copula functions to assess the spatiotemporal evolution of drought, as well as the corresponding drought risk under RCP4.5 and RCP8.5 scenarios. Furthermore, we will conduct a comprehensive analysis of SPI under the two scenarios to explicate the potential uncertainties of drought forecasting.

Conclusions
Since drought is a serious natural disaster, understanding the characteristics of drought is an essential step in establishing effective measures to alleviate drought-related problems. In this study, China's monthly gridded precipitation dataset from 1961 to 2020, with a spatial resolution of 0.5 • × 0.5 • , was used to calculate the SPI values at various time scales and identify drought events in combination with run theory. The features of drought events were described by drought duration and severity, which were fitted by exponential and gamma distributions, respectively. In addition, the Frank copula function was employed to establish the dependence between drought duration and severity. According to the formulas of return period and joint probability, the thresholds of drought duration and severity corresponding to the return periods of 5-, 10-, 20-, 50-, and 100-year under various time scales were calculated, as well as the values of joint probability. As a result, the key findings were as follows: (1) The spatial patterns of average drought duration and severity indicated that western and northern Xinjiang, western and southern Tibet, western Qinghai, central and northern Yunnan, and eastern Inner Mongolia experienced severe droughts. The spatial patterns of the average drought duration and severity were similar, suggesting that longer drought duration tended to correspond to higher drought severity and vice versa. Although droughts occurred relatively frequently in most of the three northeastern provinces, in the middle and lower reaches of the Yellow River and Yangtze River the severity was relatively mild. Notably, the spatial patterns of the average drought variables on the 12-month time scale differed remarkably from those on the 3-and 6-month time scales.
(2) In the univariate analysis, the spatial coverages of drought duration and severity at various time scales and return periods were also consistent, which further indicated that drought duration and severity were positively correlated. Drought events with longer duration and higher severity mainly occurred in western and northern Xinjiang, southwestern China, and the southern part of the middle and lower reaches of the Yangtze River. Drought events with shorter duration and less intensity were most common in the Greater Khingan Mountains, Changbai Mountains, and Pearl River Basin.
(3) The results obtained from bivariate frequency analysis, using the Frank copula function, revealed that the regions with higher P values on the 3-and 6-month time scales were northern Xinjiang, western Qinghai, southwestern China, and the middle and lower reaches of the Yangtze River. Conversely, the regions with lower P values were the southeastern coastal regions of China, Changbai Mountains, and Greater Khingan Mountains. On the other hand, for the 12-month time scale, the drought risk was found to be the highest on the Tibetan Plateau.
(4) Mainland China could be divided into six sub-regions according to joint probabilities of drought duration and severity. However, the spatial pattern of drought risk zones exhibited a high degree of similarity on 3-and 6-month time scales but differed from the results on the 12-month time scale.
Overall, regardless of average drought characteristics, joint drought probability, or drought risk regionalization, time scale played an essential role in the spatial heterogeneity of their behaviors, reflecting the complexity of drought phenomena. Additionally, our method, combining SPI and bivariate copulas for drought risk regionalization, is generalizable to other regions and can be applied to other climatic variables, although some modifications may be necessary to fit specific conditions.