Abstract

The complex hydrological events such as storm, flood and drought are often characterized by a number of correlated random variables. Copulas can model the dependence structure independently of the marginal distribution functions and provide multivariate distributions with different margins and the dependence structure. In this study, the conditional behavior of two signatures was investigated by analyzing the joint signatures of groundwater level deficiency and rainfall deficiency in Naqadeh sub-basin in Lake Urmia Basin using copula functions. The study results of joint changes in the two signatures showed that a 90–135 mm reduction in rainfall in the area increased groundwater level between 1.2 and 1.7 m. The study results of the conditional density of bivariate copulas in the estimation of groundwater level deficiency values by reducing rainfall showed that changes in values of rainfall deficiency signature in the sub-basin led to the generation of probability curves of groundwater level deficiency signature. Regarding the maximum groundwater level deficiency produced, the relationship between changes in rainfall deficiency and groundwater level deficiency was calculated in order to estimate the groundwater level deficiency signature values. The conditional density function presented will be an alternative method to the conditional return period.

HIGHLIGHTS

  • The purpose of the present study is to investigate the conditional density analysis of bivariate copulas to estimate the groundwater level deficiency using the rainfall deficiency.

  • In this regard, the diagonal section of the copulas was used to reduce the complexity of the conditional density of the pairwise variables.

  • By using conditional density and combining it with copula simulation, the dependent values can be simulated and predicted.

  • Due to the lack of application of the conditional return period, the best alternative method is to use the proposed method based on conditional density.

  • The presented equation can be used as an alarm and monitoring system.

Graphical Abstract

Graphical Abstract
Graphical Abstract

INTRODUCTION

Drought in various fields, both directly and indirectly, will have adverse effects on different parts of the community and environment. The most important direct effect of drought is its negative impact on the quantity and quality of water resources. As the global climate changes, natural disasters caused by extreme weather events are becoming more frequent and leading to increasingly serious consequences around the world (Duan et al. 2019). The sustainable and efficient use of water is of vital importance, and the uneven distribution of water in the world means that it is especially important for water-stressed countries (Luo et al. 2020). On the other hand, it is important to understand the relative contributions of rainfall and its changes in runoff to sustainable water resources management (Lyu et al. 2019). Also, the extreme values of precipitation can have an impact on the landslide (Huo et al. 2020). With a reduction in rainfall for a long time, rangelands, forests, farms and other natural resources are directly depleted and damaged. Surface water and groundwater sources, fed by atmospheric rainfall, also suffer from water shortage. During and after drought, all agricultural activities and their associated facilities (urban, rural and industrial) are affected. Using a univariate hydrological frequency analysis, the probability of an event can be calculated. The univariate hydrological frequency analysis plays an important role in estimating the return period of flood or rainfall that are used to design structures such as dams, bridges, roads, highways, wastewater treatment facilities and industrial buildings.

But, since hydrological phenomena are usually described with two or more dependent variables, multivariate statistical analysis as well as dependency analysis is required. The most important problem of probabilistic multivariate analysis is providing the dependence structure for the associated random variables (Li & Zheng 2016). The multivariate distribution functions have been widely used to model two or more hydrological dependent variables and their dependence structures (Salvadori & De Michele 2007). The multivariate hydrological analysis mainly consists of three elements: (1) demonstrating the importance and usefulness of the multivariate framework, (2) fitting multivariate distributions to model hydrological phenomena and examine relevant parameters and (3) studying multivariate return periods or other hydrological analyses and simulations (Chebana & Ouarda 2011). In recent years, several multivariate methods have been introduced in hydrological and environmental applications. The most commonly used cumulative distribution function of copulas is Gaussian, but has a limitation that should be normal marginal distribution functions. The bivariate distribution was then proposed with non-normal marginal distribution functions (such as exponential bivariate distribution (Favre et al. 2002), bivariate gamma distribution (Yue et al. 2001) and bivariate distribution of extreme values (Adamson et al. 1999)). The bivariate distributions have been applied mostly in flood and rainfall frequency analysis and include bivariate normal distribution, bivariate exponential distribution, bivariate gamma distribution, multivariate distribution and bivariate extreme values distribution (Durrans et al. 2003; Shiau 2003; Beersma & Buishand 2004; Nadarajah & Gupta 2006; He et al. 2007; Wang et al. 2018). Since 2003, researchers such as Salvadori, Favre and Genest have pioneered the use of copulas in hydrology so that Favre et al. (2004) have used copulas to model hydrological variables in two basins in Quebec, Canada. They stated that the correlations between variables in hydrology could be modeled using this method and yielded better results than traditional univariate prediction methods.

Subsequently, the concept of copulas rapidly in various hydrological fields including flood frequency analysis (De Michele et al. 2005; Shiau 2006; Chebana & Ouarda 2007; Genest et al. 2007; Chen et al. 2011; Yang et al. 2019), multivariate analysis of rainfall characteristics (Salvadori & De Michele 2007; Zhang & Singh 2007; Kao & Govindaraju 2008; Zhang et al. 2013) and multivariate analysis of drought characteristics (Salvadori & De Michele 2004; Mirakbari et al. 2010; Mirabbasi et al. 2013; Zhang et al. 2015; Abdi et al. 2017; Ramezani et al. 2019) have been studied. In recent years, numerous studies on the use of copulas in the discussion of the multivariate analysis of meteorological and hydrological parameters have been reviewed. The analyses sought to provide selected copulas in the field of study and better results and comparisons with univariate analyses.

The purpose of the present study is to investigate the conditional density analysis of bivariate copulas to estimate the groundwater level deficiency using rainfall deficiency in Naqadeh sub-basin (NSB) located in Lake Urmia Basin (LUB). The proposed method is based on the conditional density and the diagonal section of the copula functions. By using conditional density and combining it with copula simulation, the dependent values can be simulated and predicted. Also, due to the lack of the application of the conditional return period, the best alternative method is to use the proposed method based on conditional density.

MATERIALS AND METHODS

Study area

Lake Urmia is 1,267 m above mean sea level in northwest of Iran. The study area is LUB and NSB. In this study, rainfall (1973–2015) and groundwater level (2000–2015) data in NSB have been used. Groundwater level is the depth that is defined from the Earth's surface to the water table. Naqadeh aquifer is an unconfined aquifer. The geographical location of the studied sub-basin is shown in Figure 1. In each sub-basin, a rain gauge station to introduce rainfall values and a piezometric well to indicate changes in groundwater level were selected. A summary of statistical properties of the studied data is also presented in Table 1.

