2.1 Study area

The study region is part of a vast geographic unit, the steppe area, located to the west of the Algerian highlands. It is a semi-arid zone positioned between the edges of the Tell Atlas (Tlemcen and Ouarsenis mountains), to the north, and the Saharan Atlas (Ksours mountains), to the south (Fig. 1). Its coordinates are 1º 30’ E to 4º 0’ W, and 32° 20’ S to 36° 0’ N, covering a surface of 138 500 km². It is characterized by a wide endoreic expanse where discharges converge on the salt lakes lined in strings and the stream system is not well developed. Most of the wadies originate on the Tell Atlas crests and flow to the south into the Chergui salt lake.

Fig. 1 Location of the study region and rain gauging network. 

The area fits in the arid and moderate bioclimatic stage with cool winter (^{Emberger, 1942}). It constitutes a buffer zone between western coastal Algeria and Saharan western Algeria. Rainfall is not very abundant, but it often occur as violent storms on a regular basis during nine months of the year. The pluviometric mean is 318 mm spread over 47 days (^{Benslimane et al., 2015}).

Rainfall recorded on all pre-desert subregions of the western steppes is the weakest in the western Algeria region. An isohyets map drawn by the National Water Resources Agency (ANRH) shows that rainfall in this area varies from 100 to more than 300 mm yr^-1. The highest precipitation occurs in the mountainous area located to the north aside the Tlemcen and Saida mountains. Annual accumulation ranges from 200 to 300 mm on the high plateaus and the south, and up to more than 400 mm on the reliefs. The areas which received less rainfall (100 to 200 mm) are those ubicated in the salt lakes. Gross time series analysis shows that the driest spells were recorded in 1983 at the Slim station (47.2 mm), to the east; in 1985 at the Ain Skhouna station (77.7 mm), located northeast of the Chergui salt lake, and in 2004 in the north faces of the Saharan Atlas, at the Mecheria station (116 mm). The most humid years were 1971 and 1972 with about 312 mm of precipitation at the Slim station and 292 mm at the Ain Mehdi station. The most humid year recorded at the Mecheria station was 2008 (430 mm).

2.2 Data

The selection of stations was conducted based on the need of necessary information for the study regarding the length and spatial distribution of observational series of these stations within the study region. Concerning the altitude, the highest stations are those of Ain El Orak and Kherba Ouled Hellal, located at 630 masl. The rainfall gauge network of the steppe region of western Algeria includes more than 150 stations managed by the ANRH. These stations are very unequally spread from north to south and west to east on the region. The highest density of the network is found in the north; the density in the south is lower, and the salt lakes region is practically devoid of stations. For the study, 65 rainfall gauge stations that have more than 20 yrs. of observations and short gaps were selected (Fig.1). The basic characteristics of these rainfall gauge stations are presented in Table I.

2.3 L-moments

To validate the homogeneity of a region in terms of L-moments relationships, the discordance test proposed by ^{Hosking and Wallis (1993)} will be used:

ui= t(i),t3(i),t4(i)T

be u the vector containing the values t, t ₃ and t ₄ for the site i, where the exponent T refers to a vector or a matrix, and

u-= N-1∑i=1nui  the regional unweighted mean of L-moment rates for each group. The discordance measurement for the site is the. defined as follows:

In the equation below N is the sample size of each group and S ^-1 the invert of the matrix S.

^{Hosking and Wallis (1997)} proposed the criterion D _i ≥ 3 to exclude a station from the homogenous region. Inter-site relationship is used to identify the homogenous sites with a similar frequency distribution. The heterogeneity test H _i compares samples of L-moments relationships with the parameters of the Kappa distribution; it measures the heterogeneity between the sites of the same region. ^{Hosking and Wallis (1997)} proposed the following statistic:

where μ _v etσ _v are the mean and the standard deviation of N _sim of the simulated values of V ₁ .

V ₁ is calculated by the following relationship:

Simulations are done by using a flexible distribution with the regional mean of L-moments ratios 1, τ, τ ₃ , t and τ ₄ . Following ^{Hosking and Wallis (1993)}, ¹⁹⁹⁷), we used the distribution with four Kappa parameters with the following quantiles function in the simulations (^{Hosking, 1994}):

In order to get reliable values of μ _v and σ _v , the number of simulations N _Sim will be great. The region is considered as acceptably homogenous if H < 1, possibly heterogeneous if 1 ≤ H <2, and certainly heterogeneous if H ≥ 2. H1 is the homogeneity measurement in terms of L-CV; H ₂ is the homogeneity measurement in terms of L-CS, and H ₃ is the homogeneity measurement in terms of L-CK.

The statistic Z (^{Hosking and Wallis, 1991}) determines to which extent the L-skewness and L-kurtosis simulation of a fitted distribution correspond to the regional mean of L-skewness and L-kurtosis values, obtained from the data observed. In this work, the quality criterion for fitting is defined by the statistic Z^DIST depending on the various candidate distributions:

