Abstract

Drought characteristics are among major inputs in the planning and management of water resources. Although numerous studies on probabilistic aspects of meteorological drought characteristics and their joint distribution functions have been reported, multivariate analysis of groundwater (GW) drought is rarely available. In this paper, while proposing a framework for statistical analysis of disturbed hydrological systems, copula-based multivariate GW drought analysis was performed in an over-drafted aquifer. For this purpose, a 1,000-year synthetic time series of naturalized GW level was produced. GW drought was monitored via the SGI index while the multivariate GW drought probability and return period were determined via copulas. Comparison between the copula and empirical GW drought probabilities using statistical goodness-of-fit tests proved sufficient accuracy of copula models in multivariate drought analysis. The results showed strong dependence among GW drought characteristics. Generally speaking, multivariate GW drought analysis incorporates major drought characteristics and provides concrete scientific basis for planning drought management strategies.

INTRODUCTION

Drought is a natural hazard affecting society, economy, environment, agriculture, and industry. Groundwater (GW) drought is a type of drought that adversely affects the GW resources and leads to lack of public water supply (Bloomfield & Marchant 2013).

There are a number of indices for GW drought monitoring. In particular, one of the indices whose variants has been used to describe GW drought is the Standardized Precipitation Index (SPI). Although SPI is preliminarily proposed to monitor meteorological drought, adoption of its concept to GW drought monitoring has increased. In some studies, GW drought monitoring with this index was questioned because of its pitfalls and limitations (e.g. Kumar et al. 2016).

Bloomfield & Marchant (2013) proposed the Standardized GW Index (SGI) for characterizing GW drought. The SGI is calculated via a non-parametric normal score convert of GW level. They calculated the SGI in their case study and compared the SGI results with those of the SPI. They found that there is a good conformation between SGI time series and formerly independently documented GW droughts. Also, relations between SGI and SPI were dependent to a range of SPI accumulations.

Li & Rodell (2014) calculated the GW Drought Index (GWI) on the basis of GW storage simulated by the Catchment Land Surface Model (CLSM) and compared the results to those computed based on the observation boreholes data. They found that the GWI derived from CLSM had a strong relationship with the GWI calculated from the observation boreholes data. The authors also found high correlation between both GWIs with the meteorological SPI-12 and SPI-24 drought indices.

Since droughts are complicated natural events, using a single monitoring index cannot provide a comprehensive estimate (Shiau et al. 2007). In several studies, multivariate parameters were used in drought analysis in lieu of using one variable. Multivariate drought analysis provides important information for water resources management and planning (Chen et al. 2011; Ma et al. 2011). For this purpose, copula functions were adopted. In recent years, the use of copula functions in multivariate hydrologic analysis has increased significantly. Numerous studies on flood analysis (Genest et al. 2007; Shiau et al. 2007), rainfall frequency analysis (Grimaldi & Serinaldi 2006), storm analysis (Salvadori & De Michele 2006), and fitting drought characteristics (Shiau 2006) have been reported. A comprehensive description of copula types and theoretical aspects of copulas may be found in a large body of the literature. For example, the popular Archimedean copulas have been adopted in several studies (e.g. Chen et al. 2011; Dehghani et al. 2019). Genest et al. (2007) described the Metaelliptical copulas, including Normal, Student t, Cauchy, and Pearson type II, and introduced goodness-of-fit tests. In another study, Wei et al. (2015) investigated multivariate dependence concepts including affiliation, association, and lower orthant dependent in copula and introduced skew normal copula. Furthermore, Wei & Kim (2017) described skew normal copula as an asymmetric multivariate copula. They used this copula to investigate the asymmetric dependence of multivariate observation. There are other types of copulas, such as Liouville copula (McNeil & Neslshova 2010), and multivariate chi-square copula (Quessy et al. 2016). Various types of copulas have been adopted and assessed in drought frequency analysis. For example, Shiau et al. (2007) used bivariate drought characteristics including drought severity and duration via copulas and calculated the bivariate probability and return period of hydrological droughts. They found that drought severity was well correlated with drought duration. In another study, Chen et al. (2011) calculated meteorological drought characteristics including drought severity, duration, interval, and the minimum SPI. Then, they used several copulas to construct multivariate joint distributions and calculated drought probability and return period. They found that the Student copula was more suitable in drought analysis. Saghafian & Mehdikhani (2014) monitored meteorological drought via the SPI index. They calculated joint distribution functions of duration, severity and peak of drought characteristics. Six copulas were tested and the best-fitted copula was selected. It was indicated that copulas were useful to describe the correspondence between drought variables. They also calculated the copula-based three-variate drought return period. Hao et al. (2017) defined the drought in Southwest China by the 3-month Standardized Precipitation Evapotranspiration Index (SPEI). They used multivariate drought characteristics including drought duration, severity and peak via copulas and calculated the multivariate probability and return period of hydrological droughts. They showed that the Student t copula was a powerful function in drought analysis. Dehghani et al. (2019) used Archimedean copulas to model the dependence structure between the hydrological drought index (SPI) and the meteorological drought index (SHDI). They found that copulas were appropriate in drought class transition forecasting.

As pointed out, many studies have reported on multivariate meteorological drought statistical analysis. However, similar studies on GW drought are lacking. Most of the work on GW drought has been performed based on observed GW level in plains with no or minimal abstraction (e.g. Bloomfield et al. 2015) and no multivariate drought analysis has been reported to our knowledge. However, in regions with high abstractions, aquifers are intensively influenced by human activities whereas the definition of drought concerns natural variability. Thus, GW drought may be performed on the aquifer time series after removing the impact of human interventions. Moreover, the duration/frequency of GW droughts are much longer/lower than those of meteorological droughts. As a result, naturalized aquifer time series generation for GW drought frequency analysis is necessary.

Dealing with GW drought, this study was apart from previous studies in that it conducted a number of inter-related tasks involving naturalization of groundwater level, generation of GW level long time series required for GW drought analysis in an over-drafted aquifer, GW drought monitoring via SGI index, selecting the best family copula based on GW drought characteristics, multivariate GW drought analysis, and comparison with univariate analysis. For this purpose, a framework was proposed for multivariate drought analysis in a disturbed hydrological system using copulas. Based on this framework, a 1,000-year artificial time series of monthly naturalized GW level was generated via the MODFLOW-ANN modeling framework as proposed by Sanginabadi et al. (2019). They monitored GW drought using SGI as determined based on the time series of naturalized GW level while drought characteristics including duration, severity and peak were determined via frequency analysis. The joint drought duration, severity, and peak were analyzed via copula functions and multivariate return period analysis of GW drought was performed.