Table 1

Statistical properties of the studied data on the annual scale

Sub-basinMean groundwater level (m)Average annual rainfall (mm)Piezometric wellRain gauge
Naqadeh 5.22 328.1 Tazegale Naqadeh 
Sub-basinMean groundwater level (m)Average annual rainfall (mm)Piezometric wellRain gauge
Naqadeh 5.22 328.1 Tazegale Naqadeh 
Figure 1

Location of the studied sub-basin in LUB.

Figure 1

Location of the studied sub-basin in LUB.

Entropy theory

In this study, entropy theory was used to select the piezometric well in the studied sub-basin. Since there are numerous piezometers in the studied sub-basin and the distribution of the stations is uniform throughout the sub-basin, entropy theory was used instead of Thiessen polygon or weighted average, and the station most closely resembles the studied sub-basin station is selected as the top station. Entropy as a scientific concept was first used by Clausius (1850) in thermodynamics. Then, in 1988, Boltzmann presented a probabilistic interpretation within the framework of statistical mechanics. The explicit relationship between entropy and probability was developed by Plank in the early 1900s (Harmancioglu et al. 1999). For more information, refer to Tahroudi et al. (2019b). It should be noted that a simulation tool is needed to perform the entropy method and estimate the interaction between the given parameters. In fact, instead of directly using the data from each well, the interactions between other wells should be incorporated. Therefore, in this study, the support vector regression (SVR) method optimized with the ant colony algorithm was used.

Extracting water resources signatures

The values of rainfall deficiency and groundwater level deficiency mentioned in this study as water resource signatures are extracted for 60-day consecutive duration based on the proposed equation of Tahroudi et al. (2019a). The values produced are hereinafter referred to as ‘signature’. The deficiency signature is referred to as an interval where the given parameter is lower than its mean long-term. The mean individual daily long-term period (MIDLP) for rainfall and groundwater level values for each 365 days of the year is derived as follows: 
formula
(1)
Accordingly, 365 MIDLP of rainfall and groundwater level values have been calculated for each station in first step. Then, daily time series data were subtracted from their own MIDLPs to calculate the new data abbreviated by ND. For example, the daily rainfall data for 1 October 1973 were deduced from the MIDLP of 1st October. In this way, this new data set (ND) equal to the number of time series of interest was created. The least value of the year was known as the rainfall deficiency for 1-day duration. Similarly, the minimum value obtained from summing the values of the rainfall deficiency through two successive days extracted from ND was regarded as the 2-day rainfall deficiency. This procedure continued until reaching the duration when the rainfall values were less than the average. Since this study aimed to reach the consecutive deficiency periods 1- to n-day (here maximum to 60-day), and it may have faced several gaps; therefore, the filling of these gaps may be carried out using a simple linear regression method (Tahroudi et al. 2019a).

Joint frequency analysis and copula theory

Copulas were first introduced by Sklar (1959). Copulas are functions that link multivariate distribution functions to bivariate marginal distribution functions (Nelsen 2006). Copulas are able to represent the dependence structure between two or more random variables and have recently emerged as a viable and efficient method for modeling the dependency on multivariate data (e.g., Joe 1997; Nelsen 2006). The advantages of using copulas over joint distribution of multidimensional models are (1) flexibility in choosing arbitrary marginal distribution functions and their dependence structure, (2) expanding to more than two variables and (3) the separate analysis of marginal distribution functions and their dependence structure (Salvadori et al. 2007; Serinaldi et al. 2009). The hydrological use of copulas has increased in recent years. In order to provide a precise definition of copulas, Sklar's theorem has been modified. If the random variables x1, …, xn follow the arbitrary marginal distribution functions F1(x1), …, Fn(xn) then there is copula C that combines these marginal distribution functions to the joint distribution function F(x1, …, xn) in the following: 
formula
(2)

If the marginal distribution functions of Fi(xi) are continuous, copula C is unique. Conversely, if C is a k-dimension copula, F is a n-dimension distribution function and Fi(xi), …, Fn(xn) is the corresponding marginal distribution functions. The copula families are varied and mainly include the following: (1) non-symmetric elliptical copulas (normal and t); (2) Archimedean (Clayton, Gumbel, Frank and Ali-Mikhail-Haq); (3) copulas of extreme values (Gumbel, Husler-Reiss, Galambos, Tawn and t-EV) and (4) other families (Plackett and Farlie–Gumbel–Morgenstern). Among different families of copulas, Archimedean and non-symmetric elliptical copulas are more popular in hydrology.

Conditional copulas

Suppose X and Y are random variables with U1 = FX(x) and U2 = FY(y). u1 and u2 are eigenvalues. For example, the conditional distribution function X, which is Y = y, can be expressed as follows: 
formula
(3)
Similarly, an equivalent formula for the conditional distribution functions for Y with respect to X = x can be obtained. In addition, for the conditional distribution function X, YY can be expressed as follows: 
formula
(4)

Similarly, an equivalent formula for the conditional distribution function for Y with respect to Xx can be obtained. Archimedean symmetric copula parameters can be calculated using moment methods by Kendall's tau correlation coefficient. For Archimedean non-symmetric copulas, the inference functions for margins (IFM) method or the maximum pseudo-likelihood (MPL) method can be chosen (Hult & Lindskog 2002; De Michele et al. 2007; Mirabbasi et al. 2012; Bezak et al. 2017).

Kendall's tau correlation coefficient

Kendall's tau correlation coefficient is widely used to measure the dependency of non-normal multivariate distribution functions. If X1 and X2 are two random vectors, Kendall's tau correlation coefficient can be written as follows (Kruskal 1958): 
formula
(5)

Goodness of fit tests

When validating a model fit, perhaps the first and most natural idea may be the scattering plot of pairs of which and are empirical and theoretical values, respectively (Genest & Favre 2007). The root mean square error (RMSE), mean absolute error (MAE), Nash–Sutcliffe (NSE) and Akaike's information criterion (AIC) are four common methods for the proper estimation of different copulas (Zhang 2005; Ma & San 2011). AIC can be estimated by either calculating maximum likelihood or mean square error of the model (Zhang & Singh 2006). 
formula
(6)
 
