Estimation of multivariate design quantiles for drought characteristics using joint return period analysis, Vine copulas, and the systematic sampling method

Around the world, the hydrological processes have been altered and extreme hydrological events have occurred frequently (Neupane & Kumar 2015; Tenagashaw et al. 2022) due to both climate change (atmospheric circulation, precipitation, and evaporation) and anthropomorphic factors (land use, population growth, industrialization, and socioeconomic development) (Legesse et al. 2004; Kundu et al. 2017; Yang et al. 2019). Among the extreme hydrological events, drought characterized by long duration has been considered as one of the most widespread, complex, and serious disasters (Rojas 2018; Mohammed et al. 2022). It is estimated that global economic losses caused by droughts exceed 6 billion dollars (Wilhite 2000; Yang et al. 2018). Droughts have had a serious effect on agricultural production, the quality of residents’ lives, and the fragile ecological environment (Hou et al. 2019; Hu et al. 2021). How to reduce drought damage has become an urgent issue (Mishra & Singh 2010; Logar & van den Bergh 2013; Saini et al. 2021).

Drought can be divided into a meteorological, agricultural, hydrological, and socioeconomic drought (lack of precipitation, soil moisture, surface or groundwater resources, and commodities, respectively) summarized by the American Meteorological Society (Heim 2002; Sadegh et al. 2017). Hydrological drought is characterized by low surface, reservoir, and underground water levels (Hateren et al. 2019). This type of drought is considered ‘thorough’, i.e., it most accurately reflects the basin water shortage and complexity of drought resistance (Keyantash & Dracup 2002; Wu et al. 2019). Hydrological drought evaluation is an important basis for determining the construction scale of drought relief projects and formulating operational management strategies (Zhang et al. 2022). It is essential to conduct research about how to determine the drought relief project scale based on hydrological drought evaluation for reducing drought damage.

Estimation of the return period is a criterion for measuring the expected recurrence interval of a hazard and quantifying the design quantiles (usually defined as the values of a variable under a given return period) of key variables in a water system, which is a process necessary for reasonably determining the scale of a water project (Salvadori 2004; Salvadori et al. 2011a, 2011b). Univariate statistical analysis is the traditional method to estimate the return period for project scale design (Gräler et al. 2013). The project scale is determined based on the inverse function of the single variable marginal distribution function and the return period (Salvadori & De Michele 2007; Volpi & Fiori 2014). However, hydrologic events are often characterized by the joint behavior of multiple variables (Tootoonchi et al. 2021). For example, there are three main drought characteristics responsible for the multiple negative effects of drought: drought duration, peak, and severity (Yang et al. 2018). The main flood characteristics are duration, peak, and volume (Grimaldi & Serinaldi 2006). Different characteristics correspond to different properties of the natural disaster. Univariate statistical analysis is not sufficient to reflect the complex and interactive relationship among multiple variables (Sadegh et al. 2017; Ballarin et al. 2021). Hydrologic event return period analysis involving the use of multiple characteristics may be more comprehensive than that based on a single characteristic (Salvadori & De Michele 2007; Tootoonchi et al. 2021). It is vital to determine the project scale based on a multivariate joint return period involving multiple drought characteristics.

Copula is a popular tool to describe the internal complex linear or nonlinear relationships among different variables (Dodangeh et al. 2017). It is a candidate for the multivariate joint project design which is of simple, flexible, and margin-free characteristics (Gyasi-Agyei 2012). Multivariate copulas commonly used include Gaussian, Student-t, and Archimedean copulas. However, they all possess shortcomings. For example, the marginal distribution functions for the different variables of Gaussian copulas must be normal (Pereira et al. 2017). Student-t copulas do not have asymmetric dependencies (Bhatti & Do 2019). Multivariate Archimedean copulas can only model positive dependence. In addition, in three dimensions, the asymmetric schemes are only useful when two correlation structures are equal and lower than the third correlation (Grimaldi & Serinaldi 2006). Against this background, the Vine copula was proposed by Bedford & Cooke in 2001 based on pair copulas. Vine copulas decompose the multivariate joint function into a series of stagewise mixing of bivariate copulas, allowing the former to model more flexible dependence structures and reflect asymmetric dependence (Pereira et al. 2017; Lü et al. 2020; Dixit & Jayakumar 2022).

