## 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

## 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.

Sub-basin . | Mean groundwater level (m) . | Average annual rainfall (mm) . | Piezometric well . | Rain gauge . |
---|---|---|---|---|

Naqadeh | 5.22 | 328.1 | Tazegale | Naqadeh |

Sub-basin . | Mean groundwater level (m) . | Average annual rainfall (mm) . | Piezometric well . | Rain gauge . |
---|---|---|---|---|

Naqadeh | 5.22 | 328.1 | Tazegale | Naqadeh |

### 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

*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: 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

*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

*x*

_{1}, …,

*x*follow the arbitrary marginal distribution functions

_{n}*F*

_{1}(

*x*

_{1}), …,

*F*(

_{n}*x*) then there is copula

_{n}*C*that combines these marginal distribution functions to the joint distribution function

*F*(

*x*

_{1}, …,

*x*) in the following:

_{n}If the marginal distribution functions of *F _{i}*(

*x*) are continuous, copula

_{i}*C*is unique. Conversely, if

*C*is a

*k*-dimension copula,

*F*is a

*n*-dimension distribution function and

*F*(

_{i}*x*), …,

_{i}*F*(

_{n}*x*) 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.

_{n}### Conditional copulas

Similarly, an equivalent formula for the conditional distribution function for *Y* with respect to *X _{x}* 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

*X*

_{1}and

*X*

_{2}are two random vectors, Kendall's tau correlation coefficient can be written as follows (Kruskal 1958):

### Goodness of fit tests

*P*and

_{ei}*P*are empirical and theoretical probabilities, respectively;

_{i}*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.

*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

*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: 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:

*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 I

^{2}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:

*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 : I

^{2}→ I is as follows:

*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: where . It can then be shown that for each , all the family copulas {} are defined as follows: 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.

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.

Piezometer . | UTM-x
. | UTM-y
. | Ni . | ITI . | T
. |
---|---|---|---|---|---|

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 |

Piezometer . | UTM-x
. | UTM-y
. | Ni . | ITI . | T
. |
---|---|---|---|---|---|

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).

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.

Site . | Distribution parameter . | Selected distribution . | ||||
---|---|---|---|---|---|---|

10-day duration . | 30-day duration . | 60-day duration . | 10-day duration . | 30-day duration . | 60-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 |

Site . | Distribution parameter . | Selected distribution . | ||||
---|---|---|---|---|---|---|

10-day duration . | 30-day duration . | 60-day duration . | 10-day duration . | 30-day duration . | 60-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 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.

Criteria . | Copula function . | ||||||
---|---|---|---|---|---|---|---|

Clyton . | AMH . | FGM . | Frank . | Galambos . | GH . | Plackett . | |

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 |

Criteria . | Copula function . | ||||||
---|---|---|---|---|---|---|---|

Clyton . | AMH . | FGM . | Frank . | Galambos . | GH . | Plackett . | |

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.

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.

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.

Sub-basin . | R^{2}
. | Equation . |
---|---|---|

NSB | 0.92 | GDS = 0.8159exp(0.012 RDS) |

Sub-basin . | R^{2}
. | Equation . |
---|---|---|

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.