formula
(7)
 
formula
(8)
AIC values corresponding to maximum likelihood can be calculated using the following equation: 
formula
(9)
And AIC values of mean square error can also be calculated using the following equation: 
formula
(10)
where Pei and Pi are empirical and theoretical probabilities, respectively; N is the number of data, is the mean probability; m is the number of parameters and L is the maximum likelihood function value for the model. Copulas with lower AIC, MAE and RMSE and higher NSE values are better than other copulas studied. Therefore, these four criteria can play an important role in the selection of copulas.

Conditional density

An important point to consider the conditional returns period presented in various studies is the uncertainty of the results (Salvadori et al. 2007). Therefore, in this study, it has been attempted to investigate the conditional density equations of copulas to find an alternative method. In this way, the following algorithm is implemented at each stage:

At the 1st step, conditional density graph c(u, v) is plotted in the bivariate state for the studied variables. This graph is plotted for one of the values of u and v. For example, if the purpose of the analysis is to estimate the duration of the two signatures of rainfall deficiency and groundwater level deficiency, the values of v are equal to the rainfall deficiency signature. Because the rainfall deficiency signature precedes the groundwater level deficiency signature.

At the 2nd step, for each duration of rainfall deficiency signature, a graph of u is plotted.

At the 3rd step, in each graph, the maximum value of c(u, v) is selected and its corresponding values are estimated on X axis.

At the 4th step, these maximum values are in fact the duration of groundwater level deficiency signature which is affected by the corresponding duration of the rainfall deficiency signature. The copula conditional density parameter is also estimated as follows: 
formula
(11)
where C is equal to copula. The proposed method is based on the conditional density of copula functions. In this regard, the diagonal section of the copulas was used to reduce the complexity of the conditional density of the pairwise variables. Since there is a time lag between the occurrences of these two types of deficiency, the groundwater level deficiency can be estimated using the rainfall deficiency.

Copula diagonal section

The problem of finding copula C with a specific diagonal diameter has been discussed in several papers (Salvadori et al. 2007). The diagonal section of copulas in fact helps reduce the variate and makes it easier to understand the dependence structure without losing information. Therefore, diagonal copulas can be used at various stages of dependency modeling, and selection and diagnosis. Several designs are considered to simplify the dependence structure of the random pairwise variables. Some examples of such mappings are as follows: 
formula
(12)
 
formula
(13)
 
formula
(14)
where F represents the marginal distribution functions. In this study, a function of copula is considered, the result of which will be defined as the diagonal section of the copula. The diagonal section of copulas can be obtained from using the following equations: Function : I → I is defined as follows: 
formula
(15)
It is defined as diagonal section C. For each copula C, and for all I, . In addition, it also includes the following inequality: 
formula
(16)
For all values of I, . According to this definition, it is not difficult to show that any finite convex linear combination of Ci bivariate copula is itself a bivariate copula. In fact, for kN, C is calculated as follows: 
formula
(17)
That for all indices and . Then, C equals the suggested bivariate copula. Suppose C is an infinite set of copulas listed by a continuous parameter R. Now, suppose is the observation of a continuous random variable with the distribution function L, then: 
formula
(18)
It is not difficult to show that C is a copula. L usually refers to the mixing distribution of the family {} and C is defined by L as the convex sum. In the applications, it is often useful to consider a bivariate copula as restrictions to I2 on shared distribution functions whose margins are uniform rules in I. The following definition is a natural result of this fact. The problem of finding a copula C with a distinct diagonal section has been discussed in several papers (Salvadori et al. 2007). Its relation derives from the fact that if C is a copula that has two uniform variables of U and V in setting I, then the diagonal section of C contains information about the maximum behavior {U, V} and the minimum {U, V} of random variables. The diagonal section of the copula C is equal to the function : I → I, and regarding , the following properties are available: 
formula
(19)
The set of all functions with properties (1)–(4) is denoted by set D and called diagonally inside it. The equations in set D are called diagonal. The question that naturally arises is whether there is a copula for each diagonal whose diagonal section is equal to . is used to illustrate . As a first example, Bertino copula : I2 → I is as follows: 
formula
(20)
The second example of the diagonal copula : I2 → I is as follows: 
formula
(21)
If C is a copula whose diagonal section is , then C < . Moreover, if C is also symmetric, it can be concluded that C < . Kendall's values for Bertino and diagonal copulas are as follows: 
formula
(22)
 
formula
(23)
That may provide a way to fit the top copulas with the available data. Interestingly, Equations (22) and (23) have boundaries for the values of Kendall's tau of a copula C, when the values are known in quartiles C(i/4, i/4), i = 1, 2, 3 (Salvadori et al. 2007). Copulas (non-symmetric) with diagonal sections have recently been considered in studies (Salvadori et al. 2007), which has proposed a method for a family of copulas with respect to the lower and upper tail dependence coefficients. Specifically, with respect to and (the prescribed tail dependence coefficients), consider the linear function of the following: 
formula
(24)
where . It can then be shown that for each , all the family copulas {} are defined as follows: 
formula
(25)
is the upper tail dependence coefficient, and is the lower tail dependence coefficient.

RESULTS AND DISCUSSION

As mentioned, the data of rainfall and groundwater level in NSB in ULB were used in this study. The time series studied on the annual scale are presented in Figures 2 and 3. Figures 2 and 3 show the annual rainfall and the mean groundwater level in NSB, respectively. There was a lot of groundwater level data for each year. In Figure 3, mean groundwater level has been used for each year. As can be seen in Figure 2, the annual rainfall of the NSB during the statistical period is decreasing. The mean groundwater level has also increased over the studied years. This indicates that the increase of groundwater level in the study area over the studied years is due to reduced rainfall.

Figure 2

Annual rainfall values during the statistical period of 1973–2015.

Figure 2

Annual rainfall values during the statistical period of 1973–2015.

Figure 3

Mean groundwater level in the statistical period of 2000–2015.

Figure 3

Mean groundwater level in the statistical period of 2000–2015.

