Drought risk analysis based on multivariate copula function in Henan Province, China

Abstract

Global droughts have become more frequent in recent years, posing a serious threat to human production and life. In order to overcome the limitations of traditional studies in quantifying drought risk, this study uses runs theory to extract drought duration, severity and kurtosis as drought characteristic variables from standardized precipitation index (SPI) in Henan Province. Drought risk analysis models based on multivariate copula function are constructed to reveal the response relationships between drought characteristic variables. The results show that the smaller the SPI scale, the higher the sensitivity to identify drought processes. In addition, the multivariate copula function shows good fitting performance for the optimal joint distribution of drought characteristic variables, with R² values exceeding 0.9. Drought joint recurrence period is positively correlated with drought characteristic variables, and the recurrence period within the duration of 4–6 months occur frequently, indicating a higher probability of experiencing short-term droughts and cross-seasonal droughts. When the drought duration, severity and kurtosis are greater than 2.5, 2200 and 1300, respectively, the drought joint recurrence period reaches 20 months. The research results have provided new methods for drought risk analysis and data support for formulating drought mitigation strategies in Henan Province.

Keywords:

• Drought
• Standardized Precipitation Index (SPI)
• Runs theory
• multivariate copula function
• risk analysis

1. Introduction

Water resources are an invaluable and irreplaceable asset that supports sustainable development across the globe (Zhang et al. 2019). However, the water cycle on Earth is currently experiencing significant changes as a result of human activities and global climate change (Kogan et al. 2020). Almost half of the world’s countries are currently facing the threat of drought, posing a common challenge for humanity (Maeng et al. 2017; Bai et al. 2018). The 15th UN Convention to Combat Desertification’s "Drought Report" highlights that drought, which accounts for 15% of natural disasters, results in the highest number of deaths. Between 1970 and 2019, it caused approximately 650,000 deaths worldwide. Since 2000, the frequency and duration of global droughts have increased by 29%, affecting roughly 1.4 billion people. Between 1998 and 2017, the global economic losses attributed to drought reached a staggering $124 billion. The report suggests that by 2030, 700 million individuals will face the risk of displacement due to drought. The Intergovernmental Panel on Climate Change (IPCC) state (Lynn and Peeva 2021): most parts of the world will become more arid in future, and the degree of aridification will be remotely higher than the changes of the past millennium, highlighting the urgency of in-depth study of drought. Despite its severe implications, drought has become a major impediment to sustainable human development, owing to its intricate natural environment and multifarious impacts from various complex factors (Han et al. 2019). Drought risk analysis serves as a crucial metric to quantify the severity of drought to some extent and constitutes a crucial cornerstone for identifying and issuing early warnings.

Over the past decade, there has been a notable surge in research on drought risk, as evidenced by the explosive increase in the number of search results for "drought disaster risk" from 505,000 in 2016 to 6.13 million in 2020 (Liu et al. 2021). Drought and drought risk have become important research areas in the fields of natural disasters and sustainable development (Duan et al. 2022). For instance, Zhao Jing (Zhao et al. 2012) analyzed drought disaster risk and created a drought risk zoning graph using GIS, examining drought disasters in different regions and time periods. Meanwhile, Xu (Xu et al. 2015) developed a regional drought frequency analysis model based on copula functions by considering the spatial and temporal variability of drought events, and assessed and discussed the robustness of the drought frequency analysis as well as the extent of drought. In addition, Zhang (Zhang et al. 2015) constructed the joint distribution of drought characteristics by using the Archimedean Copula function. Based on the error analysis and tail dependence coefficient, the Gumbel-Hougaard copula is selected to analyse spatial distributions of the drought joint recurrence period. Wu (Wu et al. 2019) used the copula function to construct the joint distribution of drought risk management indicators in grassland ecosystems based on SPEI data, and the results indicated that the Frank Copula function could well simulate the joint dependence structure among the indicators. Zhong (Zhong et al. 2018) used on Bayesian maximum entropy method for optimal estimation of annual precipitation and calculated the annual standardized precipitation index. The hazard of drought was also assessed based on information diffusion and excess probability theory. Li (Li et al. 2020) proposed an overall framework and key technologies for quantitative evaluation of drought risk based on a scientific definition of drought risk. Richter (Richter and Semenov 2005) employed the Hadley Climate Center’s climate model and a stochastic weather generator to simulate the climate scenario for England in 2050 and evaluated the impact of future drought risk on wheat yield. Eleni (Kampragou et al. 2011) suggested that governments at all levels in the EU should shift from traditional drought crisis management to drought risk management at regional, watershed, and national levels. Kirono (Kirono et al. 2011) used global climate models (GCMs) to simulate future climate scenarios under high greenhouse gas emissions and analyzed drought characteristics in Australia from 2030 to 2050. Finally, Wilhite (Wilhite et al. 2000; Wilhite 2016; de Araujo et al. 2019) pointed out that current risk management strategies are inadequate to deal with the increasingly severe drought disasters, and argued for the need to manage drought as a risk in the context of changing climates. The copula function is flexible and versatile. On the one hand, it is possible to construct multidimensional joint distributions by considering the marginal distributions and correlation structures of different drought variables without loss of data information. On the other hand, the copula function model can account for the complex relationships between multiple variables, overcoming the limitations of traditional methods and providing a powerful framework for comprehensive drought risk analysis.