STUDY AREA AND DATA

Qazvin Plain is one of the largest aquifers in Iran, situated in the province of Qazvin along the southern Alborz Mountains (Figure 1). Qazvin Plain water withdrawal began in 1963. The aquifer is one of unconfined single layer.

Figure 1

Geographical boundary of Qazvin Plain and its sub-regions overlaid on the ground elevation map.

Figure 1

Geographical boundary of Qazvin Plain and its sub-regions overlaid on the ground elevation map.

There are 44 precipitation gauges spread over the plain. The plain's average monthly and annual precipitation is presented in Table 1.

Table 1

Mean monthly and annual precipitation (mm) over Qazvin Plain

JanFebMarAprMayJunJulAugSepOctNovDecAnnual
28.1 29.9 37.5 38.8 36.3 9.2 2.8 1.6 1.9 12.3 26.5 31.8 256.6 
JanFebMarAprMayJunJulAugSepOctNovDecAnnual
28.1 29.9 37.5 38.8 36.3 9.2 2.8 1.6 1.9 12.3 26.5 31.8 256.6 

The mean annual precipitation over the plain is 256.6 mm while the mean annual temperature is 13.5 °C. The climate of the region is arid to cold semi-arid following the Amberege climate classification. Over time, the number of deep and semi-deep wells rose such that currently the volume of GW annual extraction through 5,000 wells is estimated at 1,800 MCM, 89% of which is consumed in agriculture. Since 1964, the overall water table drawdown is about 34 m. The drawdown is more severe in Takestan, west of Abyek, and southern parts of the plain (Figure 1).

Climate time series including precipitation, temperature, and relative humidity as well as monthly stream flow and GW level time series were collected. Hydraulic conductivity and transmissivity coefficient of the aquifer were also taken from available estimates in the region. Data on the volume of abstractions, locations of wells, and bedrock level were also gathered.

Kharrod, Abharrod and Hajiarab are the main seasonal rivers running over the plain (Figure 1). The stream flows have reduced since 1999 due to dam construction over the Abharrod river and artificial barriers upstream of Kharrod and Hajiarab rivers. The average annual streamflow close to Qazvin Plain boundary associated with Kharrod, Abharrod and Hajiarab rivers are 3.66, 0.6 and 0.6 cubic meters per second, respectively.

METHODS

The framework that we propose for multivariate drought analysis in a disturbed hydrological system is depicted in Figure 2.

Figure 2

Proposed method for multivariate drought analysis in a disturbed hydrological system using copulas.

Figure 2

Proposed method for multivariate drought analysis in a disturbed hydrological system using copulas.

The main requirements in this framework are the hydro-meteorological observed data in order to produce the naturalized state of the system. Also, a long-time series of naturalized state of the system is needed for proper drought analysis. Therefore, the generation of a synthetic time series of the groundwater level is often inevitable. In the following, drought monitoring is performed based on the naturalized condition so that the drought characteristics may be determined. As multivariate drought may be conducted via copulas, selection of the best copula family and copula type was to be performed based on the type of the hydrological system, drought properties, and goodness-of-fit criteria. The procedure is described below.

Naturalization of GW level

Drought is caused by natural variability in hydro-climatological time series and must be analyzed via natural series. Since GW level time series is severely influenced by water withdrawals, the GW level should first be naturalized for drought analysis purposes. Naturalized GW level is the level that is formed due to natural meteorological and hydrological drivers in the absence of anthropogenic effects, i.e. GW abstractions. Sanginabadi et al. (2019) proposed a naturalization framework in which MODFLOW and a neural network model were used to remove the effect of aquifer abstractions in order to produce the naturalized level.

They discretized the plain into a total of 9,028 cells, 1,000 × 1,000 m2, while the boundaries were defined following the water table contour map. Parallel boundaries with the contours were determined as time-variant specified head boundaries. Others that were perpendicular to the contours were defined as no-flow boundaries. River network characteristics were also collected and input to the model. River surface water head was determined via stage-discharge curves whereas the hydraulic conductivity of the riverbeds was determined through a calibration procedure. Monthly well abstractions over the study period were collected while unregulated abstractions were determined in calibration. Using the FAO Penman–Monteith method (Allen et al. 1989), station-based monthly evapotranspiration was calculated over the calibration and validation periods and later interpolated over the plain using the Kriging technique (Chang & Teoh 1995).

The model was calibrated/validated in unsteady status. Calibration was performed from October 2014 to September 2015 against 110 boreholes' data via a manual trial-and-error procedure. Calibrated parameters were hydraulic conductivity, specific yield, streamflow recharge, permeability coefficient, transmissivity, and unregulated abstraction. Then, the model was validated over the October 2012–September 2013 period. The calibrated MODFLOW model was run 300 times for various input combinations (precipitation, streamflow, evapotranspiration, abstraction, and initial GW level), and the groundwater level was simulated in each borehole and later interpolated in each run. Next, an artificial neural network (ANN) model was trained and tested with simulated input–output series in order to develop a surrogate model that can easily simulate the average plain GW level by reducing the MODFLOW computational load. A 1,000-year period was generated by ANN by setting the abstraction to zero; thus producing the naturalized spatially-averaged GW level.

In this study, using the results provided by Sanginabadi et al. (2019), the naturalized GW level time series over a 1,000-year period was adopted for further drought analysis as discussed below.

Precipitation analysis

Monthly precipitation time series were obtained for 44 precipitation gauges over a 50-year period from 1966 to 2016. Then, Spearman's correlation coefficient was estimated for observed precipitation records. The P value of the Spearman's correlation between two successive months is presented in Table 2 and shows that the precipitation of two successive months are independent, at 5% significant level. So, synthetic monthly precipitations were generated via assuming random independent variables.

Table 2

Spearman correlation between precipitations of two successive months