The entropy monitoring network optimization method was used to select the best piezometer in terms of information transfer that well reflects groundwater level in the study area. Based on the entropy method, the piezometers in the plain were classified and ranked, and finally, the best well in terms of information transfer was introduced as the selected piezometer. There are seven piezometers in NSB, which there is sufficient accuracy to select the piezometer well given that the purpose of selecting the piezometer well in each sub-basin is to display groundwater level. The existing wells have been analyzed using entropy theory, and according to the presented indices, the existing piezometer wells have been classified according to the shortage or excess of information, and finally, according to the obtained rankings, the selected well in terms of information transfer was investigated. The study results of the entropy method in determining the selected well in NSB are presented in Table 2.

Table 2

Results of the entropy method in NSB

PiezometerUTM-xUTM-yNiITIT
Naqadeh 534053 4090609 − 0.075 0.681 4.486 
Alagoz 527544 4094723 0.056 0.889 5.619 
Darband 522946 4093595 − 0.089 0.664 4.552 
Kohnegale 521388 4092529 − 0.012 0.827 5.393 
Cheshmeh gol 519040 4098266 0.00 0.759 4.975 
Tazegale 540708 4090982 0.122 0.946 5.740 
Agabeiglou 542852 4092602 − 0.063 0.635 4.398 
PiezometerUTM-xUTM-yNiITIT
Naqadeh 534053 4090609 − 0.075 0.681 4.486 
Alagoz 527544 4094723 0.056 0.889 5.619 
Darband 522946 4093595 − 0.089 0.664 4.552 
Kohnegale 521388 4092529 − 0.012 0.827 5.393 
Cheshmeh gol 519040 4098266 0.00 0.759 4.975 
Tazegale 540708 4090982 0.122 0.946 5.740 
Agabeiglou 542852 4092602 − 0.063 0.635 4.398 

According to Ni, referred to as the net information, it can be seen that the Tazegale piezometer in NSB has the most relationship with the other wells and the information of this well represents a better condition of the plain. In fact, the information of this well indicates the state of the plain aquifer and the values of this well can be trusted with high confidence. As shown, information transfer index (ITI) from this well is higher than other wells. After determining the selected piezometer, the rainfall deficiency and the groundwater level deficiency signatures were extracted using the desired equation, completed for 60-day consecutive duration. The results of the extraction of deficiency signature values of the studied stations in NSB for 10-, 30- and 60-day consecutive durations are presented in Figures 4 and 5. In this study, the 1- to n-day durations were investigated. Since there was no deficiency in the study area for more than 60-day duration, n was equal to 60. The 10-, 30- and 60-day durations were considered as representative of the durations because the scattering pattern distribution is the same for each of 10-, 30- and 60-day durations (Figures 4 and 5).

Figure 4

Results of rainfall deficiency signature extraction at the Naqadeh rain gauge station during the statistical period of 1973–2015 for 10-, 30- and 60-day durations.

Figure 4

Results of rainfall deficiency signature extraction at the Naqadeh rain gauge station during the statistical period of 1973–2015 for 10-, 30- and 60-day durations.

Figure 5

Results of groundwater level deficiency signature extraction in the Tazegale piezometer during the period of 2000–2015 for 10-, 30- and 60-day durations.

Figure 5

Results of groundwater level deficiency signature extraction in the Tazegale piezometer during the period of 2000–2015 for 10-, 30- and 60-day durations.

All studied signatures were completed for 60-day consecutive duration. Since 2004, the increase in rainfall deficiency signature values at the Naqadeh rain gauge station is clearly visible. Increasing the values of rainfall deficiency signature in the area results in two states: first, the reduction in rainfall in LUB during the last years and second, the change in rainfall pattern and its distribution, which increases the extreme values of rainfall and flood. The study results showed that changes in groundwater level deficiency signature have been increasing in recent years. Regarding the low groundwater level in the study area, rainfall values have a direct effect on it and there is little delay between these two parameters. Increasing groundwater level deficiency regarding the extent of the aquifer covers a large part of groundwater, making it difficult to offset this reservoir deficiency. Also, with increasing groundwater level and increasing extreme rainfalls in the study area, the soil texture is damaged and the soil structure gradually lose its resilience. This will help increase the flooding potential of the basin.

Study of marginal distribution functions

After preparing the data, different statistical distribution functions are fitted to the data, and finally, using the Kolmogorov–Smirnov goodness of fit test, best distribution functions are selected. The study results of different distributions in fitting the studied data are presented in Table 3.

Table 3

Results of fitting marginal distribution functions on rainfall deficiency and groundwater level deficiency signatures

SiteDistribution parameter
Selected distribution
10-day duration30-day duration60-day duration10-day duration30-day duration60-day duration
Naqadeh α = 142.52,
β = 9.5934,
γ = 12.124,
δ = −0.16798
ξ = −4.2459 
σ = 0.81492
μ = 5.6583
γ = −32.045 
k = 0.12004
σ = 449.3
μ = 44.319 
Wakeby LN(3P) Gen-Pareto 
Tazegale k = −1.279
σ = 10.444
μ = −1.1903 
γ = −0.01447
δ = 0.22651
λ = 90.049
ξ = 8.3257 
γ = 0.2472
δ = 0.36618
λ = 173.34
ξ = 16.01 
Gen-Pareto JSB JSB 
SiteDistribution parameter
Selected distribution
10-day duration30-day duration60-day duration10-day duration30-day duration60-day duration
Naqadeh α = 142.52,
β = 9.5934,
γ = 12.124,
δ = −0.16798
ξ = −4.2459 
σ = 0.81492
μ = 5.6583
γ = −32.045 
k = 0.12004
σ = 449.3
μ = 44.319 
Wakeby LN(3P) Gen-Pareto 
Tazegale k = −1.279
σ = 10.444
μ = −1.1903 
γ = −0.01447
δ = 0.22651
λ = 90.049
ξ = 8.3257 
γ = 0.2472
δ = 0.36618
λ = 173.34
ξ = 16.01 
Gen-Pareto JSB JSB 

In this study, two, three and four-variable statistical distribution functions were used to obtain best marginal distribution functions. One of the most important stages before modeling and applying copulas is the discussion of marginal distribution functions. The choice of marginal distribution functions influences the results of the simultaneous analysis of variables. Also, the estimated return period is affected by these marginal distribution functions. For this reason, it is important to be careful in selecting the marginal distribution functions.

Study of correlation of data using Kendall's tau test

