^{1}Institute for Radiological Protection and Nuclear Safety, BP17, 92 262 Fontenay aux Roses Cedex, France
^{2}National Engineering School of Gabes, University of Gabes, Rue Omar IbnElkhattab 6029, Gabès, Tunisia
^{3}INRSETE, 490 rue de la Couronne, Quebec (QC), G1K 9A9 Canada
^{4}Institute Center for Water and Environment (iWATER), Masdar Institute of Science and Technology, Abu Dhabi, UAE
Received date: April 12, 2016; Accepted date: April 26, 2016; Published date: April 28, 2016
Citation: Hamdi Y, Chebana F, Ouarda TBMJ (2016) Bivariate Drought Frequency Analysis in The Medjerda River Basin, Tunisia. J Civil Environ Eng 6:227. doi:10.4172/2165784X.1000227
Copyright: © 2016 Hamdi Y, et al. This is an openaccess article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Visit for more related articles at Journal of Civil & Environmental Engineering
The climatology provides, for a given location or region, the time series of drought strength, the number, the mean duration, and the maximum duration of droughts of a given intensity. Similarly to most hydrological phenomena, droughts are characterized by a number of features such as their severity, duration and magnitude. Multivariate drought characterization has not been carried out in the various regions of the African continent despite the disastrous environmental, economic and social impacts of droughts. In the present paper, drought characteristics are modeled jointly in a multivariate frequency analysis (FA) framework for a data set from the Medjerda River, the principal watercourse in Tunisia. To identify drought events, the adopted threshold levels are estimated using the Flow Duration Curve (FDC) method. A sensitivity analysis to the threshold level is conducted. Results indicate that the drought features are significantly dependent and should be considered simultaneously for effective and rational modeling. Frank copula is shown to be the most appropriate copula model to represent drought features for the considered data set. The joint probabilities and bivariate return periods, based on the developed two dimensional copula models, are estimated in order to evaluate the contribution and advantages of bivariate modeling of droughts. These results are of practical relevance to hydrologists and water resources managers in Tunisia for applications in drought risk analysis and drought management, and in general for the optimal planning and management of water resources systems.
Copula; Duration; Drought characteristics; Multivariate frequency analysis; Medjerda river basin; Severity; Tunisia
Drought is a recurrent feature of the North African climate. Due to its arid and semiarid climate and to the uneven temporal and spatial distribution of rainfall in space and time, Tunisia is frequently affected by extreme hydroclimatic events such as floods and droughts. Generally, a drought event, followed by relatively severe floods, occurs in Tunisia once every 10 years [1]. These natural disasters have caused major drawbacks to the economic, social, environmental and ecological development in Tunisia throughout history. Recent severe and prolonged droughts have highlighted Tunisian’s vulnerability to this natural hazard and alerted the public, the government, and operational agencies to the many socioeconomic problems accompanying water shortage and to the need for drought mitigation measures. The impacts of these events are being aggravated by the ever increasing water demand. A number of strategies relative to water resources planning and management were developed by the Tunisian government over the last three decades. Despite these efforts, these extreme climate events continue to cause serious damages to the Tunisian economy. These challenges warn of the critical need for improved scientific investigations to meet current and future water requirements in Tunisia. Drought studies have been suffering from the lack of consistent methods for drought analysis. The first step in a drought analysis would be to define the drought event. The scientific community has only agreed on general definitions of a drought event. A drought is generally defined as a sustained period of significantly lower rainfall events, low flows and soil moisture levels relative to usual levels. Tallaksen and Van Lanen [2] defined droughts as “a sustained and regionally extensive occurrence of below average natural water availability”. During drought periods, water supply can be inadequate to meet water demand. A large number of economic and environmental activities can also be seriously affected by this lack of water availability. A multivariate approach for the analysis of the frequency of such events is presented in the present paper. In this work, droughts are related to streamflow deficits. The volume, duration and magnitude of low flows are treated as random variables. Expanding populations, increasing land and water use, lack of rainfall, high temperatures and high evaporation and evapotranspiration rates are the major factors and causes of droughts in Tunisia [3]. Severe drought events have been occurring frequently in Tunisia with considerable impacts. The last three drought events that occurred during 19871989, 19931995 and 20002002 were very severe and resulted in significant negative economic and social impacts. In general, governmental support to droughtrelated studies, mitigation and adaptation in Tunisia remains insufficient and needs to be increased. More efforts need to be devoted to statistical modeling of historical drought data in Tunisia in order to acquire a better understanding of drought characteristics, dynamics and evolution in the region. Since droughts are stochastic in nature, the proper way to explore droughts is by using probabilisticbased approaches. Following Yevjevich, who introduced the theory of runs to analyze droughts, a large number of researchers investigated the probabilistic properties of droughts [314]. Frequency estimation in hydrology is commonly used for the analysis of extreme hydrological events. The Principal aim of frequency analysis (FA) is identifying a relationship between the magnitude or the severity of extreme events and their probability of occurrence. This relationship can be obtained using the experimental probabilities and distribution functions [15]. Ouarda et al. [13] provided a review of statistical models commonly used in the estimation of low flows. A number of research efforts on drought FA in arid and semiarid regions were conducted during the last two decades. For instance, Eltahir [16] used available rainfall series in central and western Sudan to analyse drought frequency in the region. Hallack and Watkins [17] investigated how drought intensityduration frequency estimates and seasonal drought forecasts can be useful for farmers and other water managers in northern Mexico. Kim et al. [10] proposed a nonparametric approach for estimating return periods of droughts in arid regions. Note that the effectiveness of FA techniques was limited by the quality and length of available data series. Most of the previous studies focused separately on one or more of the drought characteristics in a univarite analysis framework. It has been extensively shown in the literature that extreme hydrologic events, such as droughts, are multivariate events [4,14,1826]. Yoo et al. [27] developed a practical drought frequency analysis method based on a bivariate distribution by incorporating regional drought attributes that are associated with drought frequency (e.g., duration and severity). These events are characterized by a number of random variables which are mutually correlated. These variables are often selected to be the volume Vd or severity S representing the water deficit below a selected threshold Tr, the duration D of the event and the magnitude M defined as the ratio of S over D [11,22,2831]. Therefore a better approach for describing drought events is to derive the joint distribution of these variables. The remainder of the paper is organized as follows. The multivariate drought characterization approach is presented in section 2. In this section, the definitions of a drought threshold level, droughtrelated variables and return period in drought analysis are presented. In section 3, the case study and the data base are presented. The univariate and multivariate FA procedures are presented in section 4. Multivariate drought characterisation using copulas and results are provides and discussed in section 5. Finally, a number of conclusions are presented and promising other research directions in drought modelling is suggested.
Multivariate drought: characterisation and modelling
Using the multivariate approach to deal with extreme hydrological events, such as floods and droughts, instead of the univariate analysis was justified by several investigators. In recent years, a number of research and review papers have addressed multivariate modeling of extreme events in hydrology. It was also shown in the literature that univariate hydrological FA provides a limited assessment of an extreme event probability of occurrence. In addition, the univariate analysis of event characteristics does not take into account their dependence structure and this reduces the risk estimation accuracy. The univariate approach could be justified when only one variable is significant or when variables are independent. A drought is, however, a multivariate event characterized by the variables mentioned earlier in the introduction and presented in Figure 1 (volume Vd or severity S, duration D and magnitude M). Shiau and Shen [32] suggested that a better approach for describing drought characteristics is to derive the joint distribution of D and S. A bivariate distribution is thus appropriate to describe these correlated hydrologic variables. De Michele et al. [33] defined a drought event as extreme expressions of the river flow dynamics and episodes during which the streamflow is below a given threshold. They described droughts as multivariate events characterized by two variables: D and S. In their work, a new concept of dynamic return period was introduced and formulated using the theory of Copulas and calculated via a Survival Kendall’s approach.
Figure 1: Definition sketch of drought events: This figure illustrates how drought events are defined. We extract all the events during which the streamflow is continuously below a certain threshold level T_{r}. Each drought event is characterized by its volume V_{d} (or severity S) representing the water deficit below a selected threshold T_{r}, the duration D of the event and the magnitude M defined as the ratio of S over D.
Definition of a drought threshold
In this study, a drought is defined, using the threshold level method [33], as an event during which the streamflow is continuously below a certain threshold level T_{r} as shown in Figure 1. The threshold level method has already been evaluated for its applicability to daily discharge series for streams in different climate zones and with different hydrological regimes. Shiau [12] used the Standardized Precipitation Index (SPI) method to define a drought event where the threshold corresponds to the median. They used the SPI to quantify the precipitation deficit in terms of the probability for multiple time scales. Positive or negative values of SPI indicate respectively that the observed precipitation is higher or lower than the associated median. The wet and dry conditions are, however, classified according to this median which represents a threshold level. Chung and Salas [7] converted the original flows into a sequence of zeros and ones based on a threshold flow equal to the sample mean. The value 0 denotes the dry state and 1 the wet state, depending on whether flows are smaller or greater than the mean. One of the most fundamental hydrological characteristics is the Mean Annual Runoff (MAR). It is the mean value of available annual flow totals time series. When dividing the MAR by time we obtain the longterm mean daily discharge called the Mean Daily Flow (MDF). Various streamflow drought indices may be expressed as a percentage of either MAR or MDF. The Median Flow (MF), which is the middle value in a ranked flow time series, can also represent a conservative upper bound for low flows that can be considered as a threshold level. For instance, Smakhtin [8] argued about the suitability of the MAR and the MF to distinguish between droughts and low flows. Various other streamflow drought indices may be estimated from the Flow Duration Curve (FDC) which is one of the most informative methods of displaying all discharge values including droughts and floods. FDC can be constructed using different time scales (e.g. daily, monthly and annual). It was shown in the literature that FDC based on daily flow time series provides the most reliable way of analysing event characteristics (e.g. Searcy [34]; Tallaksen et al., [2]; Smakhtin, [8]). A brief presentation of FDC is given below and a detailed description of FDC construction and interpretation can be available in a number of references (e.g. Searcy [34]; Vogel and Fennessey [35]; Smakhtin et al., [8]). FDC plots the empirical frequency of streamflow as a function of the percentage of time that the streamflow is equalled or exceeded. To construct the curve, data are ranked in a decreasing order, and for each value the probability of exceedance is computed using a probability plotting position formula. Many of the positions commonly recommended are of the form (1α)/(n2α+1) for the ith observation within a sample of size n and the value of the constant α refers to specific formulas such as α= 0.5 (Hazen), α = 0.375 (Blom), and for normal probability plots; α = 0 (Weibull and Gumbel) (see e.g. Rao and Hamed, [36]). Other rules (Gringorten, Cunnane) are tested in the present study. The corresponding results (not presented here due to space limitations) show clearly that the value of the constant α does not have a notable effect on flow duration curves. The threshold level is considered as the boundary between usually and unusually low streamflows and it is chosen based on the characteristics of the streamflow regime. In the case of perennial and intermittent streams, indices such as percentiles from the FDC, are frequently applied. For perennial streams, the more commonly used low flow indices that can be easily obtained from a FDC are between the 90 and 70 percentiles respectively denoted Q90 and Q70 [2,8,10,37,38].
Drought severity, duration and magnitude
For a given threshold level, drought events are illustrated by the shaded sets in Figure 1. As mentioned in section 2 of the present paper, each drought event is composed of a drought duration D, severity S (or drought volume V_{d} ) and magnitude M. In the literature V_{d} is also called deficit on the negative run sum. For the ith drought event, V_{d} is defined as:
(1)
Where t_{bi} and t_{ei} are respectively the starting and ending times of the ith drought event and Q_{t} is the discharge at the time t. From the selected threshold level and the available record of daily discharges and after computing all deficits and drought event durations, some of the identified droughts could be very close to each other. These represent minor drought events (events 1 and 2 in Figure 1 for instance), and therefore could be mutually dependent. This situation could occur for instance for short time scales and this is often the case for arid and semiarid regions. Note that the presence of a considerable amount of minor events could have a negative impact on the analysis. To overcome this difficulty, two assumptions are made [39]. First, we consider only events where S is larger than a given minimum value, V_{min} which can be fixed according to the case study. The second assumption is related to deficits larger than the socalled very minor ones. It can happen that during a drought event the discharge exceeds the threshold level for a very short period (Figure 1) and thus splits a practically single streamflow drought into two events (events 1 and 2 in Figure 1). In this case, it can be considered that these events represent only one deficit.
Multivariate drought characterisation using copulas
When the observed S and D are highly correlated (which is the case in the present study), the construction of the joint distribution of droughtrelated variables becomes necessary. Multivariate drought characterisation and copulas used in constructing the joint distribution are discussed in the present section. To overcome the limitations of classical dependence measures, copulas have recently received increasing attention in a number of fields (see for instance Nelsen [40]). A copula is a description and a model of the dependence structure between random variables, independently of the marginal laws. A general description of copulas theory can be found, for instance, in Nelsen (2006).
A copula is a function C: I X I→I(I=[0,1]) such that:
For all u, v ∈ I: C(u,0)=0, C(u,1)=u, C(0,v)=0, C(1,v)=v
For all u_{1}, u_{2}, v_{1}, v_{2} ∈ I such that u_{1} ≤ u_{2} and v_{1} ≤ v_{2}:
C (u_{2},v_{2})C(u_{2},v_{1})C(u_{1},v_{2})+C(u_{1},v_{1}) ≥ 0 (2)
The link between copulas and bivariate distribution functions is provided by the Sklar’s theorem [41] which states that the most general marginalfree description of the dependence structure of multivariate distributions is through its copula. We consider herein the situation with two marginal distribution functions F_{1} and F_{2} of the random variables X_{(1)} and X_{(2)} respectively. Let F_{1,2} denote the joint distribution function with marginal distributions F_{1} and F_{2}. Then, there exists a copula C such that, for all real x_{1} and x_{2}:
F_{1,2}(x_{1},x_{2})=C(F_{1}(x_{1}),F_{2}(x_{2})) (3)
Furthermore, if F_{1} and F_{2} are continuous, then C is unique. Archimedean and Extreme Value (EV) copulas represent classes of particular interest. A copula C is an EV copula if and only if there exists a realvalued convex function A on the interval [0,1] such that:
(4)
Where 0<u, v<1 and max {t, 1t} ≤ A(t) ≤ 1. The case A ≡1 corresponds to independence, and A(t)=max{t,1t} Corresponds to the copula C(u,v) = min(u,v). Statistical inference on Pickand’s function A can summarize the inference on its bivariate EV copula C. A bivariate Archimedean copula is characterized by a generator ѱ(.) that is a convex decreasing function satisfying Ѱ(1)=0 where:
C(u,v) = ψ (ψ(u) +ψ(v)), 0 < u, v <1 (5)
Copulas that belong to this class are symmetric and associative. The Frank, Clayton and Gumbel logistic models are the simplest and most popular copulas. For instance, Frank copula corresponds to:
(6)
The Frank copula function is defined as:
(7)
Note that the Gumbel copula is the only one which is simultaneously Archimedean and EV copula.
Multivariate return periods for droughts
The notion of return period is commonly used by hydrologists and civil engineers and represents an important variable for the design of water resource systems and infrastructures. The risk related to droughts is the probability that one or more droughts with a certain D, S or M_{d} occur during the life time of the project [42]. One can define the return period in drought analysis in different ways for different applications. When applied to droughtrelated variables, the concept of return period indicates the average time between the occurrences of critical events [9,43]. Douglas defined the return period in drought estimation and analysis as the average number of trials required to the first occurrence of a critical event. An alternative definition of return period is the average time between the occurrence of events with a certain S or less [44]. The return periods T_{S} and T_{DM} corresponding to S and D are respectively expressed as:
(8)
Where F_{S} and F_{D} are the univariate distributions of S and D.
The concept of return period of droughtrelated variables was introduced recently for the multivariate context [32,45]. For instance, Shiau [32] proposed a methodology that categorizes the return periods of multivariate hydrologic events as joint and conditional return periods. We can define the joint S and D return periods in two ways: (i):
(9)
The definition for the copulasbased drought events is as follows:
(10)
where SD F denotes the two dimensional distribution function with univariate distributions FS for S and FD for D. It is important to notice that various S and D combinations can result in the same value of the return period. Copulas dependence parameters can be obtained employing the method of moments based on the inversion of Spearman’s ρ and Kendall’s τ. It was shown in the literature that such approach may lead to inconsistencies [46]. More reliable results can be obtained using the maximum pseudolikelihood (MPL) approach [47]. It is based on the information contained in the moments of the first and second order of the endogenous variables. This approach consists in maximizing the log pseudolikelihood and it can provide greater flexibility than the likelihood approach in the representation of real data. The Cramérvon Mises goodnessoffit test and the associated pvalue of the estimate, for the same threshold level, are also used in the present study. In addition, the Akaike criterion (AIC) is used to discriminate between copulas. This last criterion maximises the loglikelihood over the obtained parameter and is given by:
AIC = −2ln(ML) + 2 (11)
Where ML denotes the maximum likelihood for the model. The model that possesses the minimum value of AIC should be selected.
Tunisia is the country with the smallest area in North Africa. Tunisia occupies a privileged geographic position at the crossroads of the Eastern and Western basins of the Mediterranean, between Europe and Africa (Figure 2). It is located at the northeastern tip of Africa and is bordered by the Mediterranean to the North and the East, by Libya to the South, and by Algeria to the West. Tunisia is characterised by a Mediterranean climate with hot dry summers and cool moist winters. Precipitation is very irregular and the rainfall varies considerably from the North to South. The rainy season extends from September to May, with relatively intense precipitations in the autumn. Tunisia is divided into four large geographical units, namely the Northern, Eastern, Central and Southern regions. Most of the land is either in a semiarid or arid zone. There are five bioclimatic zones in going from the most arid to the most humid based on rainfall [48]. Rainfall and temperature (especially the winter temperature) are the most important bioclimatic determinants. These are not only governed by altitude but by the degree of continentally. Inland regions have hotter summers and colder winters than coastal areas which benefit from the buffering effects of the sea. Bio climatically, therefore, the country is also divided into areas of warm, cool and cold winters. Every few years, the country experiences a largescale drought. During the last few decades, Tunisia has experienced more than three severe drought events, which may be considered as extreme. Rainfall deficits were observed in Tunisia during 1987–1989, whereas more severe droughts were recorded during the 1981–1983, 1991–1995 and 1999–2002 periods. According to Touchan et al. [49,50], the most recent drought of 1999–2002 appears to be the worst since at least the middle of the 15^{th} century. More than 20 years presented water deficit on a total of 32 years of survey (between 1972 and 2003). The area under study is a subbasin of the Medjerda river basin and is defined as the area draining into the gauging station of Jendouba (Figure 2). The drainage area of the Medjerda River basin at this location is 2420 km^{2}.
Figure 2: Location of the study region (Medjerda River basin in Tunisia) and station (Jendouba). It can be seen from this figure that Tunisia is located at the northeastern tip of Africa and is bordered by the Mediterranean to the North and the East and by Algeria to the West. The Medjerda River originates in the Atlas Mountains of eastern Algeria and then flows eastwards to Tunisia to finally enter the reach the Mediterranean Sea. The entire Medjerda catchment covers approximately 24,000 km^{2}, of which 16,300 km^{2} (32%) are located in Tunisia.
The Medjerda River basin
The Medjerda River (also known as Oued Medjerda) is the longest (with a length of 450 km) and most important river in Tunisia. It originates in the semiarid Atlas Mountains of eastern Algeria, and then flows eastwards to Tunisia (Figure 2), to finally enter the Gulf of Utica in the Mediterranean Sea. The entire Medjerda catchment covers approximately 24,000 km^{2}, of which 16,300 km^{2} (32%) are located in Tunisia. The Medjerda River is Tunisia’s principal watercourse and constitutes the water supply source for more than half of the Tunisian population (more than 5 million inhabitants out of 10 millions). It is also the major supplier of water to the country’s wheat crops [51]. The food security of the entire country relies on the Medjerda River inflows. On average, the water resources of the Medjerda River represent 1 billion m^{3}/year or 37% of the surface water in Tunisia and 22% of the country’s renewable water resources [51]. Agriculture is the primary economic activity in the interior of the basin with more than 100,000 ha irrigated with water from the Medjerda River. The River is also used for domestic water supply to several regions in the country.
Data
Daily streamflow data for the Jendouba station for the period 1966–2008 are utilized in the present study. The discharge data were provided by the General Direction of Water Resources (DGRE, Direction Générale des Ressources en Eau). Streamflow data from the Jendouba station was used by several investigators [5154]. On the basis of these daily streamflow, and using a selected threshold, drought events are identified and the corresponding features are derived. More details are given in the following section. The empirical FDC for the Medjerda River with the Hazen rule is shown in Figure 3 and corresponding low flow indices are presented in Table 1. When a threshold level of Tr = 0.387 m^{3}/s is used (which represents the 90% low flow from the daily FDC) a total of 32 drought events are detected during the period of 19662008. The total deficit for the same threshold level is 1.665106 m^{3}/s. The number of events and the total deficit corresponding to a number of threshold levels are also shown in Table 1. D, S and M are extracted from the observed drought events. M was ranked in a descending order to extract the largest observed droughts for each threshold level. Table 2 summarizes these droughtrelated variables for the 15 largest droughts using three threshold levels (the Medjerda River at Jendouba 19662008). It is shown from these Tables that since 1966 the three worst drought events, with regards to magnitude and for the Q90 threshold level, have occurred in the Medjerda River basin in 2003, 1996 and 1995 respectively. The results of Table 2 confirm the testimonies of a large number of farmers who experienced these drought events (Hamdi, 2009). It was also shown in other references [55,56] that events that took place in 2003, 19941996 and 1990 were the worst to occur in Tunisia over the last century. The drought events which occurred during the years 1987 and 2002 were pointed out by many investigators [1,5658]. The importance of the other events presented in Table 2 was noted by other investigators [49,56,59,60]. Figure 4 shows a histogram of drought characteristics (S and D) for the Medjerda River basin considering the Q90 threshold level only. Other statistical characteristics of S and D of the Medjerda River basin considering the threshold levels Q90, Q85 and Q80 are given in Table 3.
Threshold level (m^{3}/s)  Number of events  Total deficit (x 10^{7}) (m^{3})  Correlation coefficient R(D,S)  