In summary, there are several shortcomings in the current analysis of drought risk. Firstly, the identification and analysis of drought risk often focuses on a single indicator, which is limited in scope. Secondly, there is a certain degree of missing data information in the processing and transformation of drought index data, which can lead to inaccurate analysis results. In addition, there is limited research on drought risk analysis using the construction of joint distributions of three or more dimensions. Therefore, it is crucial to establish a sound mechanism for drought risk analysis, reveal the response relationship between different drought characteristic variables, and accurately evaluate the drought situation in Henan Province.

This paper aims to construct a joint distribution among drought characteristic variables using a multivariate copula function to address the limitations of drought risk analysis in Henan Province. Firstly, the calculation of SPI based on precipitation data and the trend analysis of SPI values at different scales are carried out. Secondly, three drought characteristic variables, drought duration, drought severity, and drought kurtosis, were extracted from SPI1 based on runs theory. Then, the joint distribution among the drought characteristic variables is constructed based on the multivariate copula function. Finally, by calculating the conditional probability of the joint distribution among drought characteristic variables, the joint recurrence period is determined and the analysis of drought risk in Henan Province is realized. This study overcomes the information loss between data and achieves a comprehensive analysis of drought risk in Henan Province from a multivariate perspective.

2. Study area

Henan Province, known as Zhongyuan in ancient times, is located in the central and eastern regions of China, within the middle and lower reaches of the Yellow River. Henan is one of the major provinces in China with a large population and extensive agricultural production, covering an area of 165,500 square kilometers. Figure 1 shows the research area of Henan Province, where the terrain is mainly characterized by plains, complemented by hills and mountains. Specifically, hills account for 17.7% of the basin area, mountains account for 26.6%, and plains account for 55.7% of the basin area (Luo and Li 2023).

Figure 1. Research area.

Henan Province has a diverse geological structure, which results in uneven precipitation distribution throughout the year and pronounced seasonal changes. Specifically, the winter season sees less precipitation, while the summer season experiences ample precipitation, and there are significant inter-annual variations (Wang et al. 2018). On average, the minimum precipitation is 500 mm, while the annual maximum precipitation is 900 mm. Additionally, there are notable differences in precipitation patterns across different regions, with a decreasing trend from south to north. Over the past 60 years (1951–2010), the average precipitation is 625.6 mm, and the precipitation trend rate is −4.542 mm/10a. The intergenerational characteristics of annual precipitation exhibit a downward trend from year to year, and the frequency of droughts has been as high as 58.2%. Given the constantly changing drought situation, the risk of drought is increasing every year.

2.1. Data used