Several studies have employed the Vine copulas to overcome univariate limitations in project design. For example, Gräler et al. (2013) designed hydrology characteristics of discharges composed of annual maximum peak discharge, volume, and duration. Vandenberghe et al. (2012) estimated the design quantiles of extreme discharge events including three variables: peak discharge, volume, and duration. However, there is still a limited number of studies focusing on using the Vine copulas to estimate joint design quantiles of drought in which three-dimensional multiple variables such as duration, peak, and severity are considered simultaneously. Therefore, the main purposes of this study are to (1) quantify the complex relationship among hydrological drought characteristics in three dimensions (including drought duration, peak, and severity) based on Vine copulas; and (2) compare the design quantiles of the scale of drought relief project based on the univariate and multivariate framework.

This is a case study of the Wei River Basin (WRB). The WRB is an important economic, agricultural, and cultural region. However, this region has suffered from frequent droughts (Yang et al. 2018). With recent rapid social and economic development, the water shortage has become more acute and the ecological environment has been more fragile (Chang et al. 2016). Drought is one of the most serious disasters restricting sustainable development in the WRB (Chang et al. 2016). It is vital to determine the scale of the drought relief project for reducing the drought damage in the WRB. To the best of the authors’ knowledge, this is one of the early tries to design the drought relief project scale from a multidimensional perspective especially in particular for the WRB. Overall, this study could enhance the understanding of drought events and provide support for drought relief project design.

The Wei River is the largest tributary of the Yellow River, flowing through Gansu, Ningxia, and Shaanxi three provincial-level administrative regions (Yang et al. 2018). It is a typical arid and semi-arid region, situated in Northwest China. Annual precipitation, evaporation, and average temperature are approximately 559 mm, 850 mm, and 10.6 °C, respectively. It is a water-deficient region. In addition, water demand has increased due to rapid development. Thus, droughts occur frequently (Yang et al. 2018). In this study, the WRB is taken as a case study area, and the joint drought return period and joint design quantiles of drought characteristics are estimated (to be used for better water project planning and water resource management).

Next, based on the theory of runs and predetermined truncation levels (RAP = 0), three drought characteristics are identified: the consecutive period in which the RAP is less than 0, that is, the drought duration (D); the maximal monthly water deficit (the maximal difference between the runoff in a drought event period and the average runoff), that is, the drought peak (P); and the total water deficit in a drought event period, that is, the drought severity (S) (Shiau et al. 2007).

To derive the property of a single variable, five marginal distribution functions are utilized to fit a single drought characteristic. These marginal distribution functions include Gamma, Weibull, Exponential, Rayleigh, and Log-normal. The parameters of each marginal distribution function are estimated based on the maximum likelihood algorithm via MATLAB.

To determine the marginal distribution function with the best fit, the root mean square error (RMSE) is utilized for the goodness-of-fit test (Yang et al. 2019). The RMSE reflects the difference between the theoretical and empirical cumulative probabilities. A low RMSE value indicates a good match between the theoretical and cumulative probabilities. Theoretical cumulative probabilities are derived through the marginal distribution function, equations can be found in Yang et al. (2018). Empirical cumulative probabilities of the drought characteristics are estimated based on the Gringorten formula (Gringorten 1963; Kao & Govindaraju 2007).

In this study, three widely used bivariate copulas – the Gumbel, Clayton, and Frank – are employed to derive the bivariate joint distribution of the hydrological drought characteristics. The equations are displayed in Table 1. The Frank bivariate copula is symmetric, while the Gumbel and Clayton are asymmetric. The Gumbel best represents the upper tail dependence, while the Clayton represents the lower tail dependence best (Nelsen 2006). The optimal parameters corresponding to the three bivariate copulas are first estimated using the RMSE. The joint theoretical and empirical cumulative probabilities can be found in Yang et al. (2018). Next, the copulas whose RMSE values are the lowest are considered optimal.

Table 1