MonthOctNovDecJanFebMarAprMayJunJulAugSep
P value 0.72 0.62 0.24 0.15 0.07 0.16 0.58 0.05 0.29 0.80 0.17 0.90 
MonthOctNovDecJanFebMarAprMayJunJulAugSep
P value 0.72 0.62 0.24 0.15 0.07 0.16 0.58 0.05 0.29 0.80 0.17 0.90 

Next, for each month, 48 statistical distribution functions were fitted to the monthly precipitation time series which showed that Gamma was the best distribution. Gamma distribution has been strongly recommended in previous studies due to its flexibility in describing rainfall data with just two parameters (e.g. Rosenberg et al. 2004; Piantados et al. 2009; Saghafian & Ghobadi Hamzekhani 2015).

Then, random precipitation data were generated following the approach presented by Piantados et al. (2009), such that the marginal distribution of precipitation and monthly/yearly averages were preserved. The parameters α and β were estimated via the maximum likelihood approach for the non-zero observed value in each month. The 1,000-year monthly precipitation was generated following the calculated distribution parameters.

Streamflow analysis

For generating a 1,000-year synthetic time series of monthly streamflow, the time series of historic monthly precipitation (spatial average) and monthly streamflow were considered over the 1966–1999 period. In this period, streamflow series did not indicate a trend. Then, an ANN model was trained with historic precipitation as input and monthly corresponding streamflow volume as output. Several ANN architectures and transfer functions were tested under feed-forward standard back propagation training.

Next, synthetic monthly streamflow time series for 1,000 years was generated via the trained ANN model using the 1,000-year generated precipitation as input.

Monthly evapotranspiration is an input to the MODFLOW-ANN framework over the 1,000-year period. The monthly evapotranspiration was calculated via the FAO Penman–Monteith method (Allen et al. 1989) for the long average of weather data in each station. Then, using the Kriging method the spatial distribution of evapotranspiration was determined (Chang & Teoh 1995). Then the MODFLOW and ANN model were used to determine the natural GW level (Sanginabadi et al. 2019).

GW drought index

In this study, GW drought was monitored via the SGI index. The concept of SGI is similar to the SPI which is widely used to monitor meteorological drought (Bloomfield & Marchant 2013; Bloomfield et al. 2015). The SGI is calculated by fitting monthly probability distribution on the GW level time series in each month. Then, the fitted monthly distributions are converted to a standard normal distribution, and the monthly standardized values are combined to obtain the whole SGI time series. The SGI procedure differs from the original SPI as the latter assumes a Gamma distribution over the monthly precipitation time series. The SGI framework, as outlined above, is quite similar to the SHDI hydrological drought index procedure proposed by Dehghani et al. (2014).

A GW drought starts when the SGI index falls below zero and ends when its value moves above zero. Therefore, a GW drought event is a period in which the SGI index is continuously negative. According to the definition of SGI (Bloomfield & Marchant 2013; Bloomfield et al. 2015) and its similarity with SPI, GW drought is a condition in which the natural GW level is below a certain GW level, which can be determined in the SGI calculations. In this study, different probability distributions were tried on the monthly time series of naturalized GW level. These distributions were evaluated via the Kolmogorov–Smirnov, Anderson–Darling, and chi-square tests and the best distribution was identified for each month.

Drought characteristics

Studied drought characteristics are duration, severity, and peak as defined based on the run theory proposed by Yevjevich (1967). Drought duration (denoted by D in Figure 3) is usually defined as the time between the initiation and termination time of a drought. Drought severity (denoted by S in Figure 3) is the cumulative deficiency of a single drought event below a critical (threshold) level. Drought peak is the maximum absolute value of Xt below the critical level during a drought event (minimum SGI values). Drought characteristics are shown in Figure 3 in which X0 is a threshold below which the drought occurs (here SGI = 0).

Figure 3

Drought characteristics (Kwak et al. 2016).

Figure 3

Drought characteristics (Kwak et al. 2016).

Copula functions

A copula is a multivariate distribution defined on the unit hypercube [0,1]d when the univariate marginal distributions are uniform over the (0,1) interval. Copulas are used to combine univariate marginal distributions and produce multivariate distribution function (e.g. Grimaldi & Serinaldi 2006; Chen et al. 2011).

According to Sklar (1959), if there are n random variables, X1, X2, …, Xn with marginal distribution F1(x1), F2(x2), …, Fn(xn), then the n-dimensional joint cumulative distribution function is determined as follows: 
formula
(1)
where C: [0,1]n →[0,1] is a copula.

There are many copula types varying on tail dependence and symmetry. Generally, copulas have been divided into three families: Archimedean, Elliptical, and other types such as EFGM copulas and Extreme Value copulas (Li et al. 2015). The properties of different copula types are listed in Table 3.

Table 3

Properties of Archimedean and Elliptical (here t and Normal) copulas (Forum Joint 2010)

Copula typeTNormalArchimedean
Handling tail dependence Yes No Gumbel and Clayton have tail dependence
Frank has no tail dependence 
Symmetry Symmetric in two dimensions, but generally asymmetric in higher dimensions Standard construction is symmetric 
Copula typeTNormalArchimedean
Handling tail dependence Yes No Gumbel and Clayton have tail dependence
Frank has no tail dependence 
Symmetry Symmetric in two dimensions, but generally asymmetric in higher dimensions Standard construction is symmetric 

Among all copulas, those frequently used are Gumbel, Frank, and Clayton from the Archimedean family as well as t-copula and Normal copula from the Elliptical family (Li et al. 2015).

In drought hydrology, Archimedean and Elliptical copulas have been widely used in previous studies (e.g. Genest et al. 2007; Chen et al. 2011; Ma et al. 2011; Saghafian & Mehdikhani 2014). The characteristics of the aforementioned copulas are fully presented in previous studies (e.g. Chen et al. 2011; Ma et al. 2011; Quessy et al. 2016; Dehghani et al. 2019).

If the marginal distribution functions are continuous, the copula function C is unique. Thus, the multivariate probability density function is calculated via individual probability density and copula functions, as expressed by: 
formula
(2)
 
formula
(3)
where Fk(xk) =uk for k = 1, …, n.

Copula parameter determination

The parameters of copulas may be calculated via the maximum pseudo-likelihood method (Chen & Guo 2019). In this method, the parametric marginal distributions are substituted via empirical or rank-based marginal distributions.

