Drought assessment in a south Mediterranean transboundary catchment

ABSTRACT

Predicting the impacts of climate change on water resources remains a challenging task and requires a good understanding of the dynamics of the forcing terms in the past. In this study, the variability of precipitation and drought patterns is studied over the Mediterranean catchment of the Medjerda in Tunisia based on an observed rainfall dataset collected at 41 raingauges during the period 1973–2012. The standardized precipitation index and the aridity index were used to characterize drought variability. Multivariate and geostatistical techniques were further employed to identify the spatial variability of annual rainfall. The results show that the Medjerda is marked by a significant spatio-temporal variability of drought, with varying extreme wet and dry events. Four regions with distinct rainfall regimes are identified by utilizing the K-means cluster analysis. A principal component analysis identifies the variables that are responsible for the relationships between precipitation and drought variability.

KEYWORDS:

• standardized precipitation index (SPI)
• aridity index (AI)
• multivariate statistical methods
• Mediterranean region
• drought

Editor S. Archfield Associate editor F.-J. Chang

1 Introduction

There is ever-increasing evidence that global warming is becoming more and more important as the global concentration of carbon dioxide in the atmosphere continues to rise due to anthropogenic emissions (Hsiang and Kopp 2018). Global warming causes major changes in the water cycle. It affects the availability of water resources, energy production, agriculture, hydrological regimes of rivers, and also increases the risk of extreme events (Donnelly et al. 2017, Sunde et al. 2017). Simultaneously, global water demand is also increasing, which exacerbates substantially the drought risk (IPCC (Intergovernmental Panel on Climate Change) 2014, Bhuyan et al. 2017). In this context, droughts have been portrayed as a prolonged dry weather period caused by a lack of precipitation, which leads to severe hydrological imbalances at different scales (WMO 2006, Kreibich et al. 2019).

At the regional scale, the Middle East and North Africa region have been identified as one of the hot spots that will be affected by extreme heat and increasing water demand as well as drought (Cook et al. 2016, Waha et al. 2017). According to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC), substantial changes in average rainfall are expected for this region (Bento et al. 2018). Total precipitation for this region is expected to decrease while the frequency and intensity of extreme events are expected to increase (Gara et al. 2019). Hence, based on these predictions, the succession of dry years with the concentration of prolonged drought and exceptional floods become more likely to occur in this region (Tramblay and Hertig 2018).

The above mentioned hydrometeorological alteration is a specific point of concern for Tunisia. This country is situated in the semi-arid Mediterranean region which is highly affected by climate change (Chargui et al. 2018, Moussa et al. 2019). It is characterized by large spatio-temporal variable rainfall patterns and overall scarce freshwater resources (Gader et al. 2015). Water resources are limited and unevenly distributed due to the Mediterranean and arid climate of the country. In this context, it is crucial to study the spatio-temporal dynamics of precipitation and drought.

Studies on climate variability in the Mediterranean region have focused on precipitation trends and distribution (Liuzzo et al. 2016, Blanchet et al. 2018). Statistical distributions such as the generalized extreme value distribution are used to characterize the variability (Nasri et al. 2016, Ben Khalfallah and Saidi 2018). Parametric and non-parametric statistical tests like the Mann-Kendall test are used to identify trends (Liuzzo et al. 2016, Sahany et al. 2018). Other studies have used deterministic and (geo)-statistical techniques to interpolate the space-time distribution of precipitation (Feki et al. 2012, Gupta et al. 2017). These studies allowed to assess the spatial patterns of precipitation but did not allow having a clear process understanding of the complexity of the rainfall regime and the factors that control its heterogeneity.

Other studies have focused on drought using different drought indices and assessment methods (Ionita et al. 2016, Frank et al. 2017, Hui-Mean et al. 2018). The standardized precipitation index (SPI), for instance, is used by many national meteorological and hydrological services (Svoboda et al. 2012, Ionita et al. 2016). This index is based on historical precipitation records (Eyring et al. 2016, Mathbout et al. 2018). The aridity index (AI) is also used to assess drought. This index takes into account the impacts of temperature and precipitation on drought (Bharadiya et al. 2015). Yet, the AI robustness is still poor to study drought characteristics in a warming climate (Mukherjee et al. 2018). Other studies have tried to combine different drought indices (Mukherjee et al. 2018, Yihdego et al. 2019). Finally, as drought is strongly affected by the spatial and temporal variability of the drought forcing terms, scope exists to improve the integration of the space-time variability of forcing terms in the drought assessment methodology.

In this context, this study aims to analyse and quantify the spatial and temporal dynamics of precipitation and drought in the Medjerda catchment of Tunisia. The study also aims to explore the contributing factors of the observed variability.

2 Methodology

2.1 Study area

The Medjerda is a transboundary river flowing from Algeria and to Tunisia (Fig. 1). The source of the Medjerda is located near Souk-Ahras, in the Algerian territory. From there onwards, the Medjerda and its tributaries collect the runoff water of part of the eastern Algerian and northern Tunisia territory (Gara et al. 2015). The total area of the Medjerda catchment is 23 700 km2. Almost one-third of this is located in Algeria. The river discharges nearly one billion cubic meters of water per year, which constitutes nearly 40% of the country’s wide available surface water.

Figure 1. Geographical location of the Medjerda catchment study site in Tunisia and the 41 selected rainfall stations.