Hydro-meteorological data forms the fundamental basis for analysis and calculation. In this study, the hydro-meteorological data for 18 cities in Henan Province from 1961 to 2010, including monthly precipitation, temperature, wind speed, relative humidity, and other observational data, were obtained from the China Meteorological Data Service Center (http://data.cma.cn).

In addition, the calculation and transformation of the multi-scale SPI based on the precipitation data were carried out, taking into account the current climate conditions. The analysis of drought risk in Henan Province was carried out at different time scales to improve the limitations of using only single precipitation data for risk analysis results.

3. Methodology

3.1. Calculation of the SPI

Proposed by McKee in 1993, the SPI (Belayneh et al. 2014; Li et al. 2020) describes the severity of drought with the probability of precipitation. This index applies to droughts on multiple time scales. Suppose precipitation x obeys gamma distribution. Then, the probability density of x, i.e. the SPI, can be described as: (1) $$\begin{array}{c}g(x)=\frac{1}{\beta \Gamma (\alpha )}{x}^{\alpha -1}{e}^{-\frac{x}{\beta }}(x>0)\\ \Gamma (\alpha )=\underset{0}{\overset{\infty }{\int }}{y}^{\alpha -1}{e}^{-y}dy\end{array}$$(1) where $\alpha$ is the shape factor; $\beta$ is the scale factor; $\Gamma (\alpha )$ is the gamma function. These parameters can be computed through maximum likelihood estimation (MLE): (2) $$\alpha =\frac{1+\sqrt{1+4A/3}}{4A}\beta =\frac{\overline{x}}{\alpha }A=\hspace{0.17em}\ln \hspace{0.17em}(\overline{x})-\frac{\sum _{i=1}^n\hspace{0.17em}\ln \hspace{0.17em}(x_i)}{m}$$(2)

The SPI-based drought classification standard (Table 1) was established according to the Classification of Meteorological Drought (GB/T20481-2006) issued by the National Meteorological Center, CMA.

Table 1. SPI-base drought classification.

Grade	Type	SPI value
1	No drought	−0.5 < SPI
2	Light drought	−1.0 < SPI≤−0.5
3	Moderate drought	−1.5 < SPI≤−1.0
4	Severe drought	−2.0 < SPI≤−1.5
5	Extreme drought	SPI≤−2.0

3.2. Runs theory

Based on the SPI 1, the runs theory (Zhang et al. 2015; Long et al. 2020) is used to extract three drought characteristic variables from them: drought duration, drought severity, and drought kurtosis. Three cut-off levels, namely $R_0,R_1,R_2,$ are established. Drought occurs during a period when the standardized precipitation index sequence $R_i(i=1,2\dots \dots N)$ for one or more periods is less than R₁, indicating a negative run. The specific identification procedure is as follows: Drought duration (D, monthly) refers to the length of the period of two or more drought occurrences, which constitutes a negative run. Drought severity (S, 10,000 m³) is the area of the negative run. Drought kurtosis (K, 10,000 m³) refers to the extreme value of the negative run. The runs theory diagram is shown in Figure 2.

Figure 2. Schematic diagram of runs theory.

3.3. Copula function

The copula function has a flexible and diverse form, and is easy to solve. In multivariate analysis, determining the marginal distribution of a single variable given the joint distribution of multiple variables is relatively straightforward. However, determining the joint distribution based on the marginal distributions of each variable can be challenging. Copula function theory has been proposed to address this problem. As a linking function (Li et al. 2021), it can transform the distribution between the two without compromising the data information of the random variables. According to the Sklar theorem (Afshar et al. 2016; Li et al. 2020), any multivariate joint distribution can be connected to single-variable marginal distribution functions and a copula function that describes the dependence structure between variables. The specific formula for this definition is as follows:

$F(x_1,x_2,\dots x_n)$ is the joint distribution of the marginal distribution $F_1(x_1)=P(X_1\le x_1),\dots F_n(x_n)=P(X_n\le x_n).$ Let $u_i=F_i(x_i)(i=1,2,\dots ,n),$ then there exists an n-dimensional Copula function C. For any $x\in R_n,$ we have $F(x_1,x_2,\dots x_n)=C(u_1,u_2,\dots u_n)$

3.3.1. Two-dimensional copula functions

The Copula function is a uniform joint distribution function with marginal distributions restricted to the interval [0,1]. Let F(x) be the marginal distribution function of a one-dimensional random variable X, Let G(y) denote the marginal distribution function of a one-dimensional random variable Y, and H(x,y) denote their joint distribution function. The function C serves as a link function between the marginal and joint distribution functions. For fitting purposes in this study, an Archimedean Copula function is employed (see Table 2).

Table 2. Two-dimensional Archimedean copula function.

Copula function	$C(u_1,u_2)$	$\theta$ values
Clayton	${(u_1{}^{-\theta }+u_2{}^{-\theta }-1)}^{\frac{-1}{\theta }}$	$(0,x)$
Gumbel	$\hspace{0.17em}\exp \hspace{0.17em}\{-{[{(-\hspace{0.17em}\ln \hspace{0.17em}{u}_1)}^{\theta }+{(-\hspace{0.17em}\ln \hspace{0.17em}{u}_2)}^{\theta }]}^{\frac{1}{\theta }}\}$	$[1,\infty )$
Frank	$-\frac{1}{\theta }\ln [1+\frac{(e^{-\theta u_1}-1)(e^{-\theta u_2}-1)}{(e^{-\theta }-1)}]$	$(-\infty ,\infty )\backslash \left\{0\right\}$

3.3.2. Three-dimensional copula functions

The method for constructing a three-dimensional Copula function is similar to that of a two-dimensional Copula function. However, due to the tendency for drought data to exhibit skewed distribution, the drought characteristic variable data are fitted in this paper based on the Clayton Copula and Gumbel-Hougaard Copula functions (see Table 3).

Table 3. Three-dimensional Archimedean copula function.

Copula function	$C(u_1,u_2,u_3)$	$\theta$ values
Clayton Copula	${(u_1{}^{-\theta }+u_2{}^{-\theta }+u_3{}^{-\theta }-2)}^{\frac{-1}{\theta }}$	$(0,x)$
Gumbel-Hougaard Copula	$\hspace{0.17em}\exp \hspace{0.17em}\{-{[{(-\hspace{0.17em}\ln \hspace{0.17em}{u}_1)}^{\theta }+{(-\hspace{0.17em}\ln \hspace{0.17em}{u}_2)}^{\theta }+{(-\hspace{0.17em}\ln \hspace{0.17em}{u}_3)}^{\theta }]}^{\frac{1}{\theta }}\}$	$[1,\infty )$

3.4. Goodness of fit test

Estimating Copula function parameters relies on the relationship between the rank correlation coefficient $\tau$ and parameter $\theta .$ Maximum likelihood estimation is commonly used to calculate Copula function parameters. However, it is still necessary to assess the goodness of fit of the joint distribution. Various test methods, such as the Ordinary Least Squares (OLS) and the Akaike information criterion method (AIC), can be employed to evaluate the goodness of fit. OLS is calculated as follows: (3) $$\mathit{OLS}={\left[\frac{1}{n}\sum _{i=1}^n{(F_{\mathit{enp}}-C)}^2\right]}^{1/2}(1\le i\le n)$$(3)

Where $F_{\mathit{enp}},C$ are the empirical and theoretical frequency value; n is the series length.

AIC is calculated as follows: (4) $$\begin{array}{c}MSE=\frac{1}{n-1}\sum _{i=1}^n{(F_{\mathit{enp}}-C)}^2\mathrm{ }\mathrm{ }(1\le i\le n)\\ AIC=n\hspace{0.17em}\ln \hspace{0.17em}\left(\frac{\mathit{MSE}}{n}\right)+2m\end{array}$$(4) where m is the number of parameters.

The sizes of OLS and AIC reflect the degree to which the Copula function fits the original variable data, with smaller values indicating a better model fit.

3.5. Joint recurrence period

The joint recurrence period of drought is defined as the duration between the first and second occurrences of a drought. It is worth noting that the joint recurrence period for variables that exceed or equal a specific predetermined threshold is shown below (Li et al. 2022). (5) $$T_i=\frac{N}{n\mathrm{ }[1-F(X_i)]}$$(5) where T_i is the univariate recurrence period; F(X_i) is the marginal distribution function; N is the sample length; and n is the number of occurrences greater than a given sample.

Two-dimensional joint recurrence period expression is as follows: (6) $$T_{\alpha }=\frac{N}{n\cdot P(X_1\ge x_1\cap X_2\ge x_2)}=\frac{N}{n\cdot (1-C(u,v))}$$(6)

Three-dimensional joint recurrence period expression is as follows: (7) $$T_{\alpha }=\frac{N}{n\cdot P(X_1\ge x_1\cap X_2\ge x_2\cap X_3\ge x_3)}=\frac{N}{n[1-C(u,v,w)]}$$(7) where $T_{\alpha }$ is the joint recurrence period. $C(u,v)$ and $C(u,v,w)$ are the joint distribution function. u,v and w are the marginal distribution function of each variable.

In summary, this study uses precipitation data and multi-scale SPI values for 1961–2010 in Henan Province as the basis for analysing drought risk in Henan Province from a global perspective using the method described above, overcoming the loss of information during data conversion and the limitations of the analysis results. In addition, the data used in this paper were cleaned and validated for consistency and completeness.

4. Results and discussion

4.1. Mechanisms of drought formation

Anomalies in atmospheric circulation or insufficient water vapor can lead to decreased precipitation, which is a key indicator of the onset of drought (Wang T et al. 2023; Zhang et al. 2023). Droughts are typically characterized by five distinct stages, namely gestation, initial, buffer, development and relief. Anomalies in atmospheric circulation and insufficient water vapor serve as the gestation process for drought, while reduced precipitation is the direct cause of drought and marks the beginning of the process. If the rate of evaporation increases and water vapor is replenished in time, precipitation may increase. However, if water vapor cannot be replenished in time, precipitation will continue to decrease, and the drought process will continue for a longer period of time, which represents the buffering and development period of drought. After passing through the gestation period, the initial period, the buffer period, and the development period, with the cessation of abnormal circulation, the increase in precipitation and the supply of water vapor from external sources will lead to an increase in precipitation, which will eventually lead to the relief of drought. Figure 3 shows the process of drought formation and development.

Figure 3. The process of drought formation and development.

The study of the processes that shape drought mechanisms can provide a fundamental understanding of drought. In addition, the influence factors of drought are complex, and its driving mechanism, dynamic response and propagation process are not very clear. Therefore, the study of the formation process, driving mechanism and dynamic response of drought can reveal the propagation mechanism of drought to a certain extent, and has certain guiding significance for drought early warning.

4.2. Multi-timescale analysis of drought

Given the complexity of the drought formation process and the multitude of factors involved, we calculate and transform drought characteristic indicator variables based on meteorological data from 1961–2010. Besides, we conduct a multi-scale analysis of drought in Henan Province using different SPI values (1, 3, 6, 9, 12, 24) to investigate the drought phenomenon. Table 4 shows the fitted trend lines for SPI values at different time scales. The SPI is a representation of drought at different time scales and is one of the important indicators for characterizing drought, and SPI values at different time scales have different advantages. With smaller scale SPI values being more sensitive to drought response and have a higher sensitivity for identifying drought, which can characterize drought characteristics in a short period of time. In contrast, larger scale SPI values can accurately characterize the occurrence of drought over a longer time span. Therefore, using different scale SPI values in different analyses are crucial to improve the accuracy of drought research.

Table 4. SPI values fitted to trend lines.

Trendlines SPI Scale	Fitting formula
SPI1	y = 0.027
SPI3	y=−0.0057
SPI6	y=−0.0257
SPI9	y = 0.0001x−0.0698
SPI12	y = 0.0002x−0.1081
SPI24	y = 0.0002x−0.1213

From Table 4 and Figure 4, it can be seen that the SPI values at different time scales in Henan Province show an oscillatory process that varies over time, with a consistent trend of variation, with values fluctuating within the range of [-2.5, 2.5]. In particular, there are significant differences in the oscillation and changes of SPI values at different scales, with smaller scales showing more frequent and dramatic fluctuations. This observation suggests that the SPI values with smaller scales are more sensitive in identifying the occurrence of drought in Henan Province, and can promptly detect the process of drought occurrence within a short period. As the scale of the SPI value increases, the trend line of the SPI value oscillation becomes more stable, reflecting the cumulative effect of precipitation in the early stages on the process of drought occurrence.

Figure 4. Multi-scale time series trend graph of SPI in Henan Province.

Therefore, in studies, SPI values at the appropriate scale should be selected for the problem under study. At the same time, the application of multi-timescale SPI values not only allows for the timely identification of the degree of drought in the study area, but also provides an integrated picture of the effect of cumulative precipitation on drought. This has also been verified in the literature (Mehr and Vaheddoost 2020).

4.3. Drought evolution characterization

Based on the data of drought characteristic variables, the evolution law of drought duration, severity, and kurtosis is visualized in Figure 5. It can be concluded that Zhumadian has a high frequency of drought occurrence with long duration and high intensity drought, followed by Pingdingshan, Luohe, and Zhoukou, while the drought intensity in Xinxiang, Puyang, Hebi, and Anyang is relatively lighter compared to other regions. Generally speaking, droughts are more likely to occur in central and southern Henan Province, while the frequency of droughts in northern Henan Province is relatively low. Dividing the region along the Yellow River, the frequency of drought occurrence is higher in the southern area of the Yellow River than in the northern area.

Figure 5. Graphs of the evolution of drought duration, severity and kurtosis.

The evolution of the different drought characteristic variables visually demonstrates the drought evolution pattern in Henan Province. Through the analysis of this pattern, it is possible to grasp the dynamic evolution of drought in Henan Province, which can provide more comprehensive guidance to the relevant departments in their drought prevention work.

4.4. Drought risk analysis based on a two-dimensional copula function

4.4.1. Copula function preference

In order to improve the accuracy of drought risk analysis results in Henan Province, based on the analysis results, the SPI 1 with high drought identification sensitivity is selected as the data basis. The drought duration, severity, and kurtosis of 18 meteorological stations in Henan Province are extracted based on the runs theory, and they are set as random variables X, Y, and Z. Based on the two-dimensional Archimedean Copula function, the joint distribution of drought characteristic variables is constructed, and the optimal Copula function is selected based on the OLS and AIC, and a drought risk analysis model for Henan Province is established. The calculation results are shown in Tables 5–7.

Table 5. Parameter estimation and goodness of fit results (Drought duration-severity).

Copula Function	$\theta$	Correlation test		Goodness of fit test
		Kendall	Spearman	OLS	AIC
Gumbel-Copula	2.0121	0.5030	0.6859	0.0247	−319.245
Clayton-Copula	1.5434	0.4356	0.6075	0.0720	−238.863
Frank-Copula	6.3091	0.5300	0.7284	0.0171	−346.874

Table 6. Parameter estimation and goodness of fit test results (Drought duration-kurtosis).

Copula function	$\theta$	Correlation test		Goodness of fit test
		Kendall	Spearman	OLS	AIC
Gumbel-Copula	1.2560	0.2038	0.2989	0.0165	−347.104
Clayton-Copula	0.4694	0.1901	0.2809	0.0274	−309.218
Frank-Copula	1.9595	0.2099	0.3110	0.0167	−346.204

Table 7. Parameter estimation and goodness of fit test results (Drought severity-kurtosis).

Copula function	$\theta$	Correlation test		Goodness of fit test
		Kendall	Spearman	OLS	AIC
Gumbel-Copula	2.5243	0.6038	0.7918	0.0269	−310.594
Clayton-Copula	2.4998	0.5555	0.7422	0.0776	−231.455
Frank-Copula	7.9536	0.6009	0.8019	0.0484	−266.717

According to the results are shown in Table 5, the Kendall and Spearman correlation coefficients of the joint distribution constructed by drought duration and severity fluctuate around 0.5, indicating a strong correlation between the two. In further calculations, the Frank Copula function has the smallest OLS and AIC among the three functions, with OLS = 0.0171 and AIC=−346.874. Therefore, the Frank Copula function is the most accurate choice for the joint distribution function of the two variables.

According to the results are shown in Table 6, it can be observed that the Kendall and Spearman coefficients for the joint distribution of drought duration and drought kurtosis show a slight fluctuation around 0.2. This indicates a weak correlation between the two. Additional calculations are performed to obtain a more precise selection of the joint distribution. Specifically, the OLS and the AIC for the Gumbel copula function were calculated and found to have minimum values of OLS = 0.065 and AIC = −347.104 respectively. As a result, it is recommended to use the Gumbel copula function in the joint distribution of the two variables for greater accuracy.

According to the results are shown in Table 7, the Kendall and Spearman coefficients of the joint distribution constructed by drought severity and drought kurtosis fluctuate around 0.6. In addition, the Spearman coefficient of the Frank Copula function reaches 0.8019, indicating an extremely strong correlation between drought severity and drought kurtosis. In further calculations, the OLS of the Gumbel copula function is 0.0269 and the AIC is −310.594. Both the OLS and AIC are the smallest of the three functions. Therefore, the Gumbel Copula function is selected to construct the joint distribution of drought severity and kurtosis for greater accuracy.

To further determine the fitting effect of the joint distribution of drought characteristic variables, a Q-Q plot is used to analyse the joint distribution between the variables. The Q-Q plot can intuitively express the fitting situation between the empirical and theoretical frequencies of the joint distribution. The empirical-theoretical joint distribution fitting plot for drought characteristic variables is shown in Figure 6.

Figure 6. Fitting the empirical joint distribution to the theoretical joint distribution: (a), (c), (e) are QQ graphs of drought duration-severity, duration-kurtosis, and severity-kurtosis, respectively and (b), (d), (f) are theoretical and empirical frequency fits of drought duration-severity, duration-kurtosis, and severity-kurtosis, respectively.

According to Figure 6(a,b), the Frank Copula function has a high goodness of fit between the theoretical and empirical frequencies of the joint distribution of drought duration and drought severity, with a correlation coefficient of R²=0.9985 and RMSE = 0.0098. Similarly, based on Figure 6(c,d), the Gumbel Copula function has a high goodness of fit between the theoretical and empirical frequencies of the joint distribution of drought duration and drought kurtosis, with R²=0.9990 and RMSE = 0.0089. Furthermore, Figure 6(e,f) shows that the Gumbel Copula function has a high goodness of fit between the theoretical and empirical frequencies of the joint distribution of drought severity and drought kurtosis, with R² = 0.9984 and RMSE = 0.0116. The correlation coefficients of the joint distribution of drought characteristic variables are all close to 1, and the error is less than 1%, indicating a high accuracy of different Copula functions in fitting the joint distribution of drought characteristic variables.

4.4.2. Joint probability distribution

In view of the above calculation results, the optimal joint distribution between drought duration and drought severity is the Frank Copula function with a parameter value of 6.309, while the optimal joint distributions between drought duration and drought kurtosis, and between drought severity and drought peak are both the Gumbel Copula function with parameter values of 1.256 and 2.524, respectively. On this basis, the probability density functions and empirical distribution functions of between the drought characteristic variables are shown in Figure 7. The Frank Copula function has no restrictions on the correlation between random variables (Kavianpour et al. 2018; Fan et al. 2022), making it suitable for both positively and negatively correlated variables. The image feature of the Frank Copula function is symmetric in the tails and asymptotically independent between the upper and lower tails. The probability density function graph of the Gumbel Copula function has a "J"-shaped distribution, with a high upper tail and a low lower tail, and is not symmetric between the upper and lower tails.

Figure 7. The Probability density and empirical distribution function graphs: (a), (c), (e) are graphs of the joint density function between the drought characteristic variables, respectively and (b), (d), (f) are graphs of the joint distribution function between the drought characteristic variables, respectively

Based on Figure 7(a,b), it can be observed that the probability density distribution of the two random variables of drought duration and severity have symmetrical tails and relatively stable upper and lower tail dependence, with the probability density concentrated around 4–5 months. On the other hand, Figure 7(c–f) shows that the probability density distribution of drought duration and kurtosis, drought severity kurtosis are highly sensitive to changes in the upper tail, and even subtle changes in the upper tail can be identified and captured, making it more sensitive for drought identification. Furthermore, in Figure 7(c,d), the probability density is concentrated around 3–4 months, while in Figure 7(e,f), it is concentrated around 10–15 months. In summary, the construction of a risk analysis model for drought characteristic variables using the Frank Copula and Gumbel Copula functions is of great significance for improving the accuracy of drought risk analysis in Henan Province.

Based on the joint distribution function of drought characteristic variables, it can be observed that the joint distribution of drought duration and severity, drought duration and kurtosis, and drought severity and kurtosis show consistent trends. As the probability of the marginal distribution of drought characteristic variables increases, the probability of the joint distribution also increases, and the growth trend is relatively stable. In summary, the probability values of droughts occurring in 3–5 months and 10–15 months are more concentrated. Therefore, it can be inferred that droughts in Henan Province are more likely to occur in 3–5 months and around 10–15 months, indicating a tendency towards medium-term and annual droughts, with a relatively stable development trend during the study period, and the frequency of cross-seasonal droughts is relatively high.

4.4.3. Joint recurrence period

The joint recurrence period of drought refers to a drought duration (severity, kurtosis) that is greater than or equal to the given drought duration (severity, kurtosis) value, which means that the drought of that duration (severity, kurtosis) is considered to have recurred once. In order to accurately understand the drought situation in Henan Province, the joint recurrence period and contour distribution graph were plotted based on the joint distribution function of the drought characteristic variables. The results are shown in Figure 8.

Figure 8. Graphs of the joint recurrence period of drought characteristic variables: (a), (c), (e) are joint recurrence period graphs between drought duration-severity, duration-kurtosis, and severity-kurtosis, respectively and (b), (d), (f) are joint recurrence period contour graphs between drought duration-severity, duration-kurtosis, and severity-kurtosis, respectively.

Analysis of Figure 8(a,b), it can be observed that for the joint distribution of drought duration and severity, when the drought duration and severity are less than 1 and 1000, respectively, the joint recurrence period is less than 4 months. When the drought severity is between 1000 and 1700, and the drought duration is between 1 and 1.8, the joint recurrence period is between 4 and 6 months. When the drought severity is between 1700 and 2200, and the drought duration is between 2 and 2.2, the joint recurrence period is approximately 7 months. When the drought severity is greater than 2500, and the drought duration is greater than 2.5, the joint recurrence period increases rapidly and reaches between 12 and 16 months.

From Figure 8(c,d), we know that for the joint distribution of drought duration and drought kurtosis, when drought duration and kurtosis are less than 1 and 800, respectively, the joint recurrence period is between 0 and 4 months. When the drought kurtosis is between 800 and 1000 and the drought duration is between 1 and 2, the joint recurrence period increases to between 4 and 6 months. Furthermore, when the drought kurtosis increases to between 1200 and 1400, and the drought duration is between 2 and 2.5, the joint recurrence period further increases to approximately 9 months. Finally, for drought events with higher kurtosis and longer duration, the joint recurrence period increases rapidly, reaching a peak of between 12 and 16 months.

Based on Figure 8(e,f), it can be seen that for the joint distribution of drought severity and drought kurtosis, when either the drought severity is less than 1000 or the drought kurtosis is less than 800, the joint recurrence period is less than 4 months. When the drought severity is between 1000 and 1800 and the drought kurtosis is between 800 and 1200, the joint recurrence period is between 4 and 6 months. When the drought severity is between 1800 and 2200 and the drought kurtosis is between 1200 and 1300, the joint recurrence period is 7 months. As the drought kurtosis value increases, the joint recurrence period increases rapidly and reaches 15 months.

In summary, the joint recurrence period of drought is closely related to the changes of drought characteristic variables, and the length of its period shows a positive correlation with the trend of these variables. The drought risk analysis model based on the Frank Copula and Gumbel Copula functions shows that when the drought duration is less than 1, and both the drought severity and kurtosis are less than 1000 and 800, respectively, the joint recurrence period of drought is less than 4 months, which means that droughts occur more frequently and short-term droughts are very likely to occur. When the drought duration is between 1 and 1.8, the drought severity is between 1000 and 1800, and the drought kurtosis is between 800 and 1000, the joint recurrence period of drought is 4–6 months. At this time, it is easy for medium and long-term drought to occur, and cross-seasonal drought is also likely to occur, which will significantly affect people’s lives. As the values of drought characteristic variables continue to increase, when the drought duration is greater than 2.5, the drought severity is greater than 2200, and the drought kurtosis is greater than 1300, the joint recurrence period of drought suddenly increases, which indicates that the probability of drought occurring during this period is very low.

Compared to risk analysis models based on grey theory and principal component analysis (Xu et al. 2021; Yu et al. 2021), the copula function can construct joint distributions between drought characteristic variables without loss of data information, while dealing with non-linear and non-normal dependencies between drought characteristic variables. This makes the model more applicable to drought risk analysis and provides greater accuracy when dealing with complex data. In addition, it can provide accurate data for disaster warning and economic loss estimation.

4.5. Drought risk analysis based on a three-dimensional copula function

Based on the extracted values of drought characteristic variables, with drought duration as the random variable X, drought severity as the Y, and drought kurtosis as the random variable Z, a joint distribution of drought characteristic variables is constructed using the three-dimensional Archimedean Copula function. Considering strong positive correlation between the drought characteristic variables in the aforementioned study, the three-dimensional Gumbel-Hougaard Copula function is chosen for the joint distribution construction. The results of the goodness of fit test and parameter estimation are shown in Table 8.

Table 8. Results of the goodness of fit parameter calculation.

Copula function	Duration-severity-kurtosis
	OLS	AIC	${\theta }_1$	Kendall	Spearman
Gumbel-Hougaard Copula	0.023	−325.633	1.940	0.550	0.740

According to Table 8, in the joint distribution of drought duration-severity-kurtosis constructed by the Gumbel-Hougaard Copula function, the Kendall and Spearman correlation coefficients are 0.55 and 0.74, respectively, both of which are greater than 0.5. This indicates a strong positive correlation between the drought characteristic variables. In addition, the OLS value of the joint distribution is 0.023, and the AIC value is −325.633. The joint distribution constructed by this function has a high goodness of fit.

Based on the above function formulas, a three-dimensional joint distribution is constructed, and a drought risk analysis model is established for risk analysis. The results are shown in Figure 9 and Figure 10. From the analysis of the three-dimensional joint probability trend graphs (Figure 9), we know that droughts have become more frequent in Henan Province between 1961 and 2010, with droughts occurring every year, and most of them are mild droughts. However, severe droughts occurred in 1972, 1980, 1987, 1992, and 2000. In addition, as the marginal distribution probability of any of the drought characteristic variables increases, the joint distribution probability shows a tendency to increase. When two drought characteristic variables have a certain probability value, the joint distribution probability of the three variables increases and tends to a specific value as the other variable increases, while the value of the joint distribution probability eventually tends to 1 when all three variables increase at the same time.

Figure 9. The joint probability trend.

Figure 10. Three dimensional joint recurrence period.

The analysis in Figure 10 shows a strong positive correlation between the three variables. When the joint probability of drought duration-severity is between 0 and 0.45, or the joint probability of drought severity-kurtosis is between 0 and 0.4, the drought recurrence period is 5 months, and droughts occur frequently during this period, mostly mild droughts. When the joint probability of drought duration-severity is between 0.45 and 0.78 or the joint probability of drought severity-kurtosis is between 0.45 and 0.8, the drought recurrence period is 10 months, and moderate droughts are more likely to occur during this period, but less frequent. Moreover, when the joint probability value of drought duration-severity is greater than 0.5, the joint recurrence period is at 45 months, and when the value of drought severity-kurtosis is greater than 0.7, the joint recurrence period is 35 months. Therefore, with an increase in the joint probability of drought duration-severity, duration-kurtosis and severity-kurtosis, the contour values of the joint recurrence period become larger and the spacing becomes smaller. This shows that as the joint probability of any two of the three drought characteristic variables increases to a certain extent, the drought density increases, and severe droughts are more likely to occur.

In summary, the three-dimensional copula function constructs a multidimensional joint probability distribution by considering the marginal distributions and correlation structure of arbitrary drought variables. It overcomes the limitations of traditional approaches in quantifying drought risk analysis in Henan Province. The flexible and diverse form of this function provides new means to study drought risk analysis (Wei et al. 2023). Furthermore, the traditional two-dimensional copula function has limitations in revealing the correlation between drought variables in Henan Province and drought risk. The three-dimensional copula function model can comprehensively consider more drought variables and fully reveal the response relationships between them (Seyedabadi et al. 2020).

Moreover, in the analysis of drought risk in Henan Province, the three-dimensional copula function has a better fit to the drought variables than the two-dimensional copula function. This indicates that the multivariate joint distribution has greater advantages in the analysis of drought risk in Henan Province. However, further research on parameter estimation and function selection is still needed to improve the model. Furthermore, by combining the results of the two-dimensional and three-dimensional copula function models, it is found that Henan Province is highly prone to seasonal drought, with mild drought dominating. This research result can provide theoretical support for the drought warning work of government departments.

5. Conclusions

Drought is influenced by many factors and its driving mechanisms are complex. In order to overcome the limitations of traditional studies in quantifying drought risk and to reveal the response relationships between different drought characteristic variables, this study uses runs theory to extract drought duration, severity and kurtosis from standardized precipitation index (SPI) in Henan Province and constructs drought risk analysis models based on multivariate Copula function. The results are as follows.

1. The SPI shows good sensitivity in drought identification, with higher sensitivity at smaller scales. Between 1961 and 2010, Henan Province was prone to mild drought. However, severe drought occurred in 1972, 1980, 1987, 1992, and 2000.

2. The optimal fitting functions for the joint distribution between drought duration-severity, duration-kurtosis, and severity-kurtosis are Frank copula and Gumbel copula, with R² values of 0.9985, 0.9990, and 0.9984, respectively.

3. The drought recurrence period is positively correlated with the values of drought characteristic variables. There is a higher probability of drought recurrence period occurring between 4 to 6 months. This suggests that droughts are more frequent in Henan Province, with a high probability of short-term and cross-seasonal droughts.

4. Analysing the different areas in the study region, droughts are more frequent in central and southern Henan and less frequent in northern Henan, divided along the Yellow River, droughts are more frequent in areas south of the Yellow River compared to areas north of the Yellow River.

5. In summary, the three-dimensional Copula function fits the drought variables better than the two-dimensional Copula function in the analysis of drought risk in Henan Province, which indicates that the multivariate joint distribution has greater advantages in the analysis of drought risk in Henan Province. However, further research on parameter estimation and function selection process is still needed to improve the model.

Authors’ contributions

Yunliang Wen: investigation, conceptualization, methodology, data curation, software, formal analysis, writing-original draft, visualization. Liwei Zhou: conceptualization, methodology, data curation, software, formal analysis. Ling Kang: supervision, validation, resources, writing-review & editing, funding acquisition. Hao Chen: data curation, software, formal analysis, validation. Jinlei Guo: data curation, software, visualization

Disclosure statement

No conflict of interest exits in the submission of this manuscript, and manuscript is approved by all authors for publication.

Data availability statement

The data that support the findings of this study are available on request from the corresponding author, Yunliang Wen.

Additional information

Funding

This research was funded by the National Key R&D Program of China (Grant No. 2022YFC3002704).