Distribution functions of three bivariate Archimedean copulas

Copulas type .	Distribution function C(u,v) .	Scope of the parameter .	Generator .
Gumbel
Clayton
Frank

Copulas type .	Distribution function C(u,v) .	Scope of the parameter .	Generator .
Gumbel
Clayton
Frank

View Large

In fact, there are lots of copulas that can be utilized for bivariate analysis such as total 26 copula families shown in Table 1 of Sadegh et al. (2017). The calculations of these copula families are freely available to the public. However, in the following Section 3.4, the ‘secondary return period’ (also called Kendall's return period) is utilized to compute the multidimensional joint return period which is based on the Kendall's distribution function (K_C) (Corbella & Stretch 2012; Rad et al. 2017). The general analytical expressions of K_C are known except for some special cases like that of bivariate extreme value copulas and some Archimedean copulas (Salvadori et al. 2011a). Therefore, in this study, three bivariate Archimedean copulas are chosen as candidates for the following analysis. Using more copula families to solve the multivariate analysis will be the authors’ future research purpose.

It should be noted that due to the variability of the model structure, in the three-dimensional case, there are a total of three decompositions (three variables appearing as the root nodes in the first tree). The results depend on the particular decomposition. It is particularly important to determine the root node in the first tree to further determine the relationships among the three variables. In this study, the maximum spanning tree is employed to select the pairwise copulas (Brechmann et al. 2012). The core idea is to calculate the Kendall's tau nonparametric association measure (calculations can be found in Fredricks & Nelsen 2007) of two randomly chosen variables. The variable corresponding to the maximum sum of Kendall's tau calculated from the two other variables is chosen as the root node in the first tree.

Based on the equations shown above, given an infinite combination of drought characteristics under a specific joint return period, the most likely design quantiles whose three-dimensional joint probability density function could be obtained. However, infinite data cannot be accessed. Here, the systematic sampling method is used to generate a series of drought characteristics for multivariate joint design quantiles estimation. The systematic sampling method extracts sample units at equal intervals. It is of simple, rapid, economical, flexible, ad convenient characteristics (Sethi 1965). Detailed steps to obtain the joint design quantiles based on systematic sampling and Vine copulas are as follows.

Step 1: Generate a cumulative probability distribution function for 1,000 groups with an equal interval of 0.001 for ⁠.

Step 2: Calculate the corresponding cumulative probability distribution function of based on the conditional probability (of the optimal bivariate copula (obtained in Section 3.3) to derive the under a given secondary return period condition. is directly calculated by p which is used to derive the secondary return period as shown in Equation (8).

Step 3: Further generate a cumulative probability distribution function for the root node with an equal interval of 0.001 under each generated condition. Gräler et al. (2013) pointed out that multivariate design quantiles need additional analysis as some design quantiles go beyond the standard in a single design event. In this study, the minimum cumulative probability distribution function for the root node is set to the univariate cumulative probability distribution probability under a given return period minus 0.01. The setting of 0.01 is referred to Table 2 in Gräler et al. 2013.

Step 4: Calculate the corresponding cumulative probability distribution function of and based on the conditional probability of the optimal bivariate copula to derive the and ⁠, respectively. Optimal copula type is obtained in Section 3.3 and the calculations can be found in Equations (4) and (5) and Figure 4.

Step 5: Calculate the marginal probability density function for the drought duration ⁠, severity ⁠, and peak based on the optimal marginal distribution function derived in Section 3.2.

Step 6: Calculate the copula probability density function and the three-dimensional joint probability density function of these generated combinations based on Equations (11)–(16).

Step 7: Find the maximal three-dimensional joint probability density function relating to specific joint return periods. Determine the ⁠, ⁠, and corresponding to the maximal ⁠. And calculate the drought duration, severity, and peak values based on the inverse functions of these optimal marginal distribution functions, the procedure is the same as Equation (3). These values are the multivariate joint design quantiles of duration, severity, and peak under given joint return periods.

After that, the aforementioned five marginal distribution functions are used to fit the duration, severity, and peak data. Based on the optimal parameters of each marginal distribution function, the goodness-of-fit between the empirical and theoretical cumulative probabilities are evaluated based on the RMSE values. The RMSE values indicate that in the JRB, the Exponential, Weibull, and Weibull are the most appropriate marginal distributions for the duration, peak, and severity, respectively. In the BRB, the Weibull, Weibull, and Weibull are the optimal marginal distributions for the duration, peak, and severity, respectively. In the MSWRB, the Exponential, Weibull, and Exponential distributions are optimal for the duration, peak, and severity, respectively.

Using the inverse functions of these optimal marginal distribution functions, the design quantiles can be estimated using traditional univariate drought return period analysis, i.e., using Equations (2) and (3). As Wang et al. (2019) pointed out, a basin's drought resistance capacity (representing the capacity to meet the water demand) is weak downstream of the Yellow River when the joint return period is about 6 years. As data period 1956–2015 utilized in this study and 1960–2013 data in Wang et al. (2019) is basically same, and the study area is of similarity (Wei River is the largest tributary of the Yellow River), in this study, authors focus on two drought return periods of approximately 6 years (5-year and 10-year) for further design quantiles analysis. The design quantiles of the duration, peak, and severity estimated using univariate drought return period analysis under these two drought return periods are shown in Table 2.

First, Kendall's tau correlation coefficients are calculated to determine the root node in the first tree (Table 3). Table 3 shows that the maximum sum of Kendall's tau with the two other variables corresponds to the duration, severity, and severity in the JRB, BRB, and MSWRB, respectively. Therefore, duration, severity, and severity are chosen as the root nodes for the JRB, BRB, and MSWRB, respectively.

Table 3

Kendall's correlation coefficient tau values for the three drought characteristics

Sub-basin .	Drought characteristics .	Duration .	Peak .	Severity .	.
JRB	Duration		0.6696	0.9261	1.5957
	Peak	0.6696		0.6552	1.3248
	Severity	0.9261	0.6552		1.5813
BRB	Duration		0.7198	0.9741	1.6939
	Peak	0.7198		0.7553	1.4751
	Severity	0.9741	0.7553		1.7294
MSWRB	Duration		0.6078	0.9442	1.5520
	Peak	0.6078		0.6469	1.2547
	Severity	0.9442	0.6469		1.5911

Sub-basin .	Drought characteristics .	Duration .	Peak .	Severity .	.
JRB	Duration		0.6696	0.9261	1.5957
	Peak	0.6696		0.6552	1.3248
	Severity	0.9261	0.6552		1.5813
BRB	Duration		0.7198	0.9741	1.6939
	Peak	0.7198		0.7553	1.4751
	Severity	0.9741	0.7553		1.7294
MSWRB	Duration		0.6078	0.9442	1.5520
	Peak	0.6078		0.6469	1.2547
	Severity	0.9442	0.6469		1.5911

View Large

The government of China has realized that current water conservancy projects are inadequate. Repairing existing infrastructure and initiating new projects is essential (Liu et al. 2013). In 2011, China issued the Central Document No. 1, focusing on accelerating and reforming water conservancy. In September of 2014, the Ministry of Water Resources, National Development and Reform Commission, Ministry of Finance, and Ministry of Agriculture issued the National Drought Relief Implementation Plan (2014–2016). Different types of drought relief emergency resources, such as reservoirs, water diversion projects, rubber dams, wells, and water cellars, were proposed and constructed in China to regulate water resources and mitigate the serious issue of water shortage (Liu et al. 2013; Li et al. 2018).

This study aims to compare the design quantiles calculated by univariate and multivariate frameworks for reliable drought relief project design. Besides the univariate design quantiles in Table 2 and multivariate design quantiles solved by Vine copulas in Table 4, the most likely design quantiles based on the multivariate framework solved by three-dimensional Archimedean copulas (Table 5) are also calculated for better contrastive analysis.

In practical drought relief project design, the total water shortage (which is correlated with drought severity) is the most important parameter used to determine the scale of the project. By comparing the design quantiles estimated via univariate drought return period analysis (Table 2) with the design quantiles estimated via multivariate joint return period analysis solved by Vine copulas (Table 4) and that solved by three-dimensional Archimedean copulas (Table 5), results show that under 5-year and 10-year drought return periods, severity design quantile (directly affecting the project scale) estimated by Vine copulas in the WRB considering three sub-basins is smaller than that by the univariate framework and that by three-dimensional Archimedean copulas. Possible reasons are as follows.

The univariate analysis only considers limited characteristics of the drought events rather than the complicated relationship. It may result in overestimation or underestimation of the results also pointed out by Shiau et al. (2006), Das et al. (2020), Fan et al. (2021), and Nguyen-Huy et al. (2022).

The three-dimensional copulas use one parameter to describe the relationship among multiple characteristics. Parameter estimation is relatively simple. However, this method restricts the consistency of correlation structure among variables, which may lead to insufficient correlation description among multiple elements (Kurowicka & Cooke 2006; Bedford et al. 2015). This method may be difficult to reflect the complex structure among the variables (Kurowicka & Cooke 2006; Bedford et al. 2015). It may further result in the design quantiles of three-dimensional copulas different from those derived by the univariate method, similar to the findings of Chen et al. (2018) and Shang et al. (2022).

Vine copulas decompose high-dimensional copulas into multiple bivariate copulas. It deals with the problems of high-dimensional data by reducing dimension (Pereira et al. 2017). In addition, the bivariate correlation structure is not limited. Vine copulas is more flexible than Archimedean copulas (Aas et al. 2009; Brechmann & Schepsmeier 2013). The result based on this method is more consistent with the inherent law of hydrological phenomena and can better reflect the correlation among various characteristics (Radfar et al. 2022). This method might lead to the design quantiles of severity in the WRB considering three sub-basins lower than that calculated by a univariate framework similar to the smaller design quantile of flood volume pointed out by Huang & Chen (2015).

The above results show that in traditional univariate drought return period analysis, the relationships among the three drought characteristics are not considered. Using Vine copulas provides a more comprehensive and reliable reference for drought return period analysis and design quantiles analysis. The finding that the drought severity design quantiles in the WRB calculated by three-dimensional Archimedean copulas and univariate drought return period analysis are larger than those calculated by Vine copulas indicates that future drought relief projects should be designed to have a larger total storage capacity, which is related to the water supply capacity, to reduce the water shortage and better resist drought. While drought resistance is better, higher water supply capacity will result in a larger project scale and inevitably the increase of the project investment. Further, high depreciation costs are associated with large, expensive infrastructure. For the case study area, WRB, from the characteristic's simulation of drought events perspective, the authors state that using the design quantiles calculated by Vine copulas might be more reliable to provide a basis for the construction of drought relief projects. Vine copulas calculation results not only consider the relationship among multiple variables but also achieve a balance between economic input reduction and project safety standard maintenance. This method provides a new way to estimate the drought relief project scale and provides a basis for water resources management decision-making.

It is worth noting that the most likely design quantiles may not be the most severe combinations relating to a specific return period (Sadegh et al. 2018). Relating to a given return period, all possible combinations constitute an isosurface. The design quantiles determined by multivariate framework provide a reference for engineering design. This study can be considered as a typical case study for its application in hydraulic engineering. In the process of design quantiles calculation, there are many uncertainties (detailed discussions shown in Section 5.2: Limitations and potential future work), which may lead to the deviation of the results. The authors state that the influence of various uncertainties on the design quantiles should be analyzed comprehensively in the future. By considering various uncertainties, inherent laws of hydrologic events (Volpi & Fiori 2012), project investment and other restrictive factors, the most appropriate design quantiles should be calculated scientifically and reasonably by comparing and verifying through various methods. This needs to be enhanced urgently and is our future research focus.

This study focuses on the design of hydrological infrastructure to alleviate drought. The authors assess the total storage capacity (corresponding to the water supply capacity) of the drought relief project design. Methodology in Section 3 is universal. The methodology can be applied to other regions and to various types of drought. However, the authors acknowledge that this study has several limitations.

First, in this study, the authors focus on surface water. It is also necessary to consider other water resources, such as groundwater, recycled wastewater, and transferred water, in our analysis. It is also necessary to consider other types of projects, such as those meant to increase the accuracy of drought forecasting, establish a drought-resistant decision system, and develop water-saving integrated technology.

Furthermore, this study assumes the marginal distribution function of each drought characteristic does not change during the whole data period. It means that this study assumes the data are stationary. However, with the impact of climate change and human activities, the hydrological sequence has changed (Serinaldi & Kilsby 2015). Without considering the nonstationarity may lead to the deviation and uncertainty of the results (Tootoonchi et al. 2021). Nonstationarity is not considered in this study for the following reason. The change point of the WRB is basically around 1995 (Yang et al. 2020). Considering the nonstationarity will result in the short data period before and after the change point. Shorter data are less representative of the population. This also leads to larger uncertainty and further affects management decision-making (Tong et al. 2015; Sadegh et al. 2018). If longer data are available, it is of great importance to estimate the design quantiles considering the nonstationarity under the changing environment.

Besides the nonstationarity uncertainty, multivariate design quantiles also face other uncertainties. For example, the uncertainty of univariate marginal distribution function parameter fitting (Xu et al. 2009; Sadegh et al. 2018), the uncertainty of the parameter estimation and type of the copulas (Zhang et al. 2015; Sadegh et al. 2018), uncertainty of sampling (Dung et al. 2015). Further research should focus on evaluating these uncertainties quantitatively and consider them reasonably in engineering design.

Drought causes significant economic loss worldwide. The construction of drought relief projects is an effective way to mitigate the impacts of drought. In this study, design quantiles of the water supply capacity (the key factor affecting the project scale) are estimated using return period analysis. Larger-scale projects require a large economic input, while at a smaller scale, it is difficult to reduce the water shortage. A proper return period calculation is essential for better project design. In traditional drought return period analysis, only a univariate variable-drought severity (which quantifies the magnitude of the total water shortage without considering other drought characteristics) is considered. Joint return period analysis better reflects the complex relationships among the different drought characteristics. Therefore, the main purpose of this study is to estimate the joint design quantiles of the drought characteristics using multivariate joint return period analysis for reliable design policy-making.

This study proposes a framework to achieve this goal. The first step is to generate the first conditional cumulative distribution function in Tree 2 of Vine copulas. Then, the second conditional cumulative distribution function in Tree 2 is calculated based on conditional probability. Next is to further generate the cumulative probability distribution function for the root node drought characteristic. The fourth step is to calculate the cumulative probability distribution functions of the other two characteristics still based on the conditional probability. The fifth and sixth steps are to calculate the marginal probability density functions of three characteristics, the copula probability density function and the three-dimensional joint probability density function. The final step is to find the most likely design quantiles relating to a specific return period.

The WRB (divided into three sub-basins) is chosen as the case study area to demonstrate the feasibility of the proposed framework. The results show that the design quantiles of the three drought characteristics calculated using univariate drought return period analysis are larger than those calculated using joint drought return period analysis by Vine copulas. Infrastructure designed using univariate drought return period analysis has a high total storage capacity but is expensive to build. In future drought relief project design, joint drought return period analysis calculated by Vine copulas should be employed to reduce the economic input into projects intended for drought relief while maintaining the project safety standard in the WRB.

Although the WRB is selected as the study area, our framework is highly transferable to other river basins. In addition, the proposed framework can be extended to n dimensions (n≥ 3), allowing for n-dimensional joint design quantiles analysis using a Vine copula structure that incorporates more than three event characteristics. Overall, our work increases the applicability of multivariate joint design quantiles analysis, which brings new sights for determining the project scale in drought-prone regions. However, in the future, more uncertainties (such as nonstationarity, parameter fitting, copulas type, and sampling uncertainty) should be involved. And a multi-method comparison pattern could better determine the appropriate design quantiles considering uncertainties, inherent laws, investment, and other restrictive factors that may influence design quantiles. These would be our objectives for future work.

This study is supported by the National Natural Science Foundation of China (52209030) and Fundamental Research Funds for the Central Universities (GK202207005). The authors thank the editors and reviewers for their helpful comments.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.