The results of Kendall's tau test in studying the correlation of the values of groundwater level deficiency and rainfall deficiency signatures in the study area showed that with increasing consecutive durations, Kendall's tau correlation also increased. In the study area and selected signatures, the correlation coefficient is more than 0.4. The existence of correlation is required to use copula functions. In the case of poor correlation, copulas cannot be used for multivariate analysis. For this reason, correlation is a prerequisite for copula studies. According to the study results, the correlation between the values of rainfall deficiency and groundwater level deficiency signatures is confirmed and the reliability of multivariate analysis calculations can be fully confirmed. Various studies such as Bedford & Cooke (2001), Kurowicka & Cooke (2007), Aas et al. (2009), Gräler et al. (2013), Cooke et al. (2015), Bevacqua et al. (2017), Bezak et al. (2017), Zhang et al. (2018), Brunner et al. (2019), Khan et al. (2019) and Ramezani et al. (2019) also confirm and recommend Kendall's tau test to examine the correlation of variables. The results of the correlation analysis of the studied signatures are presented in Figure 6.

Figure 6

Results of the correlation analysis of groundwater level deficiency and rainfall deficiency signatures.

Figure 6

Results of the correlation analysis of groundwater level deficiency and rainfall deficiency signatures.

Figure 6 shows that Kendall's tau correlation between groundwater level deficiency and rainfall deficiency signatures in the study area for the 10-, 30- and 60-day durations is about 0.6, 0.6 and 0.8, respectively. All three are acceptable, but the correlation between the two deficiency signatures for 60-day duration is greater than the other durations. In order to select the best copula, the comparison between different copulas and empirical copulas, and the estimation of goodness of fit tests were used. The results of copula evaluation for 60-day duration are presented in Table 4.

Table 4

Results of study and comparison of different copulas in the combined analysis of groundwater level deficiency and rainfall deficiency signatures

CriteriaCopula function
ClytonAMHFGMFrankGalambosGHPlackett
AIC − 6.65 − 6.61 − 6.25 − 6.69 − 6.66 − 6.66 − 6.67 
MAE 0.12 0.15 0.21 0.09 0.10 0.10 0.11 
NSE 0.81 0.68 0.35 0.86 0.83 0.83 0.80 
RMSE 0.13 0.16 0.23 0.11 0.12 0.12 0.13 
Teta 3.01 1.00 0.47 14.37 2.44 3.14 20.00 
CriteriaCopula function
ClytonAMHFGMFrankGalambosGHPlackett
AIC − 6.65 − 6.61 − 6.25 − 6.69 − 6.66 − 6.66 − 6.67 
MAE 0.12 0.15 0.21 0.09 0.10 0.10 0.11 
NSE 0.81 0.68 0.35 0.86 0.83 0.83 0.80 
RMSE 0.13 0.16 0.23 0.11 0.12 0.12 0.13 
Teta 3.01 1.00 0.47 14.37 2.44 3.14 20.00 

According to the results presented in Table 4, Frank copula was selected to analyze both groundwater level deficiency and rainfall deficiency in NSB. Considering the selected copula and its corresponding parameter, the bivariate analysis of these two signatures for 60-day duration is presented in Figure 7. For bivariate analysis, the combined return period in the ‘AND’ state (when both parameters are less than a threshold) and in the ‘OR’ state (when each parameter is less than a threshold) were used.

Figure 7

Return period in (a) ‘AND’ and (b) ‘OR’ states of groundwater level deficiency and rainfall deficiency signatures for 60-day duration.

Figure 7

Return period in (a) ‘AND’ and (b) ‘OR’ states of groundwater level deficiency and rainfall deficiency signatures for 60-day duration.

For 60-day duration, the results showed that if rainfall deficiency was more than 120 mm in the area and groundwater level deficiency increased to at least 1.5 m, then the return period would be 1.5 years. With a return period of 1.5- to 5-year, it is possible to estimate the increase in groundwater level due to the rainfall deficiency signature. These values are further increased in the frequency analysis in the ‘OR’ state (Figure 7(b)). Finally, considering conditional probability, it can be observed that the groundwater level change in different return periods is in the range of 0.3–2.1 m.

Study of conditional density

In this section, conditional density was used to estimate the values of groundwater level deficiency signature, taking into account rainfall deficiency signature. When using conditional probability, the returns period does not provide the exact value of the given signature. Accordingly, in this study, a new method to present conditional deficiency values was developed. In this case, the following steps are considered:

A: In the first step, the values of the studied signatures are converted to copula variables and the parameter C of each copula is estimated.

B: In the next step, by estimating the copula density values (), the conditional density equation is calculated.

C: Regarding values of the studied copula and signature, can define a curve for each deficiency signature, and for the maximum value of this curve, the expected value of the deficiency signature can be estimated by considering other signatures.

In order to estimate the conditional density of the studied signature values, all of the rainfall deficiency signature values were individually entered density functions and the resulting output was studied. Its maximum indicates the highest level of groundwater level deficiency, provided that rainfall reduction in the study area for 60-day duration. Finally, the results of the conditional density analysis of groundwater level deficiency signature are presented in Figure 8. This figure is the result of the combination of the conditional density of the selected copula and its simulation, which is simplified by the diagonal section of the copula. The provided conditional graph shows the value of the dependent parameter for different values of the independent parameter.

Figure 8

Estimation of conditional values of groundwater level deficiency based on rainfall deficiency.

Figure 8

Estimation of conditional values of groundwater level deficiency based on rainfall deficiency.

Using Figure 8, groundwater level deficiency signature for a given rainfall deficiency signature can be estimated. For example, for 50, 100 and 150 mm rainfall deficiency, groundwater levels increased by 1.8, 4 and 4.8 m, respectively, for 60-day duration. It should be noted that not all groundwater level increase in the study area is related to rainfall deficiency signature and rainfall reduction. However, regarding the correlation between rainfall deficiency signature and groundwater level deficiency signature, it can be concluded that the effective factor on increasing groundwater level in the study area is rainfall. On the other hand, groundwater level in the study area is low and rainfall can affect groundwater level.

The proposed method is more accurate than the conventional method for estimating the copula parameter due to the diagonal section of Archimedean family copulas. The diagonal section of Archimedean family copulas uniquely estimates the correlation of the data. This helps reduce variates and makes it easier to understand the underlying structure without losing information (Sungur & Yang 1996). For this reason, information obtained from this section can be used to provide typical curves for each basin. Providing the conditional density method for estimating and predicting groundwater level deficiency values by reducing rainfall (bivariate state) introduced and used for the first time, adds to the accuracy of the presented graphs.