The copula parameters are estimated via the maximization of pseudo log-likelihood function. Given a random sample: i = 1, …, n}, under {}, the procedure involves calculating the pseudo log-likelihood as follows: 
formula
(4)

in which , and

Then, copula parameter is calculated as: 
formula
(5)

R software packages were applied to calculate the copula parameters in this study.

Goodness of fit tests

The log likelihood function was used to select the best copula via the Akaike information criteria (AIC) and Bayesian information criteria (BIC) through the following equations (Akaike 1974): 
formula
(6)
 
formula
(7)
where L(θ) is the Likelihood function, k is the number of copula parameters, and n is the number of data. The copula with minimum AIC or BIC is the best copula (Akaike 1974).

Moreover, two goodness-of-fit tests based on empirical copula, including SnB and SnC criteria (Genest et al. 2009), were applied in order to assess the copulas. In this research, R software packages were used to calculate these criteria.

Drought return period

The univariate drought return period may be determined as follows (Shiau 2006; Chen et al. 2011): 
formula
(8)
 
formula
(9)
 
formula
(10)
where , and are the cumulative marginal distribution function of drought duration, severity and peak, respectively. E(I) is the expected drought inter-arrival time.

Since drought is a natural event, univariate frequency analysis may lead to erroneous results in drought risk calculations. So, drought should be considered as a multivariate phenomenon specified by drought duration, drought severity, and drought peak (Saghafian & Mehdikhani 2014). Furthermore, a system failure may occur when one or more characteristics of drought are less than one or more certain threshold(s). Thus, determination of drought risk and drought return period based on multivariate analysis is for drought prediction and decision making. Also, drought return period has been studied in previous studies (e.g. Shiau et al. 2007; Saghafian & Mehdikhani 2014; Kwak et al. 2016; Hao et al. 2017).

So, the multivariate return period of droughts are calculated in two shapes: return period for Dd or Ss or Pp, Tor, and return period for Dd and Ss and Pp, Tand, as below: 
formula
(11)
 
formula
(12)
Shiau (2006) and Chen et al. (2011) analyzed three-variate drought with the equations below: 
formula
(13)
According to Equations (1) and (11)–(13), tri-variate return period of drought is presented as follows (Chen et al. 2011; Saghafian & Mehdikhani 2014): 
formula
(14)
 
formula
(15)

RESULTS AND DISCUSSION

Groundwater level naturalization

MODFLOW and ANN models were applied to determine the naturalized groundwater level of Qazvin Plain. MODFLOW was calibrated for the October 2013–September 2014 period and validated for the October 2012–September 2013 period. Statistical goodness-of-fit indicators prove a good agreement between the observed and simulated GW level over calibration and validation periods. Mean error (ME), mean absolute error (MAE), root mean squared error (RMSE), and R2 were 0.11 m, 0.83 m, 1.07 m, and 0.98 for the calibration period and 0.24 m, 0.58 m, 0.73 m and 0.98 for the validated period, respectively. These values show a good agreement of the model over the calibration and validation period. A detailed description of the results of goodness-of-fit criteria was presented in Sanginabadi et al. (2019). Then, the best ANN model was selected based on the statistical indicators. The results showed that an ANN with eight neurons in the hidden layer and Logsig/Purelin for the transfer functions of input/output layers was the best model. The value of R2 and mean squared error (MSE) of the selected model were 0.98 and 0.001 m in the training state, 0.97 and 0.008 m in the validated state, and 0.99, 0.002 m in the test state, respectively. These values indicate that the model has good accuracy in GW level naturalization.

GW level time series naturalization over past 50 years

Using the MODFLOW-ANN model, the naturalized GW level was determined during 1966–2015. The results of the naturalized GW level (spatially averaged) and observed GW level are presented in Figure 4.

Figure 4

Observed and naturalized GW level between 1966 and 2015.

Figure 4

Observed and naturalized GW level between 1966 and 2015.

Figure 4 shows that a significant drop in the GW level occurred during the 50-year-long period due to the GW abstraction. The value of GW level drop is about 45 m during these years, i.e. an annual average of 0.9 m. More details of this issue were presented in Sanginabadi et al. (2019).

GW level time series naturalization and generation

In this study, a 1,000-year time series of naturalized GW level was generated. For this purpose, a 1,000-year time series of input data to the MODFLOW-ANN model is required. First, such data was generated for the precipitation via the procedure mentioned previously.

Results showed that the Gamma distribution was the best distribution to fit the monthly precipitation data that satisfies the Kolmogorov–Smirnov and chi-squared tests at p = 0.05 significance level. Then, based on the fitted Gamma distributions for each month, 1,000 random time series for each calendar month were generated. As Table 4 shows, the statistical characteristics of observed precipitation data are preserved in the generation process.

Table 4

Comparison between the statistical properties of observed and simulated precipitation

MonthMean
Standard deviation
ObservedSimulatedObservedSimulated
September 11.30 10.14 14.30 13.30 
October 32.20 33.46 27.90 29.71 
November 34.68 34.70 24.35 24.44 
December 28.64 29.52 17.98 19.24 
January 31.72 31.89 12.83 12.94 
February 36.32 35.36 23.14 23.24 
March 40.27 39.92 23.78 24.80 
April 39.37 38.86 25.69 25.73 
May 11.89 11.92 12.39 12.10 
June 3.87 3.77 4.97 4.64 
July 2.23 2.31 2.71 2.70 
August 3.16 3.21 4.93 5.17 
MonthMean
Standard deviation
ObservedSimulatedObservedSimulated
September 11.30 10.14 14.30 13.30 
October 32.20 33.46 27.90 29.71 
November 34.68 34.70 24.35 24.44 
December 28.64 29.52 17.98 19.24 
January 31.72 31.89 12.83 12.94 
February 36.32 35.36 23.14 23.24 
March 40.27 39.92 23.78 24.80 
April 39.37 38.86 25.69 25.73 
May 11.89 11.92 12.39 12.10 
June 3.87 3.77 4.97 4.64 
July 2.23 2.31 2.71 2.70 
August 3.16 3.21 4.93 5.17 

The results also showed that the trained precipitation-streamflow ANN model with eight neurons in the hidden layer performed best among the tested ANN models when the Purelin and Tansige represented input and output layers transfer functions, respectively. Statistical goodness-of-fit for the ANN model are presented in Table 5.

Table 5

Performance of selected ANN streamflow model