Q90  0.387  32  1.665  0.88 
Q85  0.466  48  2.943  0.87 
Q80  0.542  57  4.653  0.9 
Table 1: Threshold levels and corresponding streamflow drought characteristics.
Date  M (m^{3}/s)  S(10^{6}) (m^{3})  D(days)  Date  M (m^{3}/s)  S(10^{6}) (m^{3})  D (days)  M (m^{3}/s)  S(10^{6})(m^{3})  D (days) 

2003  0.257  1.222  55  2003  0.328  1.618  57  2003  0.398  58 
1996  0.225  1.342  69  1996  0.296  1.819  71  1996  0.344  78 
1995  0.216  1.7  91  1990  0.27  2.166  93  1990  0.338  95 
1994  0.21  1.722  95  1995  0.259  2.375  106  1987  0.326  111 
1990  0.202  1.55  89  1994  0.255  2.425  110  1981  0.269  101 
1987  0.175  1.659  110  1987  0.251  2.408  111  1980  0.254  80 
1980  0.147  0.789  62  1980  0.199  1.253  73  1995  0.253  149 
1981  0.131  1.041  92  1981  0.198  1.697  99  1991  0.245  105 
2001  0.121  0.313  30  1991  0.177  1.544  101  2009  0.233  92 
1991  0.112  0.893  92  2009  0.169  1.271  87  1998  0.225  114 
2009  0.108  0.712  76  2002  0.155  0.886  66  2002  0.214  72 
1984  0.103  0.634  71  1998  0.153  1.472  111  1984  0.201  122 
2002  0.094  0.469  58  1993  0.143  0.839  68  1992  0.199  87 
1998  0.093  0.766  95  1984  0.136  1.335  114  1993  0.193  82 
1993  0.081  0.413  59  1992  0.131  0.937  83  1989  0.182  50 
Table 2: Streamflow M, S and D for the 15 magnitudebased largest droughts using three threshold levels (Medjerda River at Jendouba 19662008).
Q90  Q85  Q80  

