## 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 Standardized GW Index (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.

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

Jan . | Feb . | Mar . | Apr . | May . | Jun . | Jul . | Aug . | Sep . | Oct . | Nov . | Dec . | Annual . |
---|---|---|---|---|---|---|---|---|---|---|---|---|

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 |

Jan . | Feb . | Mar . | Apr . | May . | Jun . | Jul . | Aug . | Sep . | Oct . | Nov . | Dec . | Annual . |
---|---|---|---|---|---|---|---|---|---|---|---|---|

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.

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 m^{2}, 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.

Month . | Oct . | Nov . | Dec . | Jan . | Feb . | Mar . | Apr . | May . | Jun . | Jul . | Aug . | Sep . |
---|---|---|---|---|---|---|---|---|---|---|---|---|

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 |

Month . | Oct . | Nov . | Dec . | Jan . | Feb . | Mar . | Apr . | May . | Jun . | Jul . | Aug . | Sep . |
---|---|---|---|---|---|---|---|---|---|---|---|---|

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 *X _{t}* below the critical level during a drought event (minimum SGI values). Drought characteristics are shown in Figure 3 in which

*X*is a threshold below which the drought occurs (here SGI = 0).

_{0}### 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).

*n*random variables,

*X*

_{1}

*, X*

_{2}

*, …, X*

_{n}with marginal distribution

*F*

_{1}(

*x*

_{1})

*, F*

_{2}(

*x*

_{2})

*, …, F*

_{n}(

*x*

_{n})

*,*then the

*n*-dimensional joint cumulative distribution function is determined as follows: where

*C:*[0,1]

*→[0,1] is a copula.*

^{n}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.

Copula type . | T . | Normal . | Archimedean . |
---|---|---|---|

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 type . | T . | Normal . | Archimedean . |
---|---|---|---|

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

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

in which , and

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

### Goodness of fit tests

*θ*) 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

*et al.*2011): 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).

*et al.*2011; Saghafian & Mehdikhani 2014):

## 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 R^{2} 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 R^{2} 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 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.

Month . | Mean . | Standard deviation . | ||
---|---|---|---|---|

Observed . | Simulated . | Observed . | Simulated . | |

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 |

Month . | Mean . | Standard deviation . | ||
---|---|---|---|---|

Observed . | Simulated . | Observed . | Simulated . | |

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.

Statistical indicator . | Training phase . | Validation phase . | Test phase . |
---|---|---|---|

R^{2} | 0.82 | 0.90 | 0.83 |

MSE (m^{3}/s) | 0.00228 | 0.00825 | 0.00909 |

Statistical indicator . | Training phase . | Validation phase . | Test phase . |
---|---|---|---|

R^{2} | 0.82 | 0.90 | 0.83 |

MSE (m^{3}/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).

Number of droughts . | Drought duration (month) . | Drought severity . | Drought peak . | ||||||
---|---|---|---|---|---|---|---|---|---|

Minimum . | Average . | Maximum . | Minimum . | Average . | Maximum . | Minimum . | Average . | Maximum . | |

222 | 2 | 27 | 233 | 0.02 | 21.7 | 262 | 0.01 | 0.79 | 3.41 |

Number of droughts . | Drought duration (month) . | Drought severity . | Drought peak . | ||||||
---|---|---|---|---|---|---|---|---|---|

Minimum . | Average . | Maximum . | Minimum . | Average . | Maximum . | Minimum . | Average . | Maximum . | |

222 | 2 | 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 5–7, 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.

Drought characteristic . | Suitable marginal distribution . | Parameters 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 characteristic . | Suitable marginal distribution . | Parameters 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 |

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

. | Duration-Severity . | Severity-Peak . | Duration-peak . |
---|---|---|---|

Pearson coefficient | 0.96 | 0.87 | 0.90 |

Tau Kendall | 0.86 | 0.86 | 0.75 |

. | Duration-Severity . | Severity-Peak . | Duration-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 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.

Copula . | Duration-Severity . | Duration-Peak . | Severity-Peak . | Duration- Severity-Peak . | ||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

parameter . | AIC . | BIC . | Parameter . | AIC . | BIC . | parameter . | AIC . | BIC . | parameter . | AIC . | BIC . | |

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 |

Copula . | Duration-Severity . | Duration-Peak . | Severity-Peak . | Duration- Severity-Peak . | ||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

parameter . | AIC . | BIC . | Parameter . | AIC . | BIC . | parameter . | AIC . | BIC . | parameter . | AIC . | BIC . | |

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 |

Copula . | Duration-Severity . | Duration-Peak . | Severity-Peak . | Severity-Duration-Peak . | ||||
---|---|---|---|---|---|---|---|---|

Pvalu_SnB . | Pvalu_SnC . | Pvalu_SnB . | Pvalu_SnC . | Pvalu_SnB . | Pvalu_SnC . | Pvalu_SnB . | Pvalu_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 |

Copula . | Duration-Severity . | Duration-Peak . | Severity-Peak . | Severity-Duration-Peak . | ||||
---|---|---|---|---|---|---|---|---|

Pvalu_SnB . | Pvalu_SnC . | Pvalu_SnB . | Pvalu_SnC . | Pvalu_SnB . | Pvalu_SnC . | Pvalu_SnB . | Pvalu_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(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*(*D* ≥ *d* 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, R^{2}, 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.

Criteria . | R^{2}
. | Root mean squared error . | Nash–Sutcliffe efficiency . | Mean absolute error . |
---|---|---|---|---|

Value | 0.96 | 0.15 | 0.72 | 0.13 |

Criteria . | R^{2}
. | Root mean squared error . | Nash–Sutcliffe efficiency . | Mean 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.

Return period (year) . | Drought duration (month) . | Drought severity . | Drought peak . | Multivariate return period (T_{or})
. | Multivariate return period (T_{and})
. |
---|---|---|---|---|---|

4 | 1.27 | 0.015 | 0.007 | 4.4 | 4.7 |

5 | 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 severity . | Drought peak . | Multivariate return period (T_{or})
. | Multivariate return period (T_{and})
. |
---|---|---|---|---|---|

4 | 1.27 | 0.015 | 0.007 | 4.4 | 4.7 |

5 | 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 (*T*_{or} and *T*_{and}). 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, *T*_{and} (Equation (15)) is greater than *T*_{or} (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.

Also, in Figures 11–13, 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.

According to Figures 10–13, 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 11–13 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 (*T*_{or} and *T*_{and}) were determined via Equations (14) and (15). According to Table 13, the return period of historic GW droughts (*T*_{and}) 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.

Drought spell . | Duration (month) . | Severity . | Peak . | T_{or}
. | T_{and}
. |
---|---|---|---|---|---|

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 spell . | Duration (month) . | Severity . | Peak . | T_{or}
. | T_{and}
. |
---|---|---|---|---|---|

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.