Statistical indicatorTraining phaseValidation phaseTest phase
R2 0.82 0.90 0.83 
MSE (m3/s) 0.00228 0.00825 0.00909 
Statistical indicatorTraining phaseValidation phaseTest phase
R2 0.82 0.90 0.83 
MSE (m3/s) 0.00228 0.00825 0.00909 

Another ANN model trained on the MODFLOW model was used to generate the naturalized GW level for 1,000 years. The results showed that, under natural conditions, the difference between maximum and minimum of naturalized GW level in Qazvin Plain was about 17 m over 1,000 years.

GW drought index

Regarding the SGI, it was found that the Burr distribution was the best fit over the generated naturalized GW level in all 12 months. Subsequently, the monthly fitted distributions were transformed to standard normal distributions and the estimated standardized SGI values in each month were combined to construct the whole SGI time series. Based on SGI time series, drought periods were determined. The characteristics of drought events identified by the SGI were determined over the 1,000 period (Table 6).

Table 6

Characteristics of drought events in Qazvin Plain over 1,000 years

Number of droughtsDrought duration (month)
Drought severity
Drought peak
MinimumAverageMaximumMinimumAverageMaximumMinimumAverageMaximum
222 27 233 0.02 21.7 262 0.01 0.79 3.41 
Number of droughtsDrought duration (month)
Drought severity
Drought peak
MinimumAverageMaximumMinimumAverageMaximumMinimumAverageMaximum
222 27 233 0.02 21.7 262 0.01 0.79 3.41 

According to Table 6, GW drought duration varied from two months to 20 years in Qazvin Plain while drought severity varied widely.

The average time between two consecutive GW droughts was 26 months. Comparing this time with the average drought duration in Table 6 indicates that the aquifer was in dry condition for about 50% of the time.

In this study, 60 distributions were fitted on drought characteristics and the best distribution was selected via the Kolmogorov–Smirnov, Anderson–Darling, and chi-square tests. The results indicated that Lognormal (3p), Burr and Exponential distributions fitted the drought duration, severity and peak, respectively, by rejecting the Kolmogorov–Smirnov, Anderson–Darling, and chi-square tests at p = 0.05 significance level. The marginal distribution parameters were calculated via maximum likelihood method and the PDFs and P-P plots for duration, severity and peak drought are shown in Figures 57, respectively. The figures indicate that the selected distributions fit the duration, severity and peak time series well. The parameters of the marginal distribution for drought duration are presented in Table 7.

Table 7

Parameters of marginal distributions

Drought characteristicSuitable marginal distributionParameters of the marginal distribution
Duration Lognormal (3p) σ = 1.86, μ = 2.02, γ = 1.87 
Severity Burr k = 0.78, α = 0.79, β = 1.95, γ = 0.017 
Peak Exponential λ = 1.27, γ = 0 
Drought characteristicSuitable marginal distributionParameters of the marginal distribution
Duration Lognormal (3p) σ = 1.86, μ = 2.02, γ = 1.87 
Severity Burr k = 0.78, α = 0.79, β = 1.95, γ = 0.017 
Peak Exponential λ = 1.27, γ = 0 
Figure 5

Probability density function and P-P plot of drought duration in Qazvin Plain.

Figure 5

Probability density function and P-P plot of drought duration in Qazvin Plain.

Figure 6

Probability density function and P-P plot of drought severity in Qazvin Plain.

Figure 6

Probability density function and P-P plot of drought severity in Qazvin Plain.

Figure 7

Probability density function and P-P plot of drought peak in Qazvin Plain.

Figure 7

Probability density function and P-P plot of drought peak in Qazvin Plain.

Figures 57 indicate that GW droughts with smaller duration, severity and peak have more probability. Also, the P-P plots in these figures show that selected distributions fit the drought characteristics well.

Copula functions

To select the appropriate copula and determine the probability of GW drought, it is necessary to evaluate the dependence between drought duration, severity and peak. According to Table 8, GW drought characteristics are strongly dependent.

Table 8

Tau Kendall and Pearson coefficient between groundwater drought characteristics

Duration-SeveritySeverity-PeakDuration-peak
Pearson coefficient 0.96 0.87 0.90 
Tau Kendall 0.86 0.86 0.75 
Duration-SeveritySeverity-PeakDuration-peak
Pearson coefficient 0.96 0.87 0.90 
Tau Kendall 0.86 0.86 0.75 

Moreover, the form of dependence between drought characteristics may be studied using the pairwise scatterplots of the ranks (Genest et al. 2007; Chen et al. 2011) as shown in Figure 8.

Figure 8

The pair wise scatterplots of the standardized ranks of drought variables.

Figure 8

The pair wise scatterplots of the standardized ranks of drought variables.

Figure 8 demonstrates positive dependence in all drought characteristics pairs including duration-severity, duration-peak and severity-peak. Duration-severity and duration-peak have upper tail dependence while severity-peak has lower tail dependence. Besides, groundwater level is a continuous variable and its structure is generally asymmetric. So, for groundwater drought analysis, it is necessary to use copula families that have the ability of positive dependence, low to high tail dependence, and asymmetry structure. Accordingly, this study took advantage of Archimedean and Elliptical families of copula to cover these specifications (Table 3). Five copulas of these families, including Frank, Gumbel, Clayton, Normal, and Student t, were applied.

The values of parameter θ, AIC and BIC for duration-severity, severity-peak and duration-peak pairs corresponding to five tested copulas are presented in Table 9, showing that for duration-peak, duration-severity, severity-peak, and duration-severity-peak the values of the AIC and BIC criteria of Frank, Student t, Normal, and Normal copulas are less than those of other copulas, respectively. Thus, these copulas are the best based on the AIC and BIC criteria. Furthermore, the goodness-of-fit test was performed for each pair of groundwater drought characteristics. The number of bootstrap samples (primary level) used for the computation of p-values was N = 1,000 in all cases. The results are shown in Table 10. Tables 9 and 10 imply that Student t, Student t, Normal, and Normal copulas were the best for duration-severity, duration-peak, severity-peak, and duration-severity-peak analysis, respectively. Hence, these copulas were selected and used to determine the dependence structure among drought characteristics.

Table 9

Copula parameters, AIC and BIC criteria