S (m^{3})  D (days)  S (m^{3}) 
D (days)  S (m^{3})  D (days) 

Minimum  1.13E+04  5  1.52E+04  6  4.80E+03  6 
Maximum  1.72E+06  110  2.42E+06  114  3.86E+06  166 
Mean  5.22E+05  46.83  6.13E+05  47.5  8.16E+05  50.8 
Median  2.78E+05  45.5  1.87E+05  39.5  3.32E+05  47 
Standard Deviation  5.79E+05  63.64  7.53E+05  36.2  9.90E+05  40.3 
Kurtosis  2.47  1.63  3.06  1.75  3.74  2.87 
Skewness  0.92  0.29  1.17  0.47  1.31  0.8 
Table 3: Statistical characteristics of S and D for different threshold levels.
In this section, we present the results of the application of the methods developed in the methodology section to the case study dataset.
Fitting drought severity and duration margins
The main objective of a FA is to relate the magnitude of the events of interest to their frequencies of occurrence through the use of probability distributions. For each selected threshold level, the number of drought events is listed in Table 1. The values of D and S corresponding to three threshold levels are extracted from observed drought events and are presented in the previous section. Univariate FA explores separately each variable in a data set. The following distributions have been commonly used for droughtrelated frequency analysis: Exponential distribution [61], Generalized Pareto [62], Gamma distribution [63], Lognormal II and Lognormal III distribution [64] and LogPearson (LP3) distribution [65]. These distributions are therefore considered in the present study for fitting S and D, and the appropriate choice for each variable is determined based on goodnessoffit tests and statistical criteria. More precisely, the Khisquared χ^{2} goodnessoffit test and the Akaike Information Criterion (AIC) are performed. Parameters of the univariate distributions are estimated using the method of moments (MOM), maximum likelihood method (MLM) and Probability Weighted moments (PWM). Results of the test and criteria mentioned above are summarized in Table 4. The χ^{2} test pvalues show that the Exponential (MLM), The Generalized Pareto (MOM and PWM), the Gamma (MLM) and the LP3 (MOM) can be selected as potential candidates (pvalue<5%) to model D. Among these distributions, the LP3 (MOM) gave the smallest value of AIC (306.738). The LP3 distribution is therefore selected to model drought duration for the Medjerda River basin. Similarly, we can deduce that the LP3 (MOM) presented an appropriate fit for drought severity (pvalue = 9.16% and AIC = 903.763). Figure 5 shows the observed D and S and obtained distribution functions for the LP3 (MOM) distribution. The closeness between the observed data and fitted distribution for both figures (S in the top panel and D in the bottom panel) indicates that the three parameters LP3 distribution fitted with MOM is appropriate to model both S and D but with different parameter values. The probability density function of the LP3 distribution with its parameters α, λ and m is given by the following equation:
AIC  c2  p  value  