where t ₄ ^R is the mean value of t ₄ relating to the data of the relevant region, and B ₄ and are the bias and standard deviation of t ₄ , respectively, and they are defined as follows:

where N _Sim is the number of regional data simulations fixed generated by using a Kappa distribution (^{Hosking and Wallis 1988}), and m is the m ^th region simulated by obtaining a Kappa distribution.

The fitting is adequate if Z^DIST is sufficiently close to zero. |Z ^DIST | ≤ 1.64 If it can be said that the fitting is reasonable.

If more than a candidate distribution is acceptable, the one with the smallest |Z^DIST| is considered as the most appropriate. Furthermore, the L-moments ratio diagram is also used to identify the best distribution by comparing its proximity to the combination of L-skewness and L-kurtosis of the L-moments ratios.

The idea behind the use of the L-moments diagram is based on the operating of unique combinations of skewness and excess coefficients, to graphically identify the function that is closest to the studied sample distribution. To choose the most pertinent probability distribution, the use of existing regional works, or the L-moments diagram, is recommended (^{Ben-Zvi and Azmon, 1997}; ^{Chen et al., 2006}). When the scatter of L-moments ratios is close to a large number of probability distributions, the mean of the L-moments ratios of observations series is placed on the L-moments diagram and the closest distribution is chosen as the best distribution (^{Kumar et al., 2003}; ^{Chen et al., 2006)}.

The candidate distributions which will be used to this effect are:

General logistic (GLO), particularly the case of Kappa distribution with the parameter h = -1.

Generalized extreme value (GEV), particularly the case of Kappa distribution with the parameter h = 0.

Generalized normal (GNO).

Gaucho, particularly the case of Kappa distribution with the parameter h = 0.5.

Generalized pareto (AMP), particularly the case of Kappa distribution with the parameter h = 1.Pearson 3 (P3).

Kappa distribution (KAP).

The condition to have a homogenous region is that all sites may be described by a probability distribution that has common distribution parameters; thus, all sites of a homogenous region have a common regional curve of frequency, which is called regional growth curve (^{L-RAP, 2011}).

After determining the most appropriate frequency distribution model for each of the three homogenous groups, the quantiles for different return periods are estimated by using the rain index method. This procedure supposes that data of maximum daily rainfalls of different sites in a homogeneous group have the same statistic distribution, except for a scale parameter specific to a site or an index factor (^{Dalrymple, 1960}). The scale factor is considered as the rain index. The quantile of a homogenous group is estimated by the following equation:

where P _i (F) is the value of the daily maximum rainfalls at station i with a non-exceedance probability F; μ ̂ _i is the sample mean to this station, and p(F) is the dimensionless quantile with exeedance probability given by F. The totality (F) values for 0 < F < 1 result from the regional growth curve. This approach was used in many countries following the example of ^{Malekinezhad and Zare-Garizi (2014)} in Iran, ^{Yurekli et al. (2009)}, in Turkey, ^{Ngongondo et al. (2011)} in Malawi, and ^{Martin (2015)} in Canada. This approach is based upon the flood-index method (^{Darrlymple, 1960}). The priori assumption that data are independent and identically distributed as per the same statistic law should be made (^{St-Hilaire et al., 2003}). The same flood index is commonly used to develop the regional models of frequency for the sites where hydrologic information is not sufficient for reliable information of extreme events (^{Cunane, 1987}; ^{Watt, 1989}; ^{Portela and Dias, 2005}; ^{Saf 2009}; ^{Nyeko-Origami et al., 2012}).

Uncertainty is detailed after quantiles estimation for a station in each group through regional analysis (Eq. 4). Uncertainty is evaluated by using the bias (BIAS) and the square root of the root mean square error (RMSE) for each homogenous group:

Pmax _i ^R and Pmax _i ^L represent the quantiles of the return period T estimated for the station i, respectively, for the GEV distribution regional and local parameters. N is the number of stations of each homogenous group.

Both MAE and RMSE express the average model prediction error in units of the variable of interest. Both metrics can range from 0 to ∞ and are indifferent to the direction of errors; they are negatively-oriented scores, which means lower values are better. The more that BIAS and RMSE values come closer to 0, both the precision of the value estimated and the method are better.

L-moments present weaknesses and strengths (^{El Adlouni et al., 2003}):

The robustness of L-moments faced with large values may exclude the information regarding extreme values (^{Bernier, 1993}). According to ^{Klemes (2000)}, this technique favors the choice for the GEV distribution. ^{Ben-Zvi and Azmon (1997)} maintain that the L-moments diagram does not provide the best fitting distribution among the acceptable ones.