CopulaDuration-Severity
Duration-Peak
Severity-Peak
Duration- Severity-Peak
parameterAICBICParameterAICBICparameterAICBICparameterAICBIC
Frank 23.56 –500 –496 13.16 –335 –332 25.68 –264 –260 18.02 –452 –449 
Clayton 4.14 –238 –234 2.56 –159 –155 7.24 –208 –204 3.80 –165 –161 
Gumbel 5.77 –477 –473 3.3 –321 –317 5.65 –261 –258 4.02 –354 –350 
Student t 0.96,4 –503 –500 0.89,3 –1,142 –1,141 .96,5 –630 –628 0.95,.89,0.96,5 –1,119 –1,120 
normal 0.95 –440 –436 0.89 –286 –283 0.97 –661 –656 .92,.95,0.8 –1,166 –1,166 
CopulaDuration-Severity
Duration-Peak
Severity-Peak
Duration- Severity-Peak
parameterAICBICParameterAICBICparameterAICBICparameterAICBIC
Frank 23.56 –500 –496 13.16 –335 –332 25.68 –264 –260 18.02 –452 –449 
Clayton 4.14 –238 –234 2.56 –159 –155 7.24 –208 –204 3.80 –165 –161 
Gumbel 5.77 –477 –473 3.3 –321 –317 5.65 –261 –258 4.02 –354 –350 
Student t 0.96,4 –503 –500 0.89,3 –1,142 –1,141 .96,5 –630 –628 0.95,.89,0.96,5 –1,119 –1,120 
normal 0.95 –440 –436 0.89 –286 –283 0.97 –661 –656 .92,.95,0.8 –1,166 –1,166 
Table 10

Goodness-of-fit tests in copula evaluation

CopulaDuration-Severity
Duration-Peak
Severity-Peak
Severity-Duration-Peak
Pvalu_SnBPvalu_SnCPvalu_SnBPvalu_SnCPvalu_SnBPvalu_SnCPvalu_SnBPvalu_SnC
Frank 0.70 0.70 0.330 0.447 0.640 0.420 0.54 0.32 
Clayton 0.001 0.005 0.002 0.005 0.005 0.003 0.005 0.005 
Gumbel 0.660 0.620 0.224 0.140 0.002 0.001 0.004 0.004 
Student t 0.760 0.700 0.480 0.490 0.60 0.460 0.62 0.50 
Normal 0.460 0.154 0.350 0.290 0.73 0.460 0.68 0.61 
CopulaDuration-Severity
Duration-Peak
Severity-Peak
Severity-Duration-Peak
Pvalu_SnBPvalu_SnCPvalu_SnBPvalu_SnCPvalu_SnBPvalu_SnCPvalu_SnBPvalu_SnC
Frank 0.70 0.70 0.330 0.447 0.640 0.420 0.54 0.32 
Clayton 0.001 0.005 0.002 0.005 0.005 0.003 0.005 0.005 
Gumbel 0.660 0.620 0.224 0.140 0.002 0.001 0.004 0.004 
Student t 0.760 0.700 0.480 0.490 0.60 0.460 0.62 0.50 
Normal 0.460 0.154 0.350 0.290 0.73 0.460 0.68 0.61 

In the next step, the joint probability distribution of GW drought characteristics was estimated. For example, Figure 9 shows the duration-severity cumulative and density probability distribution.

Figure 9

Groundwater drought duration-severity density and cumulative probability distribution. (a): Probability density function. (b): Cumulative probability function.

Figure 9

Groundwater drought duration-severity density and cumulative probability distribution. (a): Probability density function. (b): Cumulative probability function.

Figure 9(a) indicates that the probability of GW drought with duration less than about 50 months under all severity conditions is fairly high in comparison with other duration-severity combinations. Furthermore, the distance between the contours of cumulative probability in Figure 9(b) increases with increasing drought severity and duration.