S  D  S  D  S  D  
Exp. (MLM)  911.275  309.004  22.5  9  0.0004  0.1091 
G. Pareto (MOM)  Nd  307.433  24  9.5  0.0002  0.0907 
G. Pareto (PWM)  Nd  308.048  17  10.5  0.0045  0.0622 
Gamma (MOM)  906.336  311.729  13.5  14  0.0191  0.0156 
Gamma (MLM)  904.297  310.787  9.5  7  0.3579  0.2206 
LN2 (MLM)  906.295  313.607  11  16  0.0514  0.0068 
The bold character indicates the accepted distribution.
Table 4: Univariate frequency analysis results for S and D for threshold level Q90.
The LP3 distribution parameters, for S and D, are given in Table 5.
Q90  QQ85  QQ80  

S  D  S  D  S  D  
a  0.688363  1.746922  0.635839  1.837829  0.944104  2.790086 
l  0.872408  1.037528  0.952772  1.217317  1.278521  2.102812 
m  6.274639  2.049483  6.420921  2.106088  6.597664  2.255147 
Table 5: LP3 distribution parameters.
Figure 5: Observed probabilities with fitted LP3 distribution for severity S and duration D. It shows the empirical probabilities for D and S and the fitting with the LP3 distribution functions (using the method of moments MOM for parameters estimation). The closeness between the observed data and fitted distribution for both figures indicates that the three parameters LP3 distribution fitted with MOM is appropriate to model both S and D but with different parameter values.
Copulas for drought severity and duration
In this study, a copulabased distribution representing the bivariate behaviour of droughts (severityduration) in the Medjerda River basin is derived. Several types of copulas are tested to identify the most appropriate one. The candidate copulas are: the Gumbel, Clayton, Frank, Galambos, Normal, HuslerReiss and Plackett. As indicated earlier in this paper, the significant correlations observed between S and D (correlation = 0.88, Kendall’s τ = 0.71 and Spearman’s ρ = 0.89) suggest that droughtrelated variables should be modelled jointly. Copulas which consider only situations with relatively weak correlations (between 0.3 and 0.3), such as the FarlieGumbelMorgensterm copula [66], are not considered in the present paper. Results of the MPL method (for Q90 threshold level) are reported in Table 6. They indicate that relatively minor differences exist in the joint probabilities of the 7 tested copulas (with Q90 threshold level). The value of the loglikelihood function of the Frank copula is equal to 23.31 and this value is larger than those of other copulas. The minimum value of the AIC (44.62) corresponds to the Frank copula as well. This copula is then selected as the potential candidate to model droughts in the Medjerda River basin. The comparison between the obtained copula and observed droughts is presented in Figure 6. Although comparison with observed droughts is difficult to illustrate, we can still see in Figure 6 that observed data are well represented by the Frank copula. A number of authors have underlined in the literature the adequacy of Frank copula to represent drought data [67]. Figure 7 illustrates the contours of the Frank copula with the Q90 threshold level. According to the LP3 distribution function (Eq. 12) and the corresponding estimated parameters given in Table 5, S and D quantiles corresponding to return periods of T = 2, 5, 10, 20, 50 and 100 years, and to threshold levels of Q90, Q85 and Q80, are summarized in Table 7. Figure 8 is an illustration of contours of joint return periods of S and D defined by Eq. 3. Figure 8 can be used by practitioners to obtain results such as: In the case of the Medjerda River basin, for a given threshold level (for example Q90) and depending on the droughtrelated exercise, one can conclude, for instance, the following:
• The 100year drought severity varies between V_{115,100} = 1,85.10^{6} m^{3} and V_{110,100} = 1,9.10^{6} m^{3}.
• The 50year drought severity varies between V_{115,50} = 1,8.10^{6} m^{3} and V_{108,50} =1,9.10^{6} m^{3}.
• The 10year drought severity varies between V_{115,10 }= 1,49.10^{6} m^{3} and V_{95,10} = 1,9.10^{6} m^{3}.
Gumbel  Clayton  Frank  Galambos  Normal  HuslerReiss  Plackett  