De Michele et al. (2007) and Salvadori & De Michele (2007) discussed the conditional return period, and later many studies have been conducted in this regard. However, due to the impossibility of estimating the return period of the pairwise variables in the conditional state, they later stated that the conditional return period does not provide verified information about the return period affected by the pairwise variables. For this reason, the conditional return period was not discussed in this study. An alternative method for estimating the conditional values of the pairwise variables was also attempted. The conditional density method presented in this study was able to point the estimation of groundwater level deficiency values using the diagonal section of Archimedean family copulas and considering the condition of the dependent variable possibility. The most important achievement and innovation of this study in the bivariate section is the presentation of these graphs in the sub-basin. Finally, the best line was fitted to rainfall deficiency and groundwater level deficiency data, and a mathematical equation was presented for each sub-basin.

One of the problems in the coming years in the field of artificial feeding and aquifer feeding is to change the rainfall condition to the extreme values. According to the results, it can be observed that if rainfall reduced in LUB, the groundwater level of the existing sub-basins will increase. This is due to the reduction in rainfall in the coming years as well as extreme rainfall. Increasing the extreme rainfalls causes the least rainfall to become runoff and become out of reach. This increase in the extreme rainfalls also destroys the texture and structure of the soil and gradually reduces soil permeability. Finally, it will be very difficult to make up for this deficiency. Under current conditions, extreme rainfalls have a great effect on the permeability. Therefore, accumulation of rainfall and artificial feeding projects, as well as establishment and management of runoff collection systems, can greatly offset this reservoir deficiency. Using calculated conditional density and provided graphs can easily estimate the best relationship between the two groundwater level deficiency and rainfall deficiency signatures. The study results of the relationship between the two signatures are presented in Table 5.

Table 5

Estimation of groundwater level deficiency signature based on rainfall deficiency signature

Sub-basinR2Equation
NSB 0.92 GDS = 0.8159exp(0.012 RDS) 
Sub-basinR2Equation
NSB 0.92 GDS = 0.8159exp(0.012 RDS) 

GDS, groundwater level deficiency signature (m); RDS, rainfall deficiency signature (mm).

Regarding the effective parameter involvement in the studied signatures, the results showed that the use of conditional density in joint analyses and presentation of bivariate returns period provides useful information on water resources management as well as the designs. Gräler et al. (2013) also found the use of conditional density useful for estimating multivariate return periods and acknowledged that the multivariate return period provides better information than univariate one.

Now, regarding equations presented, groundwater level deficiency in the NSB can be easily estimated. This equation can be used as groundwater level monitoring in the sub-basin. These regional equations are presented using values derived from conditional density. For this reason, the effect of the rainfall deficiency signature parameter on each sub-basin has also been observed. In general, the results of the simulation and prediction of groundwater level deficiency values for 60-day duration in the studied sub-basin using conditional density and diagonal section of existing copulas showed that bivariate copulas can be well used as simulation and prediction models. Bezak et al. (2017) for the simulation of sediment values also reported that the results obtained from copulas are acceptable and more accurate than the regression method.

CONCLUSION

Conventional univariate statistical analysis cannot provide accurate information on the multivariate nature of dependent extreme events. Therefore, the univariate statistical analysis of risk associated with these events is not possible. From a scientific perspective, it is impossible to provide a general and appropriate approach for estimating multivariate design events for a large data set and designing it because the approaches available are statistically different. So far, many practical approaches have been proposed based on the concept of univariate return period. But, the concept of multivariate return period has a different meaning. Second, for hydrological and water resources variables, the best approach is the multivariate approach where bivariate conditional distribution can be preferred. In this case, return periods are not dependent on univariate case. For this reason, it should be noted that univariate, bivariate and multivariate return periods are different. For the groundwater level deficiency variable, this is also important for the following reasons:

Groundwater level is highly dependent on ground surface feeding. The ground surface feeding also depends on the surface flow and rainfall. In addition to the abovementioned, the type and texture of the soil, and the location of the aquifers are also of great importance in the feeding of the aquifers. Also, utilizing the aquifer should not be forgotten because utilizing the aquifer is also effective in this regard. It is also expected that the presence of more variables will improve the results of process modeling. In general, it should be noted that the use of such models requires a great deal of reliable data.

Simulation using copulas is a great help in checking the accuracy of the selected copulas in the analyses because by simulating and evaluating error criteria, the accuracy of copulas can be evaluated.

In this study, using copulas and their conditional density, groundwater level deficiency signature was simulated for 60-day duration using rainfall deficiency signature. The study results of Kendall's tau correlation showed high and acceptable dependency of two deficiency signatures, which is not unexpected given the groundwater level. The results of the analysis showed that the highest reduction in rainfall in the studied sub-basin for 60-day duration was in the ‘AND’ state and in the range of 90–135 mm, which is likely to occur every year. This reduction in the rainfall in the studied sub-basin will result in a 1.2–1.7 m increase in the groundwater level. The changes in groundwater level for 60-day duration were also investigated using the conditional density of copulas considering the diagonal section of copulas. The study results of the conditional density of copulas showed that by calculating and estimating conditional density and taking into account the condition of rainfall deficiency can estimate groundwater level values with high accuracy. The results showed that considering conditional density, copula values of the pairwise signatures take into account the condition of rainfall deficiency possibility in a particular state and plot the possibility of groundwater level deficiency which the maximum groundwater level deficiency is the expected deficiency. Estimating the diagonal section of bivariate copulas can also reduce this complexity and predict and analyze groundwater level values. It is also possible to estimate the value of groundwater level in the sub-basin using the presented equations, given the value of rainfall deficiency for 60-day duration. The presented equation can be used as an alarm and monitoring system. Due to the stochastic nature of the hydrological phenomena and their dependency on other parameters, the use of copulas in the analysis and simulation of hydrological values will provide better results.

ACKNOWLEDGEMENTS

The authors thank Politecnico di Milano for providing the facilities to the first author as a visiting researcher. Also, the authors thank Iran Water Resources Management Company for providing the data.

DATA AVAILABILITY STATEMENT

Data cannot be made publicly available; readers should contact the corresponding author for details.

REFERENCES