To evaluate the proposed model, three-variate empirical probabilities of historic GW drought events (P(Dd were compared with the derived copula probabilities. The empirical probability is the ratio of the number of the outcomes (droughts with to the total number of trials (total droughts). Derived copula probabilities were calculated using selected copulas and Equation (13). In order to measure the copula accuracy, R2, root mean squared error (RMSE), Nash–Sutcliffe efficiency (NSE), and mean absolute error (MAE) criteria between the empirical and copula-derived probabilities were computed as presented in Table 11. Values of the criteria in Table 11 indicate sufficient accuracy of the copula model.

Table 11

Comparison between the empirical and copula-derived probabilities

CriteriaR2Root mean squared errorNash–Sutcliffe efficiencyMean absolute error
Value 0.96 0.15 0.72 0.13 
CriteriaR2Root mean squared errorNash–Sutcliffe efficiencyMean absolute error
Value 0.96 0.15 0.72 0.13 

Furthermore, the duration, peak, severity, and return period relationships were determined. GW drought duration, severity and peak corresponding to various univariate and multivariate return periods are listed in Table 12.

Table 12

Univariate and multivariate return period of groundwater drought

Return period (year)Drought duration (month)Drought severityDrought peakMultivariate return period (Tor)Multivariate return period (Tand)
1.27 0.015 0.007 4.4 4.7 
2.7 0.32 0.1 4.7 5.2 
10 11.8 4.51 0.64 8.0 12.2 
20 33.7 17.9 1.19 14.6 24.5 
50 96.3 64.9 1.92 33.3 60.4 
100 184.5 142.2 2.46 62.9 119.1 
200 326.9 281.7 3.01 120.3 235.5 
500 638.8 626.2 3.74 286.8 578.9 
Return period (year)Drought duration (month)Drought severityDrought peakMultivariate return period (Tor)Multivariate return period (Tand)
1.27 0.015 0.007 4.4 4.7 
2.7 0.32 0.1 4.7 5.2 
10 11.8 4.51 0.64 8.0 12.2 
20 33.7 17.9 1.19 14.6 24.5 
50 96.3 64.9 1.92 33.3 60.4 
100 184.5 142.2 2.46 62.9 119.1 
200 326.9 281.7 3.01 120.3 235.5 
500 638.8 626.2 3.74 286.8 578.9 

Comparison among univariate and multivariate drought return periods in Table 12 indicates that all univariate return periods fall between the corresponding multivariate return periods (Tor and Tand). This outcome was expected according to the definition of multivariate return period given in Equations (11) and (12).

Multivariate joint distribution of GW drought involves all major drought aspects. Since for a given drought duration, severity and peak, the difference between the univariate return period and the corresponding multivariate drought is significant in Table 12, it is expected that univariate GW drought analysis leads to error in drought quantification and thus in mitigation measures. This issue may relate to the degree of correlation between the GW drought characteristics. Table 12 indicates that the multivariate joint GW drought return period grows larger as the duration, severity or peak increase. For the same value of duration, severity and peak, Tand (Equation (15)) is greater than Tor (Equation (14)); the difference significantly grows with an increase in drought duration, severity and peak.

In Figure 10 the relationship between duration, severity and return period of GW drought using three-variable copula function is shown. The uniform distance between contours in Figure 10 indicates regular changes in the drought return period with changes in drought severity and duration.

Figure 10

Duration-severity-return period (year) of groundwater drought (Tand).

Figure 10

Duration-severity-return period (year) of groundwater drought (Tand).

Also, in Figures 1113, the relationship between severity, peak, and return period of GW drought corresponding to droughts with two, five and 10-year duration via three-variate copula are shown.

Figure 11

Severity-peak-return period (year) of groundwater drought with two-year duration (Tand).

Figure 11

Severity-peak-return period (year) of groundwater drought with two-year duration (Tand).

Figure 12

Severity-peak-return period (year) of groundwater drought with five-year duration (Tand).

Figure 12

Severity-peak-return period (year) of groundwater drought with five-year duration (Tand).

Figure 13

Severity-peak-return period (year) of groundwater drought with 10-year duration (Tand).

Figure 13

Severity-peak-return period (year) of groundwater drought with 10-year duration (Tand).

According to Figures 1013, the joint return period grows larger as the duration or severity increases. Thus, GW droughts with higher severity or longer duration have greater return periods. Moreover, comparison among Figures 1113 indicates that the contour curvature of the severity-peak-return period are similar for different GW drought durations.

According to the naturalized GW SGI, the aquifer experienced four periods of GW drought over the 40-year period. The characteristics of such droughts including drought duration, severity, and peak are presented in Table 13, in which listed return periods (Tor and Tand) were determined via Equations (14) and (15). According to Table 13, the return period of historic GW droughts (Tand) vary from 12.1 to 66.3 years. The longest drought spell occurred at the end of the time series with the duration of 97 months, severity of 60.8 and a return period of 66.3 years. Also, the shortest drought has a duration of 14 months, severity of 2.2 and return period of 12.1 years. Moreover, the cumulative duration of all droughts is 196 months which is about 40% of the total study period (40 years from 1976 to 2016). This observation implies that the study region is being severely impacted by prolonged GW droughts.

Table 13

Groundwater drought spells over the 40-year period

Drought spellDuration (month)SeverityPeakTorTand
May 1977–January 1980 33 3.97 0.37 6.9 22.2 
November 1990–December 1991 14 2.2 0.31 6.3 12.1 
May 1999–August 2003 52 24.3 0.73 10.8 37.4 
March 2008–March 2016 97 60.8 1.21 19.2 66.3 
Drought spellDuration (month)SeverityPeakTorTand
May 1977–January 1980 33 3.97 0.37 6.9 22.2 
November 1990–December 1991 14 2.2 0.31 6.3 12.1 
May 1999–August 2003 52 24.3 0.73 10.8 37.4 
March 2008–March 2016 97 60.8 1.21 19.2 66.3 

The results of GW drought analysis conducted in this study indicated that the multivariate approach is more comprehensive compared with the univariate analysis. Accordingly, the multivariate estimation of probability and return period provide helpful tools in drought prediction and decision making.

CONCLUSIONS

GW drought is a natural event that has a severe impact on the society, economy, environment, and industry. The impact is much more severe in semi-arid and arid regions where most water demand is met by GW, a major supply source that is under-studied in terms of drought quantification. Furthermore, droughts have various characteristics, the main of which are duration, severity, and peak. As a result, GW drought analysis must incorporate multivariate drought aspects.

In this research, we proposed a framework for multivariate GW drought analysis in a disturbed hydrological system using copulas, case study Qazvin Plain, Iran. For this purpose, a 1,000-year long naturalized GW level time series was generated in order to monitor GW drought via an SGI index. Furthermore, duration, severity and peak drought characteristics were studied via frequency analysis. Univariate and multivariate probability and return period of historical GW drought were determined and evaluated via copula functions.

The main conclusions of this study are as follows:

  • 1.

    A substantial drop has occurred in the GW level of Ghazvin plain during the past 50 years due to substantial over-draft.

  • 2.

    Historic drought duration varied from two months to 20 years on the basis of a 1,000-year synthetic time series of naturalized GW level. On the other hand, we analyze drought duration based on the 1,000-year synthetically generated time series of naturalized GW level. The shortest GW drought had a duration of two months and the longest GW occurred with a duration of 20 years over the 1,000-year period. Such an outcome proves that the reaction of the aquifer to precipitation and streamflow comes with significant time lag, as related to the hydrodynamic character of the aquifer, high bedrock depth, and the existing GW level. The depth of the GW table is high and the influence of precipitation and streamflow on GW level occurs with a long delay. Moreover, the duration-severity-probability curve of drought indicates that the occurrence probability of droughts with duration less than about 50 months under all conditions of severity is fairly high in comparison with other duration-severity combinations. The results of multivariate GW drought analysis revealed that prolonged GW drought is to be expected.

  • 3.

    Based on the GW drought characteristics, Elliptical family copulas were found suitable copula types in multivariate GW drought analysis.

  • 4.

    Strong dependence was found among GW drought characteristics. Based on AIC, BIC and SnB and SnC criteria, Student t, Student t and Normal were the best copula models to characterize the dependence structure for the pairs of duration-peak, duration-severity and severity-peak, respectively.

  • 5.

    Comparison between the copula and empirical GW drought probabilities via statistical goodness-of-fit criteria indicated sufficient accuracy of the fitted copula models. Thus, copula was found as a powerful tool for multivariate GW drought analysis.

  • 6.

    For a given drought duration, severity and peak, the difference between the univariate return period and the corresponding multivariate drought is significant. It is expected that univariate GW drought analysis introduces errors in drought quantification and thus in mitigation measures.

  • 7.

    Multivariate analysis of GW drought based on the proposed framework evaluates GW drought more comprehensively as compared with the univariate analysis. The multivariate analysis provides essential information in drought prediction and decision making.

REFERENCES

REFERENCES
Akaike
H.
1974
A new look at the statistical model identification
.
IEEE Trans. Autom. Control
19
(
6
),
716
723
.
Allen
R. G.
,
Jansen
M. E.
,
Wright
J. L.
&
Burman
R. D.
1989
Operational estimates of reference evapotranspiration
.
Agron. J.
81
,
650
662
.
Bloomfield
J. P.
&
Marchant
B. P.
2013
Analysis of groundwater drought building on the standardised precipitation index approach
.
Hydrol. Earth Syst. Sci.
17
,
4769
4787
.
Bloomfield
J. P.
,
Marchant
B. P.
,
Bricker
S. H.
&
Morgan
R. B.
2015
Regional analysis of groundwater droughts using hydrograph classification
.
Hydrol. Earth Syst. Sci.
19
(
10
),
4327
4344
.
Chang
T. J.
&
Teoh
C. B.
1995
Use of the Kriging method for studying characteristics of groundwater droughts
.
J. Am. Water Resour. Assoc.
257
,
1001
1007
.
Chen
L.
&
Guo
S.
2019
Copulas and its Application in Hydrology and Water Resources
.
Springer Water
,
Singapore
.
Chen
L.
,
Singh
V. P.
&
Gao
S.
2011
Drought analysis based on copulas
. In
2011 Symposium on Data- Driven Approaches to Droughts
.
Purdue University
,
Lafayette
,
USA
.
Dehghani
M.
,
Saghafian
B.
,
Nasiri Saleh
F.
,
Farokhnia
A.
&
Noori
R.
2014
Uncertainty analysis of streamflow drought forecast using artificial neural networks and Monte-Carlo simulation
.
International Journal of Climatology
34
(
4
),
1169
1180
.
Forum Joint
2010
Developments in Modelling Risk Aggregation
.
Bank for International Settlements
,
Basel Francies, WB
,
Switzerland
.
Genest
C.
,
Favre
A. C.
,
Beliveau
J.
&
Jscques
C.
2007
Metaelliptical copulas and their use in frequency analysis of multivariate hydrological data
.
Water Resour. Res.
43
(
9
),
W0901
.
Genest
C.
,
Remillard
B.
&
Beaudoinc
D.
2009
Goodness-of-fit tests for copulas: a review and a power study
.
Insurance Math. Econom.
44
,
199
213
.
Grimaldi
S.
&
Serinaldi
F.
2006
Design hyetographs analysis with 3-copula function
.
Hydrol. Sci. J.
51
(
2
),
223
228
.
Hao
C.
,
Zhang
J.
&
Yao
F.
2017
Multivariate drought frequency estimation using copula method in southwest China
.
J. Theor. Appl. Climatol.
127
,
977
991
.
Kumar
R.
,
Musuuza
J. L.
,
Van Loon
A. F.
,
Teuling
A. J.
,
Barthel
R.
,
Ten Broek
J.
,
Mai
J.
,
Samaniego
L.
&
Attinger
S.
2016
Multiscale evaluation of the standardized precipitation index as a groundwater drought indicator
.
Hydrology and Earth System Sciences
20
,
1117
1131
.
Li
B.
&
Rodell
M.
2014
Evaluation of a model-based groundwater drought indicator in the conterminous U.S
.
J. Hydrol.
526
,
78
88
.
doi:10.1016/j.jhydrol.2014.09.027
.
Li
J.
,
Zhu
X.
,
Lee
C. F.
,
Wu
D.
,
Feng
J.
&
Shi
Y.
2015
On the aggregation of credit, market and operational risks
.
Rev. Quant. Finan. Acc.
44
,
161
189
.
Ma
M.
,
Song
S.
,
Ren
L.
,
Jiang
S.
&
Song
J.
2011
Multivariate drought characteristics using trivariate Gausian and Student t copulas
.
Hydrol. Process.
27
(
8
),
1175
1190
.
McNeil
A. J.
&
Neslshova
J.
2010
From archimedean to liouville copulas
.
J. Multivar. Anal.
101
(
8
),
1772
1790
.
Piantados
J.
,
Boland
J.
&
Howlett
P. H.
2009
Generating synthetic rainfall on various timescales daily, monthly and yearly
.
J. Environ. Models Assess.
14
,
431
438
.
Quessy
F. J.
,
Rivest
L. P.
&
Toupin
M. H.
2016
On the family of multivariate chi-square copulas
.
J. Multivar. Anal.
152
,
20
40
.
Rosenberg
K.
,
Boland
J. W.
&
Howlett
P. G.
2004
Simulation of monthly rainfall totals
.
Anzim J.
46
(
E
),
E85
E104
.
Saghafian
B.
&
Mehdikhani
H.
2014
Drought characterization using a new copula- based trivariate approach
.
Nat. Hazards
8
,
1391
1407
.
Saghafian
B.
&
Ghobadi Hamzekhani
F.
2015
Hydrological drought early warning based on rainfall threshold
.
Natural Hazards
79
(
2
),
815
832
.
Salvadori
G.
&
De Michele
C.
2006
Statistical characterization of temporal structure of storms
.
Adv. Water Resour.
29
,
827
842
.
Sanginabadi
H.
,
Saghafian
B.
&
Delavar
M.
2019
Coupled groundwater drought-water scarcity index for intensively over-drafted aquifers
.
J. Hydrol. Eng.
24
(
4
),
04019003
.
Shiau
J. T.
,
Feng
S.
&
Nadarajah
S.
2007
Assessment of hydrological for the Yellow River, China, using copulas
.
Hydrol. Process.
21
(
16
),
2157
2163
.
Sklar
A.
1959
N-dimensional distribution functions and their margins
.
Publ. Inst. Stat. Univ. Paris
8
,
229
231
.
(in French)
.
Wei
Z.
,
Wang
T.
&
Nguyen
P. A.
2015
Multivariate dependence concepts through copulas
.
Int. J. Approximate Reasoning
65
,
24
33
.
Yevjevich
V.
1967
An objective approach to definitions and investigations of continental hydrologic droughts. Hydrology Papers 23. Colorado State University, Fort Collins, CO
.