GOF test pvalue  0.0135  0.0015  0.1173  0.0095  0.0245  0.0135  0.0145 
Parameter (MPL)  2.89  2.73  11.88  2.19  0.89  2.88  28.58 
Loglikelihood  19.9  17.5  23.31  20.08  22.92  20.53  21.18 
AIC  37.79  33.07  44.62  38.16  43.84  39.06  40.36 
The bold character indicates the selected copula corresponding to the smallest AIC value.
Table 6: Copulas parameter estimation with the MPL method (threshold Q90).
Return period  Q90  QQ85  Q80  

(years)  
S (m^{3})  D (days)  S (m^{3})  D (days)  S (m^{3})  D (days)  
2  2.79E+05  43  2.53E+05  41  3.77E+05  41 
5  1.09E+06  81  1.28E+06  83  1.60E+06  86 
10  1.49E+06  96  1.89E+06  102  2.41E+06  111 
20  1.71E+06  104  2.26E+06  113  3.00E+06  129 
50  1.83E+06  109  2.50E+06  121  3.49E+06  147 
100  1.86E+06  111  2.58E+06  124  3.71E+06  156 
Table 7: Return periods of S and D separately.
Similarly, for future severe streamflow droughts in the Medjerda River, one can use Figure 8 to determine the return periods corresponding to the values of the severity and duration. Finally, it can be concluded that in terms of return values estimates the bivariate estimation procedure performs comparably to the univariate procedure. For instance, the 10, 50 and 100year severity estimates (1, 49.10^{6} m^{3}; 1, 83.10^{6} m^{3} and 1, 86.10^{6} m^{3}) obtained when using the univariate FA (Table 7), belong to the ranges estimated by the Multivariate analysis [1,49  1,9].10^{6}; [1,25  1,9].10^{6} and [1,8  1,9].10^{6} m^{3} respectively. However, in addition to the fact that the univariate estimation does not consider the variation in the dependence structure of the phenomenon, it gives us single severity values regardless of events durations. It means that two drought events with the same severity and two different durations (20 days and 100 days for instance) will have the same probability of occurrence. This illustrates the major limitation of the univariate framework and how the information it provides is practically erroneous.
A bivariate drought frequency analysis study in the Medjerda River basin, Tunisia, is presented. A statistical model based on copulas is developed in this study to compute the joint drought severity S and drought duration D probability. Drought severity and duration are important characteristics and are often used as design variables for water management infrastructures. However, S and D are usually analyzed separately. Since the literature has shown that droughts cannot be modeled with separate univariate characteristics, bivariate probabilistic proprieties are derived in this study. Droughts in the Medjerda River basin in Tunisia are used as case study. The basic data is obtained from specific threshold levels which are considered as boundaries between usually and unusually low streamflow values. Streamflow drought indices and threshold levels are estimated employing the Flow Duration Curve. A very significant relationship exists between the observed D and S with a correlation coefficient greater than 85%. The construction of a joint distribution for these droughtrelated variables has therefore become necessary. The threeparameter LogPearson (LP3) distribution function using the method of moments (MOM) gave the best fit to both D and S. Copulas are employed to construct the dependence structure for D and S. Several types of copulas are tested to select the most appropriate one. The method of momenttype and the maximum pseudolikelihood approach are employed to estimate the copula parameters. The Cramérvon Mises goodnessoffit test, the Kfunction graphical test and the Akaike criterion (AIC) are used to discriminate between copulas. It was found that the Frank copula provides the best fit in comparison to other copulas, for the observed drought data of the Medjerda River basin. The contours of the Frank copula and the joint probabilities are also investigated and presented. It was demonstrated in this paper that copulas can be easily applied to construct the dependence structure of multivariate high correlated droughtrelated variables. The advantages of using the multivariate approach to model drought characteristics are also discussed.
Financial support for this study was graciously provided by the Natural Sciences and Engineering Research Council of Canada (NSERC). The authors would like to acknowledge the contribution of the DGRE (Direction Générale des Ressources en Eau) in Tunisia for providing the streamflow level series.3
All the authors declare they have no conflict of interest.
Make the best use of Scientific Research and information from our 700 + peer reviewed, Open Access Journals