REFERENCES
Aas
K.
,
Czado
C.
,
Frigessi
A.
&
Bakken
H.
2009
Pair-copula constructions of multiple dependence
.
Insurance: Mathematics and Economics
44
(
2
),
182
198
.
Abdi
A.
,
Hassanzadeh
Y.
,
Talatahari
S.
,
Fakheri-Fard
A.
&
Mirabbasi
R.
2017
Regional bivariate modeling of droughts using L-comoments and copulas
.
Stochastic Environmental Research and Risk Assessment
31
(
5
),
1199
1210
.
Adamson
P. T.
,
Metcalfe
A. V.
&
Parmentier
B.
1999
Bivariate extreme value distributions: an application of the Gibbs sampler to the analysis of floods
.
Water Resources Research
35
(
9
),
2825
2832
.
Bedford
T.
&
Cooke
R. M.
2001
Probability density decomposition for conditionally dependent random variables modeled by vines
.
Annals of Mathematics and Artificial Intelligence
32
(
1–4
),
245
268
.
Beersma
J. J.
&
Buishand
T. A.
2004
The joint probability of rainfall and runoff deficits in the Netherlands
.
World Water and Environmental Resources Congress 2004
,
June 27–July 1, 2004, Salt Lake City, UT. American Society of Civil Engineers (ASCE), pp. 1–10
.
Bevacqua
E.
,
Maraun
D.
,
Hobæk Haff
I.
,
Widmann
M.
&
Vrac
M.
2017
Multivariate statistical modelling of compound events via pair-copula constructions: analysis of floods in Ravenna (Italy)
.
Hydrology and Earth System Sciences
21
(
6
),
2701
2723
.
Bezak
N.
,
Rusjan
S.
,
Kramar Fijavž
M.
,
Mikoš
M.
&
Šraj
M. J. W.
2017
Estimation of suspended sediment loads using copula functions
.
Water
9
(
8
),
628
.
Brunner
M. I.
,
Furrer
R.
&
Favre
A.-C.
2019
Modeling the spatial dependence of floods using the Fisher copula
.
Hydrology and Earth System Sciences
23
(
1
),
107
124
.
Chebana
F.
&
Ouarda
T. B.
2007
Multivariate L-moment homogeneity test
.
Water Resources Research
43
(
8
),
1
18
.
Chebana
F.
&
Ouarda
T. B.
2011
Multivariate quantiles in hydrological frequency analysis
.
Environmetrics
22
(
1
),
63
78
.
Chen
L.
,
Singh
V. P.
,
Shenglian
G.
,
Hao
Z.
&
Li
T.
2011
Flood coincidence risk analysis using multivariate copula functions
.
Journal of Hydrologic Engineering
17
(
6
),
742
755
.
Cooke
R. M.
,
Kurowicka
D.
&
Wilson
K.
2015
Sampling, conditionalizing, counting, merging, searching regular vines
.
Journal of Multivariate Analysis
138
,
4
18
.
De Michele
C.
,
Salvadori
G.
,
Canossi
M.
,
Petaccia
A.
&
Rosso
R.
2005
Bivariate statistical approach to check adequacy of dam spillway
.
Journal of Hydrologic Engineering
10
(
1
),
50
57
.
De Michele
C.
,
Salvadori
G.
,
Passoni
G.
&
Vezzoli
R.
2007
A multivariate model of sea storms using copulas
.
Coastal Engineering
54
(
10
),
734
751
.
Duan
W.
,
Hanasaki
N.
,
Shiogama
H.
,
Chen
Y.
,
Zou
S.
,
Nover
D.
,
Zhou
B.
&
Wang
Y.
2019
Evaluation and future projection of Chinese precipitation extremes using large ensemble high-resolution climate simulations
.
Journal of Climate
32
,
2169
2183
.
doi:10.1175/JCLI-D-18-0465.1
.
Durrans
S.
,
Eiffe
M.
,
Thomas
W.
&
Goranflo
H.
2003
Joint seasonal/annual flood frequency analysis
.
Journal of Hydrologic Engineering
8
(
4
),
181
189
.
Favre
A. C.
,
Musy
A.
&
Morgenthaler
S.
2002
Two-site modeling of rainfall based on the Neyman-Scott process
.
Water Resources Research
38
(
12
),
43
41
.
Favre
A. C.
,
El Adlouni
S.
,
Perreault
L.
,
Thiémonge
N.
&
Bobée
B.
2004
Multivariate hydrological frequency analysis using copulas
.
Water Resources Research
40
(
1
),
1
12
.
Genest
C.
&
Favre
A. C.
2007
Everything you always wanted to know about copula modeling but were afraid to ask
.
Journal of Hydrologic Engineering
12
(
4
),
347
368
.
Genest
C.
,
Favre
A. C.
,
Béliveau
J.
&
Jacques
C.
2007
Metaelliptical copulas and their use in frequency analysis of multivariate hydrological data
.
Water Resources Research
43
(
9
),
1
12
.
Gräler
B.
,
van den Berg
M.
,
Vandenberghe
S.
,
Petroselli
A.
,
Grimaldi
S.
,
De Baets
B.
&
Verhoest
N.
2013
Multivariate return periods in hydrology: a critical and practical review focusing on synthetic design hydrograph estimation
.
Hydrology and Earth System Sciences
17
(
4
),
1281
1296
.
Harmancioglu
N. B.
,
Fistikoglu
O.
,
Ozkul
S. D.
,
Singh
V. P.
&
Alpaslan
N.
1999
Water Quality Monitoring Network Desighn
.
Kluwer
,
Boston, USA
, p.
299
.
He
H.
,
Zhou
J.
,
Yu
Q.
,
Tian
Y. Q.
&
Chen
R. F.
2007
Flood frequency and routing processes at a confluence of the middle Yellow River in China
.
River Research and Applications
23
(
4
),
407
427
.
Hult
H.
&
Lindskog
F.
2002
Multivariate extremes, aggregation and dependence in elliptical distributions
.
Advances in Applied Probability
34
(
3
),
587
608
.
Huo
A.
,
Yang
L.
,
Peng
J.
,
Cheng
Y.
&
Jiang
C.
2020
Spatial characteristics of the rainfall induced landslides in the Chinese Loess Plateau
.
Human and Ecological Risk Assessment: An International Journal
1
16
.
doi:10.1080/10807039.2020.1728517
.
Joe
H.
1997
Multivariate Models and Multivariate Dependence Concepts
.
Chapman and Hall/CRC
,
Boca Raton, FL
.
Khan
F.
,
Spöck
G.
&
Pilz
J.
2019
A novel approach for modelling pattern and spatial dependence structures between climate variables by combining mixture models with copula models
.
Kruskal
W. H.
1958
Ordinal measures of association
.
Journal of the American Statistical Association
53
(
284
),
814
861
.
Kurowicka
D.
&
Cooke
R.
2007
Sampling algorithms for generating joint uniform distributions using the vine-copula method
.
Computational Statistics & Data Analysis
51
(
6
),
2889
2906
.
Li
F.
&
Zheng
Q.
2016
Probabilistic modelling of flood events using the entropy copula
.
Advances in Water Resources
97
,
233
240
.
Luo
P.
,
Sun
Y.
,
Wang
S.
,
Wang
S.
,
Lyu
J.
,
Zhou
M.
,
Nakagami
K.
,
Takara
K.
&
Nover
D.
2020
Historical assessment and future sustainability challenges of Egyptian water resources management
.
Journal of Cleaner Production
363
,
1
11
.
doi:10.1016/j. jclepro.2020.121154
.
Ma
J.
&
Sun
Z.
2011
Mutual information is copula entropy
.
Tsinghua Science and Technology
16
(
1
),
51
54
.
Mirakbari
M.
,
Ganji
A.
&
Fallah
S.
2010
Regional bivariate frequency analysis of meteorological droughts
.
Journal of Hydrologic Engineering
15
(
12
),
985
1000
.
Mirabbasi
R.
,
Fakheri-Fard
A.
&
Dinpashoh
Y.
2012
Bivariate drought frequency analysis using the copula method
.
Theoretical and Applied Climatology
108
(
1–2
),
191
206
.
Mirabbasi
R.
,
Anagnostou
E. N.
,
Fakheri-Fard
A.
,
Dinpashoh
Y.
&
Eslamian
S.
2013
Analysis of meteorological drought in northwest Iran using the Joint Deficit Index
.
Journal of Hydrology
492
,
35
48
.
Nadarajah
S.
&
Gupta
A. K.
2006
Intensity-duration models based on bivariate gamma distributions
.
Hiroshima Mathematical Journal
36
(
3
),
387
395
.
Nelsen
R. B.
2006
An Introduction to Copulas, ser. Lecture Notes in Statistics
.
Springer
,
New York
.
Ramezani
Y.
,
Tahroudi
M. N.
&
Ahmadi
F
, .
2019
Analyzing the droughts in Iran and its eastern neighboring countries using copula functions
.
IDŐJÁRÁS
123
(
4
),
435
453
.
Salvadori
G.
&
De Michele
C.
2007
On the use of copulas in hydrology: theory and practice
.
Journal of Hydrologic Engineering
12
(
4
),
369
380
.
Salvadori
G.
,
De Michele
C.
,
Kottegoda
N. T.
&
Rosso
R.
2007
Extremes in Nature: An Approach Using Copulas
, Vol.
56
.
Springer Science & Business Media
,
Linz, Austria
.
Serinaldi
F.
,
Bonaccorso
B.
,
Cancelliere
A.
&
Grimaldi
S.
2009
Probabilistic characterization of drought properties through copulas
.
Physics and Chemistry of the Earth, Parts A/B/C
34
(
10–12
),
596
605
.
Shiau
J.
2003
Return period of bivariate distributed extreme hydrological events
.
Stochastic Environmental Research and Risk Assessment
17
(
1–2
),
42
57
.
Sklar
M.
1959
Fonctions de repartition an dimensions et leurs marges
.
Publications de l'Institut Statistique de l'Université de Paris
8
,
229
231
.
Sungur
E. A.
&
Yang
Y.
1996
Diagonal copulas of Archimedean class
.
Communications in Statistics – Theory and Methods
25
(
7
),
1659
1676
.
Tahroudi
M. N.
,
Pourreza-Bilondi
M.
,
Ramezani
Y. J. T.
&
Climatology
A
, .
2019a
Toward coupling hydrological and meteorological drought characteristics in Lake Urmia Basin, Iran
.
Theoretical and Applied Climatology
138
,
1511
1523
.
Tahroudi
M. N.
,
Siuki
A. K.
&
Ramezani
Y.
2019b
Redesigning and monitoring groundwater quality and quantity networks by using the entropy theory
.
Environmental Monitoring and Assessment
191
(
4
),
250
.
Wang
W.
,
Dong
Z.
,
Si
W.
,
Zhang
Y.
&
Xu
W.
2018
Two-dimension monthly river flow simulation using hierarchical network-copula conditional models
.
Water Resources Management
32
(
12
),
3801
3820
.
Yang
J.
,
Zhang
H.
,
Ren
C.
,
Nan
Z.
,
Wei
X.
&
Li
C. A.
2019
A cross-reconstruction method for step changed runoff series to implement frequency analysis under changing environment
.
International Journal of Environmental Research and Public Health
16
(
22
),
4345
.
Yue
S.
,
Ouarda
T. B. M. J.
&
Bobee
B.
2001
A review of bivariate gamma distributions for hydrological application
.
Journal of Hydrology
246
(
1–4
),
1
18
.
Zhang
L.
2005
Multivariate Hydrological Frequency Analysis and Risk Mapping
.
Louisiana State University and Agricultural & Mechanical College
,
Baton Rouge, LA
.
Zhang
L.
&
Singh
V.
2006
Bivariate flood frequency analysis using the copula method
.
Journal of Hydrologic Engineering
11
(
2
),
150
164
.
Zhang
L.
&
Singh
V. P.
2007
Bivariate rainfall frequency distributions using Archimedean copulas
.
Journal of Hydrology
332
(
1–2
),
93
109
.
Zhang
Q.
,
Li
J.
,
Singh
V. P.
&
Xu
C. Y.
2013
Copula-based spatio-temporal patterns of precipitation extremes in China
.
International Journal of Climatology
33
(
5
),
1140
1152
.
Zhang
D.-D.
,
Yan
D.-H.
,
Lu
F.
,
Wang
Y.-C.
&
Feng
J.
2015
Copula-based risk assessment of drought in Yunnan province, China
.
Natural Hazards
75
(
3
),
2199
2220
.