Bivariate Drought Frequency Analysis in The Medjerda River Basin, Tunisia

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. *Corresponding author: Fateh chebana, INRS-ETE, 490 rue de la Couronne, Quebec (QC), G1K 9A9 Canada, Tel: 1-877-326-5762; E-mail: fateh.chebana@ete.inrs.ca Received April 12, 2016; Accepted April 26, 2016; Published 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/2165-784X.1000227 Copyright: © 2016 Hamdi Y, et al. This is an open-access 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.


Introduction
Drought is a recurrent feature of the North African climate. Due to its arid and semi-arid climate and to the uneven temporal and spatial distribution of rainfall in space and time, Tunisia is frequently affected by extreme hydro-climatic 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 socio-economic 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 1987-1989, 1993-1995 and 2000-2002 were very severe and resulted in significant negative economic and social impacts. In general, governmental support to drought-related 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 [3][4][5][6][7][8][9][10][11][12][13][14]. 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 semi-arid 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 non-parametric 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,[18][19][20][21][22][23][24][25][26]. 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 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 [11,22,[28][29][30][31]. 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 V d 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.

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 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. We extract all the events during which the streamflow is continuously below a certain threshold level Tr. Each drought event is characterized by its 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.  [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-α)/(n-2α+1) for the i th 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 i th drought event, V d is defined as: Where t bi and t ei are respectively the starting and ending times of the i th 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 semi-arid 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 so-called 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 drought-related 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 1 , u 2 , v 1 , v 2 ∈ I such that u 1 ≤ u 2 and v 1 ≤ v 2 : The link between copulas and bivariate distribution functions is provided by the Sklar's theorem [41] which states that the most general marginal-free 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 : 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 real-valued convex function A on the interval [0,1] such that: Where 0<u, v<1 and max {t, 1-t} ≤ A(t) ≤ 1. The case A ≡1 corresponds to independence, and A(t)=max{t,1-t} 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: 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: The Frank copula function is defined as: 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 drought-related 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 D corresponding to S and D are respectively expressed as: Where F S and F D are the univariate distributions of S and D.
The concept of return period of drought-related 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): The definition for the copulas-based drought events is as follows: 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 pseudo-likelihood (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 pseudo-likelihood and it can provide greater flexibility than the likelihood approach in the representation of real data. The Cramér-von Mises goodness-of-fit test and the associated p-value 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: Where ML denotes the maximum likelihood for the model. The model that possesses the minimum value of AIC should be selected.

Case Study: Basin and Data
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 north-eastern 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 semi-arid 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 large-scale 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 15th 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 sub-basin 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 .

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 semi-arid 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 m3/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 [51][52][53][54]. 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 T r = 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 1966-2008. The total deficit for the same threshold level is 1.66510 6 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 drought-related variables for the 15 largest droughts using three threshold levels (the Medjerda River at Jendouba 1966-       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, 1994-1996 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,[56][57][58].

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

Case Study: Results
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 drought-related frequency analysis: Exponential distribution [61], Generalized Pareto [62], Gamma distribution [63], Lognormal II and Lognormal III distribution [64] and Log-Pearson (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 goodness-of-fit tests and statistical criteria. More precisely, the Khi-squared χ 2 goodness-of-fit 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 p-values show that the Exponential (MLM), The Generalized Pareto (MOM and PWM), the Gamma (MLM) and the LP3 (MOM) can be selected as potential candidates (p-value<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 (p-value = 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: The LP3 distribution parameters, for S and D, are given in Table 5.

Copulas for drought severity and duration
In this study, a copula-based distribution representing the bivariate behaviour of droughts (severity-duration) 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, Husler-Reiss 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 drought-related variables should be modelled jointly.  Events number  The bold character indicates the accepted distribution. Copulas which consider only situations with relatively weak correlations (between -0.3 and 0.3), such as the Farlie-Gumbel-Morgensterm 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 log-likelihood 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 drought-related exercise, one can conclude, for instance, the following:
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 100-year 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    The bold character indicates the selected copula corresponding to the smallest AIC value.   -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.

Concluding Remarks
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 drought-related variables has therefore become necessary. The threeparameter Log-Pearson (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 moment-type and the maximum pseudo-likelihood approach are employed to estimate the copula parameters. The Cramér-von Mises goodness-of-fit test, the K-function 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 drought-related variables. The advantages of using the multivariate approach to model drought characteristics are also discussed.