On the other hand, sample variability affects L-moments less (^{Vogel and Fennessey, 1993}) and they are more robust in the presence of outsiders (^{El Adlouni et al., 2003}).

2.4. Ascendant hierarchical classification

The ascendant hierarchical classification is a mode of classification that consists in aggregating first the most similar individuals, then slightly less similar observations or groups of observations, and so on until the trivial grouping of the whole of the sample.

Let there be a set of n individuals characterized by p variables X ₁ ,X ₂ ,X ₃ ,……….,X _p that we want to group in k classes or subclasses as homogeneous as possible according to defined criteria. To do this, according to each criterion, the closest individuals that are connected to each other must be searched. Gradually, all the individuals following a hierarchical tree or dendrogram are grouped together. This dendrogram presents the composition of the different classes as well as the order in which they are formed.

Individuals in the dendrogram are organized hierarchically according to the distances that separate them. The type of distance between individuals is selected depending on the data studied. Among others, the Euclidean distance, which is the type of distance most commonly used, must be distinguished, since it is simply a geometric distance in a multidimensional space. The Euclidean distance is given by the following equation:

where I _i is the individual of day i, I _j is the individual of day j, and X _ik is the observation of day i at the station k.

Performing an ascending hierarchical classification therefore consists in partitioning all the elements of the population into subsets such that each subset is well differentiated from the others. For this, an index of aggregation (distance between groups or between an individual and a group) is fixed taking into account aggregation criteria such as single linkage (minimum distance), complete linkage (maximum distance), and Ward’s method, to minimize the sum of squares of all the pairs of classes that can be formed at each step.

By cutting the tree (dendrogram) at a significant jump in the aggregation index, a good score can be obtained because the individuals grouped below the cut are close to each other, while those grouped after the cut are far apart.

3.1. Test for stationarity and independence

The software HYFRAN (^{Bobée et al., 1999}) allows to achieve a full frequential analysis for any variable for which one has observations: independent (absence of self-correlation) and identically distributed (homogeneity, stationarity and absence of outsiders). Thus, the study of the series independence and stationarity precedes the series distribution stationarity of the maximum daily rainfalls. Statistic tests of stationarity (Kendall) and independence (Wald-Wolfowitz) are applied to each series. Mann-Kendall is a non-parametric statistic test to detect trends in the time series devoid of seasonality. Both tests used are accepted for all stations at a 5% threshold (Table II). From these results, representativeness of annual maximum daily rainfalls series is accepted.

3.2. Discordancy, homogeneity, and goodness-of-fit tests

Regional estimation of hydrometeorological variables is necessary to specify the asymptotic behavior of the annual maximum daily rainfalls distribution, to reduce the sampling influence on short series, and particularly to remedy the data shortage in sites devoid of measuring stations (^{Muller, 2006}). The first stage for this procedure is the decomposition of the study region into homogeneous groups of stations. The methods commonly used in hydrology to create homogeneous groups of stations are generally based on the determination of regional indices and multivariate analysis (^{Muller, 2006}; ^{Meddi and Toumi, 2015}). The variables belonging to a homogeneous region result from the same population; thus, they follow the same probability distribution and share the same parameters. To validate the homogeneity of a region in terms of L-moments, the statistic homogeneity test proposed by ^{Hosking and Wallis (1993)} will be used. Once the regions limits are fixed, the homogeneity of the sites within the region can be validated. This step often consists in calculating the coefficients of variations, asymmetry and kurtosis for every site in the region and then comparing their variability with that of the homogeneous regional model (^{Schaefer, 1990}; ^{Cong et al., 1993}; ^{Alila, 1999}; ^{Sveinsson et al., 2000}).

Regional homogeneity was tested first by using discordance (D _i ) and heterogeneity (H _i ) tests. The results showed heterogeneity on the whole area with regard to L-Cv and L-Cs. If the variable studied (namely annual maximum daily rainfalls) does not follow the same probability distribution with the same parameters, it does not belong to the same population (Table III). Therefore, the assembly of stations in homogenous groups is necessary. To this end, we used ascending hierarchical classification, which highlighted three distinct groups (Fig. 2). Station distributions in homogeneous groups (Fig. 2) are due to the fact that such groups are naturally separated under the effect of their exposure to winds; the stations that coincide with this direction of exposure to winds receive more rainfalls. This is the case of stations in group I. On the other hand, stations in the second group are situated under shelter, which makes them receive less rainfalls despite the fact that they are located to the north, where the mountain barrier hinders the advance of dominating humid air masses, which follow a north-west direction causing abundant rainfalls. Finally, the third group is located at the meeting point of two mountain chains, in a zone characterized by flat topography where rainfall is generally stormy.

Fig. 2 Distribution of stations with homogeneous groups. 