The Medjerda catchment in Tunisia is divided into three distinct parts. The first part, referred to as the High Medjerda, initiates from the Algerian border to the city of Bousalem. The second part, the Middle Medjerda, includes all the tributaries and extends until Medjez el Bab. Finally, the Lower Medjerda ends at the river mouth of Kalaat Landalous. From a geographic point of view, the catchment belongs to the mountainous Tell region, characterized by cold and rainy winters, with a very heterogeneous hydrological regime. The climate of the Medjerda catchment belongs to the humid, sub-humid and semi-arid bioclimatic Mediterranean stages.

The relief in the Medjerda catchment is characterized by a large variability, having an altitude ranging between 0 and 1637 m. The slopes of the catchment are generally high in the upstream mountainous parts exceeding sometimes 40%. The slopes become weaker in the downstream part with values less than 1% in most parts of its lower valley.

2.2 Data collection and processing

The overall methodology of data processing methodology is summarized in Fig. 2. The rainfall database collected from the General Directory of Water Resources in Tunisia (DGRE) (http://www.agridata.tn/organization/dgre) was used. From 172 available rainfall stations, 41 stations were selected having monthly and annual observed rainfall data for 41 years, from 1972/73 to 2012/13. To obtain reliable time series, we only considered stations with missing value ratio less than 15% (Table 1). It was also the ambition to have a well-distributed spatial coverage of the catchment (Fig. 1).

Table 1. Descriptive statistics of annual rainfall for the raingauges in Medjerda catchment. Z: altitude (m a.s.l.), Min: minimum annual rainfall (mm/year), Max: maximum annual rainfall (mm/year), SD: standard deviation of annual rainfall (mm/year); CI: coefficient of irregularity; CV: coefficient of variation.

Station	Z (m)	Gaps (%)	Min (mm/year)	Max (mm/year)	Mean (mm/year)	Variance	SD (mm/year)	CI
Ain Beya Oued Ghezala	330	14.63	559.0	1458.3	961.5	55,363.9	235.3	2.6
Ain Debba	470	12.20	790.8	2214.3	1293.6	123,277.4	351.1	2.8
Ain Guesil 1	563	2.44	178.4	631.0	362.3	10,911.4	104.5	3.5
Ain Kerma 1	601	14.63	190.7	683.5	355.1	9536.5	97.7	3.6
Ain S’koum	990	14.63	174.9	513.2	353.8	7412.7	86.1	2.9
Ain Tabia	416	4.88	259.2	776.8	427.4	17,945.9	134.0	3.0
Ain Tounga SE	110	14.63	237.0	714.6	450.9	16,159.3	127.1	3.0
Ain Zeligua	853	14.63	222.4	644.1	412.2	10,113.6	100.6	2.9
Akouat Gare	350	0.00	244.7	722.0	422.4	12,576.6	112.1	3.0
Ben Metir 2 SM	525	12.20	230.9	684.2	428.4	13,982.3	118.2	3.0
Borj El Amri	55	14.63	638.6	1811.9	1173.9	91,980.0	303.3	2.8
Bou Salem DRE	138	7.32	241.5	799.1	483.9	21,043.0	145.1	3.3
Cherfech CRGR	59	2.44	258.7	698.6	430.6	13,848.5	117.7	2.7
Cite Mellegue SM	256	14.63	609.4	1704.3	1004.0	71,154.9	266.7	2.8
Dehmani Municipalite	622	4.88	242.6	789.7	420.3	14,243.3	119.3	3.3
Gardimaou DRE	195	9.76	257.8	722.7	454.3	14,943.1	122.2	2.8
Jantoura	390	14.63	240.8	718.4	456.9	12,095.7	110.0	3.0
Kef CMA	491	0.00	236.4	727.3	434.1	12,786.1	113.1	3.1
Kef. BIRH	620	2.44	315.1	814.6	534.8	17,315.1	131.6	2.6
Krib Ferme Cossem	447	9.76	201.4	599.7	403.5	11,617.2	107.8	3.0
Ksar Tyr Les Allobro	206	2.44	240.8	694.9	409.5	13,887.6	117.8	2.9
Ksour Ecole	720	14.63	353.8	775.3	522.1	14,560.6	120.7	2.2
Maktar PF	900	9.76	234.2	644.1	405.5	12,395.5	111.3	2.8
Mejez El Bab PF	142	4.88	144.1	479.6	312.1	5297.6	72.8	3.3
Oued Mellegue K 13	324	4.88	194.1	789.0	418.5	14,598.1	120.8	4.1
Oued Rmil	300	12.20	222.4	681.8	420.5	12,887.2	113.5	3.1
Oued Tine	406	14.63	311.9	975.8	622.9	31,417.8	177.3	3.1
Porto Farina Ghar Melah	10	4.88	264.8	868.9	532.4	20,527.1	143.3	3.3
Sakiet Sidi Youssef	803	4.88	227.3	679.7	408.0	10,300.1	101.5	3.0
Sers Agricole	501	12.20	220.0	665.9	404.1	11,002.5	104.9	3.0
Sers Delegation	501	14.63	304.8	806.9	505.8	18,216.3	135.0	2.6
Sidi Bou Rouis SM	412	2.44	247.3	799.3	526.4	19,034.6	138.0	3.2
Sidi Hamada	690	9.76	231.6	662.9	409.6	11,702.0	108.2	2.9
Sk El Arba (Jendouba) SM	143	14.63	246.9	670.2	410.2	10,765.8	103.8	2.7
Sk El Arba (Jendouba) SE	143	9.76	284.8	695.0	459.5	12,933.5	113.7	2.4
Slouguia	65	7.32	245.0	704.2	433.1	13,114.4	114.5	2.9
Tajerouine Ferme D’etat	511	7.32	155.9	552.4	343.7	9017.0	95.0	3.5
Teboursouk SM	440	0.00	343.8	918.6	580.3	18,854.1	137.3	2.7
Tibar SM	365	12.20	361.5	817.4	553.2	17,586.1	132.6	2.3
Zaafrane UCP	530	12.20	234.8	811.6	481.7	17,392.8	131.9	3.5
Zouarine Gare	571	14.63	151.3	463.8	260.1	5417.1	73.6	3.1

Figure 2. Flowchart of the methodology.

The homogeneity of rainfall series for long periods was assessed to ensure the quality of the data series and to identify the breaks of homogeneity. The Hydraccess software (Vauchel 2005) was used to obtain homogeneous and continuous time series for the 41 years and 41 stations. The Regional Vector Method (RVM) developed by ORSTOM/IRD was used to obtain the standardized index of monthly and annual rainfall for a homogeneous region (Taibi et al. 2017).

Hiez (1977) defined the regional vector (RV) as a time series of annual or monthly rainfall or hydrometric indices, extracted from the most likely information available from a set of observed stations grouped by region. The calculation of the RV is based on two hypotheses: (a) the pseudo-proportionality hypothesis, ensuring concomitant variations in rainfall at all stations; and (b) the hypothesis that most probable information can be achieved using a frequency criterion.

To have a continuous time series for each station, an index Z_i, the ratio of the annual rainfall to the average annual rainfall, was calculated. If Z_i is greater than one, the year is humid and if it is less than one then the year is dry (Brunet-Moret 1979, Sambou et al. 2018).

Since the RV has been assessed, the missing data of the individual stations can be estimated on a monthly as well as on an annual scale. The extension consists of reconstructing the rainfall values by multiplying the vector index of a year by the extended average of the station over the study period.

The Pettitt test was chosen to validate the performance of the RVM based homogenized and interpolated data (Pettitt 1979). This test is known for its robustness in detecting changes in stationarity, as demonstrated in previous studies (Kulkarni and Pardeshi 2019, Kundu and Mondal 2019, Li et al. 2019).

It consists of cutting the main series of N elements into two sub-series at each time t between 1 and N – 1. The main series has a singularity at time t if both sub-series (series (X_i), i = 1, t and (X_j), i = t + 1, N) have different distributions:

(1) $D_{ij}=\mathrm{s}\mathrm{g}\mathrm{n}\left(X_i-X_j\right)\left\{\begin{array}{c}1\mathrm{i}\mathrm{f}{X}_i-X_j>0\\ 0\mathrm{i}\mathrm{f}{X}_i-X_j=0\\ -1\mathrm{i}\mathrm{f}{X}_i-X_j<0\end{array}\right.$(1)

where sgn(X_j − X_i) is the sign function, and X_i and X_j are the data values in the time series i and j (j > i), respectively.

The variable U_t,_N is subsequently defined as:

(2) $U_{t,N}={\sum }_{i=1}^t{\sum }_{j=t+1}^N{D}_{ij}$(2)

The variable K_N is defined by the absolute maximum of U_t,_N for t from 1 to N – 1:

(3) $K_N=\left|\mathrm{m}\mathrm{a}\mathrm{x}\left(D_{ij}\right)\right|$(3)

The null hypothesis H₀ of non-failure is tested utilizing a non-parametric test. For a given risk value α, H₀ is rejected if a certain probability (Prob (K_N > k)) calculated from the studied time series is less than α.

2.3 Descriptive statistics

To characterize the spatial and temporal variability of rainfall in the Medjerda catchment, the descriptive statistics for annual rainfall for the 41 stations were used. Use was made of the variance (V), the standard deviation (SD), the coefficient of variation (CV), the interval of variation (IV) and the coefficient of irregularity (CI). These techniques are defined as follows:

(4) $V=\frac{1}{N}{\sum }_{i=1}^n{n}_i{x_i}^2-{\bar{x}}^2$(4)

where V is the variance, $\bar{x}$is the sample mean, x_i is the ith element from the sample, n_i is the number of the ith element from the sample and N is the number of elements in the sample.

(5) $\mathrm{S}\mathrm{D}=\sqrt{V}$(5)

(6) $\mathrm{C}\mathrm{V}=\frac{\sigma }{\bar{x}}$(6)

where $\sigma$ is the square of the sample variance and $\bar{x}$ is the sample mean.

(7) $\mathrm{I}\mathrm{V}=P_{\mathrm{m}\mathrm{a}\mathrm{x}}-P_{\mathrm{m}\mathrm{i}\mathrm{n}}$(7)

where P_max and P_min are, respectively, the maximum and minimum annual precipitation (mm/year).

(8) $\mathrm{C}\mathrm{I}=\frac{P_{\mathrm{m}\mathrm{a}\mathrm{x}}}{P_{\mathrm{m}\mathrm{i}\mathrm{n}}}$(8)

2.4 Trend testing

2.4.1 Mann-Kendall test

The Mann-Kendall trend test was presented by Mann (1945) and then further developed by Kendall (1975). It is a non-parametric test generally used to identify the trend in meteorological time series data (Abdulla and Al-Qadami 2019, Panda and Sahu 2019). The null hypothesis H₀ of this test is that there is no trend. The alternative hypothesis means that there is a significant trend. The Mann-Kendall test is based on the calculation of Kendall’s tau (τ) coefficient measuring the association between two variables. More specifically, it measures the correlation of rank between the two variables. In the case of the trend test, the first series is a growing time indicator with ever-increasing ranks, which simplifies the computation process. Kendall’s S statistic and its variance are given by:

(9) $S={\sum }_{i=1}^{n-1}{\sum }_{j=i+1}^n\mathrm{s}\mathrm{g}\mathrm{n}\left(x_j-x_i\right)$(9)

(10) $\mathrm{v}\mathrm{a}\mathrm{r}\left(S\right)=\frac{n\left(n-1\right)\left(2n+5\right)}{18}$(10)

where n is the number of data in the series, x_i and x_j (j > i) are the observations which are supposed to be independent and sgn(x_j − x_i) is the sign function, as defined in eq (1)(1) $D_{ij}=\mathrm{s}\mathrm{g}\mathrm{n}\left(X_i-X_j\right)\left\{\begin{array}{c}1\mathrm{i}\mathrm{f}{X}_i-X_j>0\\ 0\mathrm{i}\mathrm{f}{X}_i-X_j=0\\ -1\mathrm{i}\mathrm{f}{X}_i-X_j<0\end{array}\right.$(1) .

The independence hypothesis is verified by the Wald-Wolfowitz test, which is a non-parametric test that tests that the elements of the sequence are mutually independent. It is used to check whether there is a significant autocorrelation between the observations.

2.4.2 Sen’s slope estimator test

Sen’s slope estimator test is a non-parametric method for estimating the slope of the trend in the sample of N pairs of data (Sen 1968). The method is based on the calculation of Q_i, given by:

(11) $Q_i=\frac{x_j-x_k}{j-k}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{i}=1,\dots ,\mathrm{N}$(11)

where x_j and x_k are the data values at times j and k (j > k), respectively.

The N values of Q_i are ranked from the smallest to the largest. The slope median defines the Sen slope estimator and is calculated as follows:

(12) $Q_{\mathrm{m}\mathrm{e}\mathrm{d}}=\left\{\begin{array}{c}Q\left[\left(N+1\right)/2\right]\mathrm{i}\mathrm{f}N\mathrm{i}\mathrm{s}\mathrm{o}\mathrm{d}\mathrm{d}\\ \frac{Q\left[N/2\right]+Q\left[\left(N+2\right)/2\right]}{2}\mathrm{i}\mathrm{f}N\mathrm{i}\mathrm{s}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{n}\end{array}\right.$(12)

Sen’s slope estimator, which expresses the trend and slope magnitude in the entire time series, is used in several studies for its robustness for estimating changes (Atta-ur-Rahman and Dawood 2017, Hui-Mean et al. 2018).

2.5 Drought indices

Two climatic indices were applied to study the variability of the rainfall regime in the Medjerda catchment: the SPI used for the study of droughts (Louvet et al. 2011) and the AI, used to characterize the complex phenomenon of aridity (Sahin 2012).

The SPI was presented first by McKee et al. (1993) to characterize precipitation deficits. It reflects the effect of drought on the availability of water resources. The World Meteorological Organization recommends the SPI for dry weather monitoring as it only requires precipitation for calculations (Santos et al. 2017, 2019). This index is calculated as follows:

(13) $\mathrm{S}\mathrm{P}\mathrm{I}=\frac{P_i-P_i}{\mathrm{S}\mathrm{D}}$(13)

where SD is the standard deviation of annual precipitation; P_i is the annual precipitation of a given year i (mm/year); and P_m is the mean annual precipitation for a station (mm/year).

Many studies have used the SPI for the analysis of the spatio-temporal variability of drought (Caloiero and Veltri 2019, Dhurmea et al. 2019, Gara et al. 2019, Ouatiki et al. 2019). The SPI is often used as a key in drought classification schemes (Table 2).

Table 2. Standardized precipitation index (SPI) and aridity index (AI) classification.

SPI	Degree of drought	AI	Type of climate
SPI > 2	Extreme humid	0 < AI < 5	Hyper-arid
1 < SPI < 2	Severe humid	5 < AI < 7.5	Arid
0 < SPI < 1	Moderate humid	7.5 < AI < 10	Steppe
–1 < SPI < 0	Moderate dry	10 < AI < 20	Semi-arid
–2 < SPI < – 1	Severe dry	20 < AI < 30	Sub-humid
SPI < – 2	Extreme dry	AI < 30	Humid

The AI is defined as the ratio of precipitation to average temperatures according to the equation (Martonne 1926):

(14) $\mathrm{A}\mathrm{I}=\frac{P}{T+10}$(14)

where P the annual average rainfall (mm/year) and T is the annual average temperature (°C).

The annual average temperature data for the observed period were provided by the National Institute of Meteorology (INM). The AI is also often used as a key in climate classification schemes (Table 2).

2.6 Spatial data treatment

Geostatisticss allow analysis of the spatial and temporal structure of hydrometeorological and climatic variables (Slimani et al. 2007). Geostatistics is built on the theory of random variables and is used to study the spatial variability of a quantitative variable defined at any point in space. Geostatistics proves to be very suitable for mapping natural phenomena (Webster and Oliver 2007, Tripathi et al. 2015). The semi-variogram is used with ordinary kriging (OK) to predict the random variable for unobserved locations (Varouchakis et al. 2018). In this study, a semi-variogram analysis and OK were used to study and map the annual mean precipitation.

The cross-validation method was used to evaluate the robustness of the spatial prediction model (Gerstmann et al. 2016). The cross-validation results allow the calculation of the following model’s goodness-of-fit parameters: the mean standardized error (MSE), which should be close to zero; and the root mean square standardized error (RMSSE), which should be close to one. These parameters were calculated as follows:

(15) $\mathrm{M}\mathrm{S}\mathrm{E}=\frac{{\sum }_{i=1}^n\left(\hat{z}\left(s_i\right)-z\left(s_i\right)\right)/\hat{\sigma }\left(s_i\right)}{n}$(15)

(16) $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{S}\mathrm{E}=\sqrt{\frac{{\sum }_{i=1}^n{\left[\left(\hat{z}\left(s_i\right)-z\left(s_i\right)\right)/\hat{\sigma }\left(s_i\right)\right]}^2}{n}}$(16)

where $\hat{z}\left(s_i\right)$ are the measured values and $z\left(s_i\right)$ are the predicted values at time/place i.

2.7 Data reduction and station classification

Drought clustering was performed using the K-means clustering technique. This technique was introduced by MacQueen (1967) and consists of a hierarchical classification approach commonly used to define climate divisions (Li et al. 2015).

To identify a station class with a similar rainfall pattern, the K-means clustering method was performed by considering the selected stations as quantitative variables and years as individuals.

To choose the number of clusters, the K-means was launched with different class values. The method allowed calculating the variance of the different clusters. The optimal number of classes was selected from the minimization of the intra-class distance. The classification that was used to obtain an intra-class variance minimization solution is the Wilks lambda criterion. This method is sensitive to scale effects.

Drought data were further analysed with principle component analysis (PCA). When the rainfall regime represents a large variability, the PCA can serve as a tool for analysing and understanding the influence of different factors.

A PCA based on a Spearman correlation matrix was performed by admitting the selected stations as individuals. The following attributes were considered as explicative variables: the latitude of the station (UTM coordinate), the longitude (UTM coordinate), the altitude (m a.s.l.), the distance to the sea (km), the average monthly precipitation (mm/month), the average annual precipitation (mm/year), the SPI average; and the AI average. The explicative variables were projected in the factorial plan allowing to highlight their correlations.

To determine the optimal number of factors to be retained for the interpretation of the results, the cumulative percentage of variability represented by the factorial axes is considered. The method further looks for the presence of an inflection point when representing the factorial axes with the accumulation of eigenvalues.

3 Results and discussion

3.1 Overview of the used data

The RVM was applied to correct and complete the original precipitation time series. The results of the Pettitt test on the corrected and completed precipitation time series show that the 41 series of annual rainfall are homogeneous. This indicates that the reconstructed series can be used to study the spatio-temporal variability of annual rainfall.

The series of the annual index Z_i is shown in Fig. 3. This index illustrates the overall annual variability of rainfall in the Medjerda catchment. 2002/03 appears as a very rainy season, while 1993/94 is the driest season. This result is consistent with several previous studies for other Tunisian regions (Aouissi et al. 2018, Jemai et al. 2018).

Figure 3. Temporal variation of Z_i (i.e. the ratio of the regional value of the annual rainfall to the average annual rainfall) of the precipitation in the Medjerda catchment.

The SD and CV of the annual rainfall events indicate high inter-annual variability of rainfall in the Medjerda catchment (Table 1). The SD is always larger than 72 mm/year and sometimes reaches more than 350 mm/year. The average annual rainfall of the Medjerda catchment over the selected period is 490 mm. The highest averages are observed at the stations located in the northwestern part of the study area. The CI of the annual rainfall is, on average, equal to 3. It ranges from 2.2 at Maktar PF to more than 4.1 at Oued Rmil. This reflects a spatial and temporal irregularity of the rainfall regime in this region. This result is consistent with other rainfall analysis studies carried out in catchments close to the Medjerda catchment (Ennajah et al. 2015, Jemai et al. 2018).

The analysis of the ratio of extreme values by the annual average for each station confirms the high temporal variability of annual precipitation. During wet years, these ratios can exceed 190%. On the other hand, this ratio may not exceed 50% for dry years. The ratio of the annual average precipitation at the level of each station to the annual average precipitation of all 41 stations illustrates the important spatial variability of annual rainfall. This ratio is 51% at Zouarine Gare and 256% at Ain Debba.

The box plots of the mean annual precipitation over the 41 raingauges for the period from 1972/73 to 2012/13 are given in Fig. 4. The box plot confirms the important spatial and temporal variability of the annual precipitation in the catchment. This important space-time variability of annual precipitation is consistent with observations in the Merguellil catchment in Tunisia (Chargui et al. 2018).

Figure 4. Boxplots of the annual precipitation in the Medjerda catchment from 1972/73 to 2012/13. The boxplots were created with sigmaplot 11.0. The upper and lower limits represent the 75th and 25th percentiles and the horizontal bars at the ends of the lines outside the boxes represent the 90th and 10th percentiles. The horizontal line in the box represents the median. Black circles represent measures that fall outside the 90th and 10th percentile boundaries.

The time course of the mean annual precipitation displays exceptional humid seasons such as 1972/73, 1973/74, 1982/83 and 1985/86. During these seasons, annual precipitation exceeds 1000 mm/year.

3.2 Trend analysis

The Mann-Kendall test and Sen’s slope estimator test were applied to detect possible temporal trends in the annual rainfall series for all the 41 rainfall stations (Table 3). The null hypothesis H₀ of the Mann-Kendall test (i.e. there is no trend) has been accepted if the p value exceeds the confidence level α = 0.05). Results show that there is no significant trend in annual rainfall for the majority of the selected stations. Significant trends are only observed at the Akouat Gare and Dehmani municipalite stations, with trend magnitudes of 3.2 and 2.8 mm/year, respectively. The absence of significant trends in annual rainfall is consistent with the conclusions drawn from many other studies showing that in Tunisia rainfall trends are not always detected (Chargui et al. 2018, Feki et al. 2018).

Table 3. Results of the MK test and the Sen slope estimator test of annual rainfall for the raingauges in the Medjerda catchment.

Station	Kendall’s tau, τ	S	P value	Sen’s slope
Ain Beya Oued Ghezala	0.066	54	0.554	2.345
Ain Debba	0.046	38	0.680	2.721
Ain Guesil 1	0.157	129	0.150	1.602
Ain Kerma 1	0.111	91	0.312	0.8
Ain S’koum	0.178	146	0.104	1.639
Ain Tabia	0.093	76	0.402	1.479
Ain Tounga SE	0.039	32	0.729	1.157
Ain Zeligua	0.163	134	0.136	2.003
Akouat Gare	0.268	220	0.013*	3.224
Ben Metir 2 SM	0.131	107	0.234	2.25
Borj El Amri	0.090	74	0.415	5.659
Bou Salem DRE	0.068	56	0.539	1.353
Cherfech CRGR	0.083	68	0.454	1.053
Cite Mellegue SM	0.090	74	0.415	2.789
Dehmani Municipalite	0.229	188	0.035*	2.782
Gardimaou DRE	0.149	122	0.175	2.474
Jantoura	0.146	120	0.183	2.406
Kef CMA	0.153	125	0.164	1.906
Kef. BIRH	0.117	96	0.288	2.394
Krib Ferme Cossem	0.171	140	0.119	2.48
Ksar Tyr Les Allobro	0.129	106	0.240	2.047
Ksour Ecole	0.146	120	0.183	2.398
Maktar PF	0.140	115	0.200	1.929
Mejez El Bab PF	0.151	124	0.168	1.251
Oued Mellegue K 13	0.166	136	0.130	2.517
Oued Rmil	0.146	120	0.183	2.068
Oued Tine	0.149	122	0.175	3.306
Porto Farina Ghar Melah	0.180	148	0.099	3.261
Sakiet Sidi Youssef	0.173	142	0.114	2.252
Sers Agricole	0.207	170	0.057	2.835
Sers Delegation	0.117	96	0.288	2.561
Sidi Bou Rouis SM	0.088	72	0.428	1.852
Sidi Hamada	0.120	98	0.278	1.679
Sk El Arba (Jendouba) SM	0.171	140	0.119	2.246
Sk El Arba (Jendouba) SE	0.137	112	0.214	2.087
Slouguia	0.111	91	0.312	2.114
Tajerouine Ferme D’etat	0.105	86	0.342	1.234
Teboursouk SM	0.056	46	0.616	0.942
Tibar SM	0.200	164	0.067	3.236
Zaafrane UCP	0.093	76	0.402	1.671
Zouarine Gare	0.071	58	0.525	0.51

*: trend significant according to 5% significance level.

3.3 Geostatistical analysis

The variogram of annual rainfall in the Medjerda catchment is modeled employing a Gaussian model (Fig. 5 and Table 4). The variogram model was used to predict rainfall using OK. The quality of the variogram model was thereby evaluated using a cross-validation goodness-of-fit of predicted annual precipitation. The RMSSE of predicted annual precipitation was close to one and the MSE was close to zero, suggesting that the variogram model is appropriate. The graphical representation of the observed versus predicted annual precipitation confirms this goodness-of-fit (Fig. 6).

Table 4. Omnidirectional variogram parameters of the annual precipitation applied to raingauges of the Medjerda catchment.

Model type	Parameter			Cross-validation
	Nugget (mm^2/year^2)	Sill (mm^2/year^2)	Range (m)	MSE (mm/year)	RMSSE (mm/year)
Gaussian	3300	35,000	83,290	–0.0133	1.0412

Nugget: The height of the jump of the semi-variogram at the discontinuity at the origin.

Sill: Limit of the variogram tending when the lag tends to infinity.

Range: The distance in which the difference of the variogram from the sill becomes negligible.

The Gaussian variogram model: γ(h) = C [1 – exp(–3(h/a)2)], with γ(h): the semi-variogram; C: sill; a: range; h: lag distance.

Figure 5. Omnidirectional variogram of inter-annual rainfall in the Medjerda catchment.

Figure 6. Graphical goodness-of-fit of the observed and predicted rainfall of the 41 raingauges in the Medjerda catchment.

Variographic analysis shows that the variance of rainfall increases to a certain distance. This spatial dependence may be due to the relief in the catchment. Indeed, the study catchment has an accentuated topography in the upper valley. This is consistent with several studies showing that the spatio-temporal variability of rainfall in Tunisia is influenced by the relief (Slimani et al. 2007, Gader et al. 2015). Yet, altitude is not the only factor that influences this variability (Feki et al. 2012).

The results of the variography were used for rainfall mapping using OK (Fig. 7). The mean annual rainfall map confirms that the rainfall in the Medjerda catchment is characterized by high spatial variability, varying between 300 and 1200 mm/year. It is found that rainfalls are higher than 1000 mm/year in the northwest and gradually decrease when moving to the southeastern part of the catchment. This high spatial variability of precipitation is typical for the Mediterranean regions (Rodrigo Comino et al. 2016, Rysman et al. 2016, Cid et al. 2017).

Figure 7. Mean annual rainfall map of the Medjerda catchment.

3.4 Drought indices

Mapping of the mean inter-annual SPI, calculated per decade over the period 1973–2012, reflects the temporal dynamics of precipitation (Fig. 8). The 1973–1982 decade is characterized by high temporal variation of the rainfall regime varying from extreme humid to extreme dry. The 1983–1992 and 1993–2002 decades are characterized by a moderate drought over the majority of the catchment except in some locations in the southwest.

Figure 8. Spatial distribution of inter-annual SPI during four decades in the Medjerda catchment.

The regionalization of the average values of the SPI over different periods highlights the large spatio-temporal variability of the rainfall in the Medjerda catchment. During the last decade of 2003–2012, an increase in the rainfall occurred. The rainfall regime became humid all over the catchment. This decade is characterized by a succession of floods in the Medjerda catchment such as the floods of 2003, 2006, 2009, 2011 and 2012 (Ben Khalfallah and Saidi 2018). The period 2000–2007 presented four floods whereas the period 1960–2000 had an average flood frequency of 1.75 events per 10 years (Jebari et al. 2012). Hence, the hydrological regime is subjected to considerable temporal alteration at the decadal scale, even if no trends were observed at the multi-decadal scale. The observed hydrologic alteration could be explained by climate change, which is well described for Mediterranean catchments (Zielhofer et al. 2015).

Using the AI, we observe that the Medjerda catchment can be decomposed into three climate zones from north to south (Abid et al. 2018) (Fig. 9). The AI ranges from 7 to 49 and is lower in the southwest and higher in the north of the Medjerda catchment. The mapping of the AI over four decades shows that the majority of the catchment situates in the semi-arid bioclimatic stage. The northern part of the catchment belongs to the sub-humid and humid region. This part of the catchment exhibits a larger temporal variability.

Figure 9. Spatial distribution of inter-annual AI during four decades in the Medjerda catchment.

The bioclimatic stratification of the study area, as illustrated in the AI maps, is largely controlled by the topography. The catchment exhibits a transition between the northern region with high altitudes characterized by a humid and sub-humid climate to the southern region with medium and low altitudes characterized by a semi-arid to arid climate.

3.5 Multivariate statistical analysis

We further explored the spatial dynamics of the mean annual rainfall through the multivariate cluster and ACP analysis. The K-means method on the annual rainfall was carried out for different class numbers ranging from 3 to 6. The curve of the interclass distance according to the class number, showed an inflection point in class 4, corresponding to the optimal number of classes to be considered. The identified clusters are shown in Fig. 10. The clusters can be described as follows:

Figure 10. Clusters of rainfall stations differentiated by means of the K-means clustering method on the mean annual rainfall of the Medjerda catchment.

3.5.1 Moderate dry and humid stations

This class encompasses the following stations: Ain Debba, Jantoura, Ben Metir 2 SM, Ain Beya Oued Rhezalla and Porto Farina Ghar Meleh. These gauges are located in the northern part of the Medjerda catchment relatively close to the sea. They are characterized by a moderate drought to humid. This group is influenced by the north-west, north and north-eastern rainfall, which is reinforced by the presence of the sea.

3.5.2 Moderate dry and semi-arid stations

This group encompasses the following gauges: Dehmani municipality, Ksour Ecole, Ksar Tyr the Allobro, Medjez El Bab PF, Oued Rmil, Oued Tine, Sers Agricole, Sers Delegation, Slouguia, Sk El Arba (Jendouba) SE, and Bou Salem DRE. The behavior of these stations is close to that of the previous stations, but they exhibit less precipitation and are situated more southeast in the Medjerda catchment. This cluster is characterized by a large variation of altitude, decreasing from 720 m a.s.l. at Ksour Ecole to 65 m a.s.l. at the Slouguia stations.

3.5.3 Hilly and mountainous stations

The third cluster contains Ain Guesil 1, Ain Kerma 1, Ain S’koum, Ain Zeligua, Akouat Gare, Oued Mellegue K 13, Tajerouine Ferme d’Etat and Zouarine Gare. These stations are located in the hilly and mountainous region of the Tells Mountains. The altitudes of these stations exceed 400 m a.s.l. and they are influenced by continental characteristics. The precipitations range varies between 260 and 422 mm/year, respectively, at Zouarine Gare and Akouat Gare.

3.5.4 Stations in the Medjerda Valley

This cluster includes the following stations, located in the Medjerda Valley and within its main tributaries: Ain Tabia, Ain Tounga, Borj El Amri, Cherfech CRGR, Cite Mellegue SM, Gardimaou DRE, Kef BIRH, Kef CMA, Krib Ferme Cossem, Maktar PF, Sakiet Sidi Youssef, Sidi Bou Rouis SM, Sidi Hmada, SK El Arba (Jendouba) SM, Teboursouk SM, Tibar SM and Zaafrane UCP. This cluster encompasses the largest group of gauges with similar behaviour. The gauges of this cluster have an altitude varying between 55 and 900 m a.s.l. and are characterized by an average annual rainfall between 420 and 580 mm/year. Moreover, there has been a decrease in the rainfall of the north towards the south, which is explained by the presence of the Khemirs and Mogods chains bordering the Medjerda Valley in the north and the Tells Mountains in the south.

To explore the relationships between mean precipitation, drought indices, and possible other explaining factors, a PCA analysis was performed. The PCA allows summarizing the correlation between average monthly precipitation, mean inter-annual precipitation, climate indices, position in the catchment and altitude. Cluster analysis combined with PCA analysis has proven to be effective in analyzing rainfall distribution patterns according to topographic and other climatic parameters (Rau et al. 2017). The first two components of the PCA explain 70% of the total variance. The projection of the variables within the biplot composed of the first and second principal component (PC1 and PC2) is given in Fig. 11. It is observed that inter-annual rainfall is strongly influenced by the wet seasons (autumn and winter) and that it is strongly correlated with latitude but weakly with longitude. These results are in agreement with similar studies performed on Mediterranean catchments. It illustrates the important seasonality of rainfall (Deitch et al. 2017, Corona et al. 2018). The precipitation of the summer months is the most correlated with the altitude and the proximity of the sea. The first principal component, PC1, is mainly explained by monthly winter precipitations, SPI and the distance to the sea, while PC2 is determined by the monthly summer precipitation, AI, and latitude. A strong negative correlation is observed between the SPI and the rainy months, and between AI and the summer months. The SPI should decrease strongly during years characterized by winter precipitation shortage, however, it increases when precipitation is high during rainy months (Fischer et al. 2019); this could explain the negative correlation. The aridity index is related to high temperature, mainly in summer, and to less precipitation during dry years (Zhao et al. 2019).

Figure 11. Projection of the selected variables in the biplot composed of the first (PC1) and second principal components (PC2). X_UTM, Y_UTM and Z (m): coordinates of the selected station in the Universal Transverse Mercator coordinate system, Zone 32 Carthage North; SPI: annual mean SPI; AI: annual mean AI; P_m: mean inter-annual rainfall (mm); Sea: distance to the sea (km); January, …, December: average monthly rainfall (mm).

These results are consistent with those found by Ellouze et al. (2009), who showed that the SPI is explained by annual rainfall variability, topography and seasonality in southern Tunisia.

The representation of gauges in the factorial plan shows the dispersion of the stations around the two axes (Fig. 12). The stations are labelled according to the clusters identified by the K-means method. The stations belonging to the moderate dry and humid stations have the greatest inertia with the first axis while that of the hilly and mountainous stations have the lowest values of inertia with this axis.

Figure 12. Projection of the rainfall stations in the biplot composed of PC1 and PC2. The rainfall stations are labelled according to the clusters identified by the K-means clustering algorithm.

4 Conclusion

Drought trends and variability in precipitation patterns were identified for the Medjerda catchment in Tunisia. We used data from 41 ground-truth precipitation gauges collected for 41 years that were first pre-processed using the RVM method. We also explored factors controlling the spatio-temporal dynamics of drought and rainfall. The pre-processing of the precipitation data utilizing the RVM allowed generating homogeneous time series that could be further used to assess droughts, drought trends, and variability. The analysis of rainfall data and derived drought indices revealed remarkable spatio-temporal variability of drought. The analysis revealed also the occurrence of significant drought in major parts of the catchment. The variographic analysis showed that abounded Gaussian variogram model could be used to fit the annual rainfall semi-variogram. The variogram could subsequently be used to map the annual precipitation. A cluster analysis on the data of the 41 stations allowed identifying four distinct clusters that partially could be discriminated based on geographical attributes such as the geographical position, the altitude, and the distance to the sea. These results were also consistent with the results of the principal component analysis.

The obtained results make it possible to improve the knowledge of rainfall and drought variability in this specific catchment and the Mediterranean region more generally. The adopted methodology allows a better understanding of the space–time dynamics of rainfall and droughts and the factors that control this variability. Understanding the space–time dynamics is a prerequisite for hydrological risk management and mitigation in a climate-changing context.

Acknowledgements

The authors would like to thank the Tunisian General Directorate of Water Resources (DGRE) and the National Institute of Meteorology (INM) for providing the data. We gratefully acknowledge the editors and the anonymous reviewers of this paper, whose comments helped to improve it.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

The authors are grateful to the Laboratory of Water Science and Technology (STE) of the National Agronomic Institute of Tunisia (INAT) and to Sustainable Management of Water and Soil Resources (UR17AGR03) of High School of Engineering of Medjez el Bab (ESIM) – research programme 1. We received funding from the project “Adaptation of water resources management in the Medjerda watershed to the challenges of climate change” through Project 11 – Axis 2 of the Wallonia Brussels International and Tunisia Joint Commission 2016–18. Part of this research study was accomplished in GERU, UCLouvain with the funding of mobility programme provided by the Ministry of Higher Education and Scientific Research in Tunisia; Sustainable Management of Water and Soil Resources (UR17AGR03) of High School of Engineering of Medjez el Bab (ESIM) – research program 1; the project ‘Adaptation of water resources management in the Medjerda watershed to the challenges of climate change’through the project 11 – axis 2 of the Wallonia Brussels International and Tunisia Joint Commission 2016–18.; mobility program provided by the Ministry of Higher Education and Scientific Research in Tunisia.