To evaluate the degree of homogeneity of each group, series of 1000 simulations of maximum daily precipitation as per the GEV distribution were carried out. Results of the homogeneity tests for each group in terms of L-moments ratio are presented in Table III. To validate the homogeneity of the three groups in terms of L-moments, the D _i discordance and heterogeneity tests proposed by ^{Hosking and Wallis (1993)} were taken into account. The discordance test indicates if rainfall stations in the same group are significantly different from others. This statistic is calculated for the three groups. The results show that this test is concluding in the case of the three groups. The critical value of D _i = 3 is not exceeded (Table III).

Furthermore, homogeneity tests H ₁ , H ₂ and H ₃ estimate the degree of homogeneity for a given group to determine if a specific group is homogeneous (Table III). These confirmations are necessary to choose the statistic model of regional fitting. The results of homogeneity tests for each group in terms of L-moment ratios are presented in Table III, where it is shown that group I (14 stations) is homogeneous in terms of L-Cv, L-Cs and L-Ck (H < 1) (Fig. 3).

Fig. 3 L-moments diagram for stations in group 1. 

Group II is homogeneous in terms of L-Cv, L-Cs and L-Ck (H < 1). This group consists of 23 stations which are represented in Figure 4.

Fig. 4 L-moments diagram for stations in group 2. 

Group III is homogeneous in terms of L-Cs and of L-Ck (H < 1), but probably heterogeneous in terms of L-Cv. This group consists of 18 stations, as shown in Figure 5.

Fig. 5 L-moments diagram for stations in group 3. 

After testing for homogeneity in the three groups, the best model for the fitting of maximum daily precipitations was selected by using the Z^DIST statistic as suggested by ^{Hosking and Wallis (1993)}, ¹⁹⁹⁷). In addition, the L-moments diagram for each of the three groups was plotted in order to choose the most appropriate regional distribution for each group. Z^DIST statistic (Table IV) and diagrams (Figs. 3, 4 and 5) show that GEV law represents the best model for annual maximum daily precipitations fitting for the three groups. This result corroborates with the results obtained by the local approach (^{Habibi et al., 2013}; ^{Boucefiane et al., 2014}).

4.1. Design rainfall depth estimation through index-rainfall approach

The regional statistic analysis of maximum daily rainfalls allows drawing regional distribution curves for each group of stations described as being homogeneous in terms of what is shown in Figure 6. Thus, the rainfall index determined for each of the three homogeneous groups is used in Eq. (4) to calculate the maximum daily rainfall for different return periods for the stations that have short observation series, by using the p(F) values of written down in Table V. The main assumption of the rain index method was checked, verifying that the statistic distribution of maximum daily rainfalls in the three homogeneous groups is similar.

Fig. 6 Regional growth curves 

To calculate the maximum daily rainfall from one of the three groups for a given return period, the average maximum daily rainfall of the series of station i should be multiplied by the corresponding quantile function of Table V using Eq. (9).

This approach will make it possible to estimate rare frequency quantiles in stations with short or incomplete observation series. Also, these quantiles can be calculated for sites without measurement stations by using the maximum average daily rainfall that can be obtained by interpolation.

The estimation of these precipitations frequency is essential to water development and in the design of hydraulic structures. The estimation of quantiles for maximum rainfall is required for calculating flood protection works. The region is subject to repeated floods caused by heavy rainfalls. As an example, the events of October 2011 caused the death of 10 persons and huge material damage in El Bayadh city.

These numbers show the importance of estimating extreme rainfall in development studies and the regional approach as a solution to the lack of rainfall data and the inadequate quality of information measured at some stations.

4.2. Validation of use of regional growth curve for quantiles estimation

Reliability of the regional method for quantile estimation is validated by bias (BIAS) and root mean square error (RMSE) for each return period (Table VI).

In terms of bias, quantiles estimated from regional information are rather close to those locally estimated. For low return periods (T < 20 years), the bias is practically low. Beyond this threshold, the bias still remains acceptable (< 20%). The root mean square error is below 25% for the quantiles relating to the return periods below 100 yrs.; but, for higher return periods, estimations for daily maximum precipitations should be dealt with caution. These results are comparable to those by ^{Onibon et al. (2005)} in Canada, ^{Benhattab et al. (2014)} in Algeria, and ^{Malekinezhad and Zare-Garizi (2014)} in Iran. Table VII presents the deviations of local estimation to regional estimation of quantiles for the stations of Oued Taria, Mecheria and Ain Skhouna. For T ≤ 10 years, the difference is almost negligible. Beyond this threshold, regional estimate introduces a difference characterized by underestimation or overestimation of quantiles. Ain Skhouna and Oued Taria stations illustrate overestimation. The deviations on large return periods are due to the regional information effect on the estimation of L-C_V and L-C_S. When the regional estimation of the latters gives values above or below those locally estimated, the regional model tends to overestimate or underestimate the quantiles associated with the large return periods.