# Atmospheric Predictors for Annual Maximum Precipitation in North Africa

Bouchra Nasri, Yves Tramblay, Salaheddine El Adlouni, Elke Hertig, Taha Ouarda

► To cite this version:

Bouchra Nasri, Yves Tramblay, Salaheddine El Adlouni, Elke Hertig, Taha Ouarda. Atmospheric Predictors for Annual Maximum Precipitation in North Africa. Journal of Applied Meteorology and Climatology, 2016, 55 (4), pp.1063 - 1076. 10.1175/JAMC-D-14-0122.1. hal-01871879

HAL Id: hal-01871879

https://hal.umontpellier.fr/hal-01871879

Submitted on 29 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Atmospheric Predictors for Annual Maximum Precipitation in North Africa

BOUCHRA NASRI

Canada Research Chair on the Estimation of Hydrometeorological Variables, Eau Terre Environnement Research Centre,
Institut National de la Recherche Scientifique, Québec, Québec, Canada

YVES TRAMBLAY

Unité Mixte de Recherche Hydrosciences, Institut de Recherche pour le Développement, Montpellier, France

SALAHEDDINE EL ADLOUNI

Département de Mathématique et de Statistique, Université de Moncton, Moncton, New Brunswick, Canada

ELKE HERTIG

Institute of Geography, University of Augsburg, Augsburg, Germany

TAHA B. M. J. OUARDA

Canada Research Chair on the Estimation of Hydrometeorological Variables, Eau Terre Environnement Research Centre,
Institut National de la Recherche Scientifique, Québec, Québec, Canada, and Masdar Institute of Science
and Technology, Abu Dhabi, United Arab Emirates

(Manuscript received 5 May 2014, in final form 21 July 2015)

## ABSTRACT

The high precipitation variability over North Africa presents a major challenge for the population and the infrastructure in the region. The last decades have seen many flood events caused by extreme precipitation in this area. There is a strong need to identify the most relevant atmospheric predictors to model these extreme events. In the present work, the effect of 14 different predictors calculated from NCEP-NCAR reanalysis, with daily to seasonal time steps, on the maximum annual precipitation (MAP) is evaluated at six coastal stations located in North Africa (Larache, Tangier, Melilla, Algiers, Tunis, and Gabès). The generalized extreme value (GEV) B-spline model was used to detect this influence. This model considers all continuous dependence forms (linear, quadratic, etc.) between the covariates and the variable of interest, thus providing a very flexible framework to evaluate the covariate effects on the GEV model parameters. Results show that no single set of covariates is valid for all stations. Overall, a strong dependence between the NCEP-NCAR predictors and MAP is detected, particularly with predictors describing large-scale circulation (geopotential height) or moisture (humidity). This study can therefore provide insights for developing extreme precipitation downscaling models that are tailored for North African conditions.

## 1. Introduction

Heavy precipitation events are causing extensive damage to the populations and infrastructure of the countries located in the southern part of the Mediterranean basin. The last decades saw several deadly flood

events caused by extreme precipitation, including the 2001 flood near Algiers, Algeria, causing more than 700 fatalities (Argence et al. 2008), the 1969 floods in the region of Kairouan, Tunisia, with 150–400 fatalities (Poncet 1970), or the 1995 flood in the Ourika valley, Morocco, with more than 200 fatalities (Saidi et al. 2003). To better mitigate the impacts of these extreme events, there is a need to evaluate their predictability on different time scales. In particular, it is necessary to estimate their return periods in a climate change context

<footnote>Corresponding author address: Bouchra Nasri, INRS-ETE, 490 Rue de la Couronne, Québec, QC G1K 9A9, Canada.
E-mail: bouchra.nasri@ete.inrs.ca</footnote>
since several countries experienced an increased vul-
nerability to these events during the last decade (Di
Baldassarre et al. 2010). Several recent studies have
focused on seasonal precipitation and its extremes in the
Mediterranean region, with the objective of identifying
the associated large-scale patterns and the relevant
predictors (Knippertz et al. 2003; Xoplaki et al. 2004;
Martín-Vide and Lopez-Bustins 2006; Toreti et al. 2010;
Tramblay et al. 2011; Kallache et al. 2011; Hertig et al.
2012; Tramblay et al. 2012a,b; Hertig et al. 2013;
Ouachani et al. 2013; Donat et al. 2014). Indeed, to re-
solve the mismatch of scales between general circulation
models and the locations of interest for impact studies,
there is a need to develop downscaling techniques tai-
lored for extreme precipitation (Fowler et al. 2007;
Maraun et al. 2010a). To overcome the limitations of
climate models in reproducing extremes (Sillmann et al.
2013), several studies have used covariates in non-
stationary extreme precipitation frequency analysis
(e.g., Vrac and Naveau 2007; Aissaoui-Fqayeh et al.
2009; Beguería et al. 2010; Friederichs 2010; Kallache
et al. 2011; Tramblay et al. 2011; Maraun et al. 2010b;
Ouachani et al. 2013; El Adlouni and Ouarda 2009;
Cannon 2010; Ouarda and El Adlouni 2011). However,
even if several authors have shown the efficiency of at-
mospheric humidity and moisture flux as predictors for
daily rainfall modeling and downscaling (Cavazos and
Hewitson 2005; Bliefernicht and Bárdossy 2007;
Tramblay et al. 2011, 2013), the best predictors may
differ from one site to another (Kallache et al. 2011;
Hertig et al. 2013; Chandran et al. 2016). In addition, it
should be noted that the above studies have mostly ap-
plied polynomial dependence (linear or quadratic) be-
tween the covariate and the variable of interest. In the
present work, the nonstationary generalized extreme
value (GEV) model with B-spline dependent function
(Nasri et al. 2013; Chavez-Demoulin and Davison 2005)
is applied. B-spline functions are piecewise polynomial
functions that have certain advantages. A smoothing
B-spline basis is independent of the response variable
and depends only on the following information: (i) the
extent of the explanatory variable, (ii) the number and
position of the knots, and (iii) the degree of the B-spline.
These advantages make it a suitable option for use in the
GEV model with covariates to explain the effect of co-
variates on the response variable. Therefore, the goal of
this study is to identify relevant large-scale predictors
influencing the annual maximum precipitation at coastal
stations in the southern part of the Mediterranean re-
region using the GEV B-spline model. This model will
help describe the predictors' influence on the pre-
cipitation records within the study period. A number of
studies on the impact of climate variability on extreme

precipitation in the Mediterranean region have em-
ployed atmospheric–oceanic teleconnection indices
such as the North Atlantic Oscillation (NAO; Wanner
et al. 2001), the Mediterranean oscillation (Conte et al.
1989), or the western Mediterranean oscillation (WEMO)
(Knippertz et al. 2003; Vicente-Serrano et al. 2009).
However, in Morocco Tramblay et al. (2012a) observed a
possible dependence of precipitation extremes with these
indices only at 2 stations (Larache and Tangier) out of 10.
In the present study, we tested the effect of different
predictors measured with NCEP–NCAR reanalysis
data to evaluate their influence on the maximum annual
precipitation (MAP) time series. To our knowledge, no
studies have previously examined these extreme precip-
itation events and their relationship to large-scale atmo-
spheric influences in this area at the daily time step with
rain gauge data, mainly because of the limited access to
the data. Since reanalysis data are available at a spatial
scale similar to that of global circulation models (GCMs),
this study provides the first step toward the development
of extreme precipitation downscaling methods that are
tailored for North African conditions.

## 2. Datasets

In the present study, we collected long daily pre-
cipitation time series maintained by the governmental
hydrological services of Algeria [Agence Nationale des
Ressources Hydrauliques (ANRH)], Morocco [Di-
rection de la Recherche et de la Planification de l'Eau
(DRPE)], and Tunisia [Direction Générale des Re-
sources en Eau (DGRE)]. The daily data of the Melilla
station located in northern Morocco were obtained
from the European Climate Assessment and Dataset
(ECAD; http://eca.knmi.nl). The data from each station
were carefully scrutinized, in particular to look for shifts,
absurd values, and missing data (Tramblay et al. 2013).
The stations that were subsequently selected had less
than 5% missing days between September and May. The
years with more than 5% missing days during this period
were discarded. Figure 1 illustrates the geographic lo-
cation of all stations, and Table 1 presents a description
of the selected stations with long precipitation records.
Reanalysis data from the National Centers for Envi-
ronmental Prediction (NCEP; Kalnay et al. 1996; Kistler
et al. 2001) are used to compute several large-scale
predictors. Various variables were extracted to be tested
in the model; the selection of covariates is based on the
previous studies of Cavazos and Hewitson (2005),
Kallache et al. (2011), Tramblay et al. (2011), and Hertig
et al. (2013). The NCEP–NCAR reanalysis data have
been selected over more recent products such as
ERA-Interim because the time span of the NCEP–
<figure><img src="image_2.png" /><figcaption>FIG. 1. Geographic location of all stations (three selected stations in Morocco, one in Algeria, and two in Tunisia).</figcaption></figure>

NCAR reanalysis is larger and encompasses the whole period of observations, up to the present. The advantages of recent reanalysis products are manifold, including new atmospheric and assimilation systems and finer grid spacing. However, they cover only the recent period (from 1979 to present for MERRA, CFS reanalysis, or ERA-Interim; Hofer et al. 2012). Thus the advantage of the choice of the NCEP-NCAR reanalysis is to cover the whole period where observations are available using a single reanalysis product in order to identify possible associations with large-scale climate dynamics. It must be noted that our goal is not to check the adequacy of the particular NCEP-NCAR reanalysis product but to evaluate if the distribution of extreme precipitation can be related to large-scale predictors, if the same predictors are valid for different stations, and at which time scales the large-scale forcings are relevant. For these objectives, the bias of a given reanalysis product compared to other products is of little relevance, and because of the strong interannual variability of precipitation in North Africa, it is important to evaluate the relationships with large-scale dynamics over long time periods to obtain robust results. The gridded NCEP-NCAR data (six grid cells covering the

study area) were interpolated by the inverse distance method to the station locations in order to provide individual descriptor sets for each station. The selected variables include the following:

* geopotential height at 500 and 850 hPa (geopot_500 and geopot_850),

* vertical velocity at 500 and 850 hPa and at surface (omega_500, omega_850, and omega_surf),

* potential temperature at surface (ptemp_surf),

* precipitable water content (pwater_surf),

* relative humidity at surface (rhum_surf),

* specific humidity at 500 and 850 hPa (shum_500 and shum_850),

* mean sea level pressure (slp_surf),

* surface temperature (temp_surf),

* zonal wind at surface (uwind_surf), and

* meridional wind at surface (vwind_surf).

The homogeneity of these covariates has been assessed by following Pettitt (1979) and by using the modified version of the standard normal homogeneity test (SNHT) by Khaliq and Ouarda (2007). Indeed, the gradual introduction of satellite data into reanalysis products can introduce an artificial changepoint leading

<table><caption>TABLE 1. Description of the selected stations with long records of precipitation.</caption><thead><tr><th>Station</th><th>Country</th><th>Lat</th><th>Lon</th><th>Alt (m)</th><th>Record length (yr)</th><th>Starting year</th><th>Ending year</th></tr></thead><tbody><tr><td>Algiers</td><td>Algeria</td><td>36.74°</td><td>3.06°</td><td>140</td><td>47</td><td>1951</td><td>2005</td></tr><tr><td>Larache</td><td>Morocco</td><td>35.18°</td><td>-6.15°</td><td>5</td><td>51</td><td>1942</td><td>2011</td></tr><tr><td>Tangier</td><td>Morocco</td><td>35.77°</td><td>-5.80°</td><td>5</td><td>33</td><td>1972</td><td>2006</td></tr><tr><td>Melilla</td><td>Morocco</td><td>35.29°</td><td>-2.94°</td><td>47</td><td>46</td><td>1907</td><td>2009</td></tr><tr><td>Gabès</td><td>Tunisia</td><td>33.88°</td><td>10.10°</td><td>4</td><td>57</td><td>1950</td><td>2009</td></tr><tr><td>Tunis</td><td>Tunisia</td><td>36.83°</td><td>10.23°</td><td>66</td><td>58</td><td>1950</td><td>2009</td></tr></tbody></table>
to the false detection of trends or homogeneity breaks (Sterl 2004). The Pettitt and SNHT tests agree only on a significant changepoint at the 5% level in relative humidity, in 1957 for SNHT and in 1963 for the Pettitt test. Therefore, no changepoints are detected in the beginning of the 1980s following the introduction of satellite data.

These covariates are considered at different time steps. In the first case, we considered the maximum observed daily precipitation during the extended winter season (October–March) and the simultaneous daily covariate in the reanalysis data associated with this extreme rainfall event. This gives one observation of maximum daily winter rainfall and its associated covariate for each year (hereafter case 1). In case 2, we considered the maximum winter precipitation and the average value of each covariate 5 days before the date of the annual maximum rainfall during the winter. These two cases can be considered as dealing with the short-term effect of covariates on MAP. In case 3, we calculated the 30-day average of the covariate before the date of maximum winter precipitation. Finally, in case 4, we considered the maximum daily winter precipitation and the value of each covariate for the entire season (October–March average). These last two cases can be considered as dealing with the long-term effect of covariates on MAP.

### 3. Methods

For modeling extreme rainfall events, we used the GEV distribution (Coles 2001). The role of the GEV distribution is to describe a sample that follows a maximum of distributions introduced by Fisher and Tippett (1928). The GEV distribution is flexible and has been the subject of several theoretical studies and applications for modeling extreme flood, precipitation, and wind events (El Adlouni et al. 2007; Hundecha et al. 2008). The development of stationary GEV distribution models for univariate extreme value analysis can be found in the literature (Coles 2001; Olsen et al. 1999). The use of this distribution in the frequency analysis of extreme events is based on a number of specific hypotheses concerning the variable of interest. Indeed, the observations must be independent and identically distributed. However, the stationarity assumption is often not met for observed hydroclimatic datasets (Khaliq et al. 2006). In this case, the distribution parameters and the distribution itself could be changing in time. Therefore, it is essential to develop the GEV model in the multivariate space, where extreme events can be associated with other variables. To model the relationship between the covariates and the extreme variable of interest, we can

use the GEV B-spline approach (Nasri et al. 2013). This approach has been developed to describe the association of an external covariate with the variable of interest. The estimation of the parameters of the GEV B-Spline model is done in a Bayesian framework to obtain the posterior distribution by applying Markov chain Monte Carlo (MCMC) algorithms.

#### a. The GEV distribution

The GEV distribution is characterized by three parameters: the location $\mu$, scale $\sigma$, and shape $\xi$ parameters. Depending on the value of the shape parameter we have three types of extreme value distributions—namely, the Gumbel ($\xi = 0$), Fréchet ($\xi > 0$), and Weibull ($\xi < 0$). Considering a sample $Y = (y_1, \dots, y_n)$, the GEV distribution function is as follows:

$$ 
\begin{aligned}
F(y, \mu, \sigma, \xi) &= \exp\left\{-\left[1 + \xi\left(\frac{y-\mu}{\sigma}\right)\right]^{-1/\xi}\right\} \quad \xi \neq 0 \\
F(y, \mu, \sigma, \xi) &= \exp\left[-\exp\left(-\frac{y-\mu}{\sigma}\right)\right] \quad \xi = 0.
\end{aligned}
\quad (1) 
$$

This classical GEV distribution is based on the stationarity assumption and does not consider the dependence of extreme events on other variables. In the following section, the nonstationary GEV approach is presented to consider the effect of a covariate on extreme values.

#### b. The nonstationary GEV B-spline model

In the nonstationary case of a GEV distribution, the parameters of the GEV distribution are assumed to change in time or depend on a covariate. In the present form of the GEV, parameters $\sigma$ and $\xi$ are assumed to be constant. Having a random variable $Y$ that follows the GEV($\mu_x$, $\sigma$, $\xi$) and a vector of $p$ covariates given by $\mathbf{X} = (X_1, X_2, \dots, X_p)$, the location parameter of the GEV is written as follows:

$$ \mu_x = \sum_{i=1}^{p} f_i(X_i) = f_1(X_1) + f_2(X_2) + \dots + f_p(X_p), \quad (2) $$

where $f_i$ is a function that represents the relationship between the parameter and the covariates $X_i$. This function can be described by the following B-spline function:

$$ f_i(x_i) = \beta_{0,i} + \sum_{j=1}^{m} \beta_{j,i} B_{j,i,d}(x_i), \quad (3) $$

where $B_{j,i,d}(x)$ is a polynomial function of degree $d$ and $m$ is the number of control points (Nasri et al. 2013). Therefore, Eq. (2) can be rewritten as follows:
$$ \mu_x = \sum_{i=1}^{p} f_i(x_i) = \sum_{i=1}^{p} \left[ \beta_{0,i} + \sum_{j=1}^{m} \beta_{i,j} B_{j,i,d}(x_i) \right]. \quad (4) $$

The predictors' interaction can be expressed in our model by using multivariate B-spline functions (de Boor 2001). These functions allow considering the correlation between the predictors. In this study 14 predictors are used. Consequently, in order to simplify the model we did not consider the interaction between predictors.

### c. Parameter estimation

In this study, the estimation of the parameters of the GEV B-spline model is carried out in a Bayesian framework. In the Bayesian approach, the unknown parameters are not constant and are considered as random variables with a prior distribution $\pi(\theta)$.

Bayes's theorem therefore gives the following definition of the posterior distribution of these parameters:

$$ f(\theta | y) = \frac{f(y | \theta) \times \pi(\theta)}{f(y)}, \quad (5) $$

where

$$ \theta = (\mu_x, \sigma, \xi) = (\beta, \sigma, \xi) \quad \text{and} \quad \beta = \begin{pmatrix} \beta_0 \\ \beta \end{pmatrix}. \quad (6) $$

According to Nasri et al. (2013), we choose a multivariate normal distribution as a prior for a location parameter $\beta \sim N(0, \Sigma_\beta \times I)$, a noninformative prior for scale parameter $1/\sigma$, and a beta distribution as prior for a shape parameter $[\beta(6, 9)]$.

The posterior distribution is written as follows:

$$ f(\theta | y) \propto \frac{1}{\sigma} \left\{ 1 + \xi \left[ \frac{y - \sum_j (1 B) \begin{pmatrix} \beta_0 \\ \beta \end{pmatrix}}{\sigma} \right] \right\}^{-(1/\xi)-1} \exp \left\{ - \left\{ 1 + \xi \left[ \frac{y - \sum_i (1 B) \begin{pmatrix} \beta_0 \\ \beta \end{pmatrix}}{\sigma} \right] \right\}^{-(1/\xi)} \right\} \\ \left[ 2\pi \det(\Sigma_\beta) \right]^{-(k/2)} \exp \left( - \frac{\|\beta\|^2}{2\sigma^2} \right) \frac{1}{\sigma_\xi \sqrt{2\pi}} \exp \left[ - \frac{(\xi - 0.1)^2}{2\sigma_\xi^2} \right] \frac{1}{\sigma}. \quad (7) $$

The GEV B-spline model, which takes into account nonstationarity and nonlinearity, offers a great flexibility and takes into account the heavy-tailed character of the extreme distribution.

The posterior distribution is estimated by the Metropolis-Hasting algorithm (see appendix).

To select the number of knots (kt; 1 kt = 0.51 m s⁻¹) and the degree of B-spline functions used in this study, we compared several combinations of degrees and knots using the maximum likelihood method. The following algorithm explains how these parameters are chosen:

(i) set $d \in [1, 10]$ and $m \in [1, 10]$,

(ii) calculate Eq. (7) for all combinations of $(d, m)$, and

(iii) choose values of $(d, m)$ that maximize Eq. (7).

In this case, we apply the B-spline functions with 3° and 3 kt. This choice was found to be optimal for the majority of the stations' data.

### d. Validity of the model with covariates

To validate the influence of covariates on the variable of interest, the log likelihood of the GEV B-spline (M1) model and the stationary GEV (M0) model

(without covariates) are compared using the test of deviance:

$$ D = 2[l(M1) - l(M0)], \quad (8) $$

where $l$ is the maximum log likelihood function for model M. Large values of $D$ indicate that model M1 is more adequate at representing the data than model M0. The $D$ statistic is distributed according to a chi-square distribution $\chi^2$, with $\nu$ degrees of freedom, where $\nu$ is the difference between the number of parameters of the M1 and M0 models. For a given $\alpha$ confidence level, we reject H0 hypothesis (H0: M1 and M0 are similar) when $D \ge \chi_{1-\alpha}^2$. This statistic is often used to compare two models when one model is a special case of the other (M0 $\subset$ M1; Coles 2001; El Adlouni and Ouarda 2009). This test accounts for differences in model complexity to avoid overfitting.

### e. Quantile estimation

The MCMC algorithm also produces the conditional quantile distribution for an observed value $x_0$ of the covariate $X_t$. Indeed, for each iteration $i$ of the MCMC algorithm $i = 1, \dots, N$, the quantiles corresponding to
<figure><img src="image_3.png" /><figcaption>FIG. 2. Variation of all MAP series vs time for selected stations.</figcaption></figure>

the nonexceedance probability $1-p$, $x_{p,x_0}^{(i)}$, and the parameter vector $[\mu_{x_0}^{(i)}, \sigma_{x_0}^{(i)}, \xi_{x_0}^{(i)}]$ are estimated using the inverse of the cumulative distribution function of the GEV distribution:

$$y_{p,x_0}^{(i)} = \mu_{x_0}^{(i)} - \frac{\sigma}{\xi^{(i)}} \{1 - [\log(1-p)]^{-\xi^{(i)}}\}, \quad (9)$$

where $\mu_{x_0}^{(i)}$ is the position parameter conditional on the particular value $x_0$ of X.

# 4. Results

## a. Tests for independent and identically distributed random variables

In the first step of using the nonstationary GEV model (in this case GEV B-spline) we checked stationarity, homogeneity, and independence using the Mann-Kendall (Mann 1945), Mann-Whitney (Wilcoxon 1945), and Wald-Wolfowitz tests (Wald and Wolfowitz 1940), respectively, for MAP series for each station. The results of these tests showed that all MAP series are non-stationary, except for the Algiers station where a negative trend is observed and is significant at the 5% level, as previously observed by Reiser and Kutiel (2011). However, all the time series of MAP respect the hypotheses of homogeneity and randomness. Figure 2 shows the variation of all MAP series versus time, and Fig. 3 shows the monthly frequency of occurrence of annual maximum daily precipitation.

## b. Predictors from reanalysis data

We selected the 14 NCEP-NCAR covariates extracted from reanalysis (see section 2) and tested our models with these covariates considering the four time scales (cases 1-4). The negative log likelihood and deviance between the GEV B-spline model and the stationary GEV model are analyzed to detect the influence of NCEP-NCAR predictors on extreme rainfall for each of the four cases. To avoid overfitting, each covariate is considered separately in the GEV B-spline model. This allows evaluating whether each covariate leads to a better fit than the stationary GEV model. As there are 14 covariates for each case, the results are presented in Table 2 only for the significant covariates on MAP at each station at the 5% and 10% significance levels, according to the test of deviance. We note that all 14 covariates, depending on the station, are selected into nonstationary GEV models that better reproduce extreme precipitation than a standard stationary model, with both short-term and long-term associations. Overall, a similar number of significant covariates is selected for the four cases tested (i.e., from daily to seasonal averages of covariates), with 11 covariates identified for case 1, 13 for case 2, 16 for case 3, and 13 for case 4. This shows that all considered covariates may have an impact on extreme daily precipitation at different time steps, from daily values to seasonal averages.

It is observed that the geopotential height (geopot_500 and geopot_850) generally affects rainfall at all
<figure><img src="image_4.png" /><figcaption>FIG. 3. Monthly frequencies of occurrence for daily MAP in each selected station.</figcaption></figure>

stations excluding Melilla (station in northern Morocco). For the stations of Tangier and Larache, the geopotential heights have a short-term association with MAP (case 1 and case 2). On the other hand, for the Algiers station, these variables have an influence

generally at the seasonal time scale (case 4), and for the stations of Tunisia (Gabès and Tunis) these variables influence MAP in both the short and long term (cases 1, 2, and 4). Humidity predictors (rhumSurf, shum_500, and shum_850) generally influence precipitation at all

<table><caption>TABLE 2. The significant covariates at 5% and 10% significance levels for each station.</caption><thead><tr><th>Station</th><th>Predictors significant at 5% (case)</th><th>Predictors significant at 10% (case)</th><th>Station</th><th>Predictors significant at 5% (case)</th><th>Predictors significant at 10% (case)</th></tr></thead><tbody><tr><td rowspan="2">Algiers</td><td>SlpSurf (case 4)</td><td rowspan="2">pwaterSurf (case 1)</td><td rowspan="2">Melilla</td><td>omega_Surf (case 1)</td><td rowspan="2">rhumSurf (case 4)</td></tr><tr><td>geopot_500 (case 4)</td><td>omega_500 (case 2)</td></tr><tr><td rowspan="7">Gabès</td><td>geopot_500 (case 1)</td><td rowspan="2">geopot_850 (case 2)</td><td rowspan="7">Tangier</td><td>omega_850 (case 2)</td><td>geopot_850 (case 1)</td></tr><tr><td>geopot_850 (cases 1, 2, and 4)</td><td>rhumSurf (case 2)</td><td>temp_Surf (case 4)</td></tr><tr><td>pwaterSurf (case 3)</td><td>slpSurf (cases 1 and 4)</td><td>omega_850 (case 2)</td><td>uwindSurf (case 4)</td></tr><tr><td>slpSurf (cases 1 and 4)</td><td>shum_500 (case 3)</td><td>omega_850 (case 2)</td><td>temp_Surf (case 4)</td></tr><tr><td>shum_500 (case 3)</td><td>shum_850 (case 1)</td><td>uwindSurf (case 2)</td><td>temp_Surf (case 4)</td></tr><tr><td>shum_850 (case 1)</td><td>vwindSurf (case 3)</td><td>omega_850 (case 2)</td><td>geopot_500 (case 4)</td></tr><tr><td>geopot_500 (case 1)</td><td>pwaterSurf (case 2)</td><td>pwaterSurf (cases 1 and 2)</td><td>rhumSurf (case 1)</td></tr><tr><td rowspan="9">Larache</td><td>temp_Surf (case 2)</td><td rowspan="9">pwaterSurf (case 4)</td><td rowspan="9">Tunis</td><td>temp_Surf (case 4)</td><td>shum_850 (case 1)</td></tr><tr><td>omega_500 (case 2)</td><td>shum_500 (case 2)</td><td>geopot_850 (case 2)</td></tr><tr><td>shum_850 (case 2)</td><td>shum_850 (case 2)</td><td>omega_500 (case 3)</td></tr><tr><td>geopot_500 (case 3)</td><td>slpSurf (case 3)</td><td>omega_850 (case 4)</td></tr><tr><td>temp_Surf (case 3)</td><td>omega_850 (case 4)</td><td>omega_500 (case 3)</td></tr><tr><td>pwaterSurf (case 3)</td><td>omega_850 (case 4)</td><td>temp_Surf (case 4)</td></tr><tr><td>rhumSurf (case 3)</td><td>omega_850 (case 4)</td><td>temp_Surf (case 4)</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>FIG. 4. Contributions of the 14 NCEP-NCAR reanalysis covariates on the two principal components (F1 and F2) in selected station (results for case 1). The numbers in the parentheses represent the percentage of explained variance for the represented axes (F1 and F2).</figcaption></figure>

stations, excluding Algiers. For stations in Morocco, these predictors appear in almost all cases (cases 1 and 3 for the Tangier station, case 4 for the Melilla station, and cases 1, 2, and 3 for the Larache station). For stations in Tunisia, these predictors have both short- and long-term effects on MAP time series at the stations of Gabès (cases 1, 3, and 4) and Tunis (cases 1, 2, and 3). Velocity predictors (omega_500, omega_850, and omega Surf) have more effects on precipitation in Morocco. We see a stronger influence of these predictors on rainfall in Morocco. Wind predictors (uwind Surf and vwind Surf) influence the MAP only at the Tangier station. Overall, we note the small influence of wind covariates on precipitation extremes at all stations. The potential temperature at the surface (ptemp Surf) influences MAP at stations located in Morocco for the long term (cases 3 and 4). The surface temperature influences the MAP in Morocco at different time scales (case 3 for Melilla, case 1 for Larache, and case 4 for Tangier stations). The precipitable water content has an influence on MAP at all stations, usually only for the short-term cases (1 and 2) at all stations. The mean sea level pressure influences MAP in only the Tunis, Gabès, and Algiers stations, generally in the long-term cases 3 and 4.

### c. Principal analysis of components for NCEP-NCAR predictors

After the analysis of the dependence of MAP with individual covariates, the possible relationships are also investigated in a multivariate context. Principal component analysis (PCA; Preisendorfer 1988) is used to that end. The reason for using PCA is to take into consideration the common signals in multivariate datasets. PCA represents a method for dimensionality reduction. PCA has been used for this purpose in many other studies (e.g., Wetterhall et al. 2005; Maraun et al. 2010a). The objective of this analysis is to summarize as much information as possible by transforming interrelated variables into new components (principal components) that are uncorrelated with each other. In this study, we first applied PCA on the 14 covariates for each station. Figures 4 and 5 show the results of the projections of the 14 covariates on the first and second components (F1 and F2, respectively) for the first and last case in each station. A number of criteria, such as the Kaiser criterion (Kaiser 1960), can be used for the selection of the factorial axis. The Kaiser criterion lies on the factorial axis choices, where their eigenvalues are greater than 1. In the present study, we noticed that for all factors that
<figure><img src="image_6.png" /><figcaption>FIG. 5. As in Fig. 4, but for case 4.</figcaption></figure>

have an eigenvalue greater than 1, those are generally factors 1 and 2. This justifies the choice of two factorial axes. It can be seen that at all stations, there are significant correlations between the covariates, depending on the case (1–4) considered for temporal aggregation. To avoid overfitting, each component is considered separately in the GEV B-spline model to evaluate if it provides a better fit than a stationary GEV model. The results show, first, that most variables contribute to the formation of the components F1 and F2, with some covariates having a larger contribution, such as the geopotential height (geopot_500 and geopot_850), velocity (omega_500, omega_850, and omega Köppen-Geiger), and humidity (rhum Köppen-Geiger, shum_500, and shum_850). Predictors such as uwind Köppen-Geiger and vwwind Köppen-Geiger contribute more in stations close to the Mediterranean coast such as Tangier, Tunis, and Gabès. We then applied the GEV B-spline model of the MAP series for each station and each case using F1 and F2 as covariates. Next, we calculated the deviance between the results of the GEV B-spline and GEV0 models to investigate the influence of these components on MAP data. Table 3 shows the results of deviance with a threshold of 5% and 10%. According to these results one can see that, at all stations, there is at least one component (F1 or F2) that influences the MAP

series. For stations in Morocco (Larache, Melilla, and Tangier), we note that the component that contains more information about the geopotential height, humidity, velocity, and wind has the largest influence on MAP. For the station in Algeria (Algiers), the MAP is more influenced by the geopotential height predictors rather than by others. For stations in Tunisia (Tunis and Gabès), we can see that the influence of geopotential

<table><caption>TABLE 3. The results of the deviance for PCA.</caption><thead><tr><th>Station</th><th>Significant component</th><th>Case</th></tr></thead><tbody><tr><td>Algiers</td><td>F1 at 5%</td><td>4</td></tr><tr><td>Larache</td><td>F1 at 5%</td><td>1</td></tr><tr><td rowspan="2">Tangier</td><td>F1 at 5%</td><td>1</td></tr><tr><td>F2 at 5%</td><td>1</td></tr><tr><td rowspan="3">Melilla</td><td>F2 at 5%</td><td>2</td></tr><tr><td>F1 at 10%</td><td>3</td></tr><tr><td>F2 at 5%</td><td>4</td></tr><tr><td rowspan="2">Tunis</td><td>F1 at 10%</td><td>3</td></tr><tr><td>F2 at 5%</td><td>3</td></tr><tr><td rowspan="4">Gabès</td><td>F1 at 10%</td><td>1</td></tr><tr><td>F2 at 5%</td><td>1</td></tr><tr><td>F1 at 5%</td><td>2</td></tr><tr><td>F2 at 5%</td><td>2</td></tr><tr><td>F2 at 5%</td><td>4</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>FIG. 6. Example of nonstationary and stationary median for each station using the first or the second principal component analysis as covariates.</figcaption></figure>

height, velocity, temperature, and sea surface pressure on MAP is important.

### d. Quantile estimation

We can also see the impact of the covariates on the estimated quantile level for each of the models. In the case of the GEV B-spline model, quantile values depend not only on the nonexceedance probability $1 - p$ but also on the covariate values. This allows computing quantiles on a seasonal or annual basis, depending on the values of the covariates. To demonstrate the covariates' impact on quantile values, we show some quantile estimation examples for each station. Figure 6 displays a nonstationary quantile estimation example for each station for the 2-yr return period (nonexceedance probability = 0.5), which represents the median value of MAP.

For each station, we observed different values of the 2-yr quantiles estimated with the GEV B-spline model since quantile values are dependent on covariates. In contrast, the GEV0 model provides just one estimate for the 2-yr quantile (e.g., for the Algiers station, the median precipitation value corresponding to the 2-yr quantile is 100 mm for the GEV B-spline and 64 mm for the GEV0 model, and for the Tunis station, the stationary quantile is equal to 50 mm and the median of the nonstationary quantiles is equal to 70 mm). According to this figure, we notice that the covariate-dependent quantile values are more flexible and allow

reaching more extreme data values, unlike the stationary quantiles that do not take into consideration the interannual climatic variability. The estimated quantiles show the advantage of incorporating additional information into nonstationary models.

## 5. Conclusions

In this work, the influences of climatic variables such as geopotential height, pressure, or temperature on maximum annual daily precipitation have been studied at six stations located in North Africa with long precipitation time series. A total of 14 variables were computed from NCEP-NCAR reanalysis data. To study the influence of these covariates at the different stations, the GEV B-spline model (Nasri et al. 2013) was used. The originality of this model, as opposed to other nonstationary models, is that it takes into consideration the nonstationary and the nonlinear temporal fluctuations of covariates. Nonstationary models, such as the GEV1 (linear dependence) and the GEV2 (quadratic dependence), define in advance the form of dependence between the variable of interest and the covariates. On the other hand, the GEV B-spline model takes into consideration all continuous dependence forms between the covariates and the variable of interest.

The results of this study are divided into two parts. In the first part, the possible dependencies between the
maximum annual precipitation and each of the individual climatic covariates were considered. The GEV B-spline model was used to detect these dependencies, and the deviance likelihood ratio test was used to identify the nonstationary models with covariates that provide an improvement in comparison to stationary models in each station. In the second part, the combined dependencies were analyzed using principal component analysis of the different atmospheric predictors. From the results of the principal component analysis, we analyzed the influence of the combined variables using the two principal components (F1 and F2) for each station in the GEV B-spline model. Our results indicate that no single combination of atmospheric predictors is optimal for stations. The relevant covariates may vary from one station to another and also depend on the considered time scale, from daily to annual averages. These results are consistent with the fact that extreme precipitation is a process exhibiting a high spatiotemporal variability between different locations. Given this variability, it must be noted that the covariates describing the moisture flux in the atmosphere (relative or specific humidity) or in atmospheric circulation (pressure and geopotential heights) are often selected in the different stations as valid predictors. During winter, when most of the annual maximum precipitation occurs, geopotential height might be more important because of the southerly position of the extratropical westerlies. In other seasons thermodynamic predictors like humidity may gain significance because of the convective nature of precipitation in these seasons.

The present work provides a first step prior to the development of statistical downscaling methods tailored for extreme precipitation in North Africa. The next step would be to use GCM outputs to first validate the method in the present climate, with the covariates that are correctly reproduced in historical climate simulations, and then to make future projections. However, in this case the use of the nonstationary GEV model with B-spline functions would probably be less appropriate because of some limitations: (i) the introduction of several covariates within these types of models increases the number of hyperparameters, which increases the number of parameters to estimate as well as the estimation errors; (ii) the interactions between the predictors make the model much more complex since we need to take into consideration multivariate spline functions (de Boor 2001) or use some decisional model, such as an artificial neural network as in Cannon (2010); and (iii) this type of model allows the description of the impact of covariates on the variable of interest and is not able to use them for prediction outside this period of study.

Consequently, an alternative to this type of model is quantile regression methods (Buchinsky 1998). Unlike linear regression, which results in the estimation of the conditional mean for the response variable given certain values of predictor variables, quantile regression aims at estimating either the conditional median or other quantiles of the response variable. Quantile regression was considered by Jagger and James (2006) for wind speed and by Friederichs and Andreas (2008) for precipitation, based on several climatic covariates. Future work can focus on the comparison of extreme value models and the quantile regression approach to distinguish the relative benefits of the use of these two types of models for downscaling purposes.

*Acknowledgments.* The datasets were provided by the Agence Nationale des Ressources Hydrauliques (Algeria), Direction de la Recherche et de la Planification de l'Eau (Morocco), Direction Générale des Ressources en Eau (Tunisia), and European Climate Assessment and Dataset. Special thanks are given to H. Ben-Mansour, R. Bouaicha, L. Behlouli, K. Benhattab, R. Taibi, and K. Yaalaoui for their helpful contribution to database collection. The authors are indebted to editor Thomas Mote and to two anonymous reviewers whose comments helped considerably improve the quality of the manuscript.

# APPENDIX

## MCMC Algorithm for GEV B-Spline Model

The basic idea of the MCMC method is, for each parameter, to construct a Markov chain with the posterior distribution being a stationary and ergodic distribution. After running the Markov chain of size *N* for a given burn-in period *N*₀, one obtains a sample from the posterior distribution *f*(*θ*|*y*). One popular method for constructing a Markov chain is via the Metropolis-Hastings (MH) algorithm (Metropolis et al. 1953; Hastings 1970). We simulated the realizations from the posterior distribution by way of a single-component MH algorithm (Gilks et al. 1996). Each parameter was updated using a random-walk MH algorithm with a Gaussian proposal density centered at the current state of the chain. Some techniques to assess the convergence of the MCMC methods, such as the Raftery and Lewis diagnostic (Raftery and Lewis 1992, 1995) and subsampling methods (El Adlouni et al. 2006), make it possible to determine the length of the chain and the burn-in time. In all cases, the convergence methods indicated that the Markov chains converged within a few iterations. In this study, we
considered chains of size $N = 15000$ and a burn-in period of $N_0 = 8000$ runs. In every case, a sample of $N - N_0 = 7000$ values is collected from the posterior of each of the elements of $\theta$.

The principal steps of the MH algorithm can be summarized as follows:

(i) Initialization: assign initial value $\theta^{(0)}$ and choose an arbitrary proposal probability density $Q(\theta^* | \theta)$. In this case we propose a multivariate normal distribution.

(ii) For each iteration $t$, generate $\theta^*$, a candidate for the next sample, by picking from the distribution $Q(\theta^* | \theta_t)$.

(iii) Calculate the acceptance ratio, given by $\alpha(\theta^*, \theta_t) = [\pi(\theta^* | y)/\pi(\theta_t | y)]$.

(iv) If $\alpha \ge 1$, then the candidate is more likely than $\theta_t$; automatically accept the candidate by setting $\theta_{t+1} = \theta^*$. Otherwise, accept the candidate with probability $\alpha$; if the candidate is rejected, set $\theta_{t+1} = \theta_t$ instead.

REFERENCES

Aissaoui-Fqaych, I., S. El Adlouni, T. B. M. J. Ouarda, and A. St-Hilaire, 2009: Non-stationary lognormal model development and comparison with the non-stationary GEV model. *Hydrol. Sci. J.*, **54**, 1141–1156, doi:10.1623/hysj.54.6.1141.

Argence, S., D. Lambert, E. Richard, J. P. Chaboureau, and N. Söhne, 2008: Impact of initial condition uncertainties on the predictability of heavy rainfall in the Mediterranean: A case study. *Quart. J. Roy. Meteor. Soc.*, **134**, 1775–1788, doi:10.1002/qj.314.

Beguería, S., M. Angulo-Martínez, M. Vicente-Serrano, J. I. López-Moreno, and A. El-Kenawy, 2010: Assessing trends in extreme precipitation events intensity and magnitude using non-stationary peaks-over-threshold analysis: A case study in northeast Spain from 1930 to 2006. *Int. J. Climatol.*, **31**, 2102–2114, doi:10.1002/joc.2218.

Bliefernicht, J., and A. Bárdossy, 2007: Probabilistic forecast of daily areal precipitation focusing on extreme events. *Nat. Hazards Earth Syst. Sci.*, **7**, 263–269, doi:10.5194/nhess-7-263-2007.

Buchinsky, M., 1998: The dynamics of changes in the female wage distribution in the USA: A quantile regression approach. *J. Appl. Econ.*, **13**, 1–30, doi:10.1002/(SICI)1099-1255(199801/02)13:1<1::AID-JAE474>3.0.CO;2-A.

Cannon, A. J., 2010: A flexible nonlinear modelling framework for nonstationary generalized extreme value analysis in hydroclimatology. *Hydrol. Process.*, **24**, 673–685, doi:10.1002/hyp.7506.

Cavazos, T., and B. C. Hewitson, 2005: Performance of NCEP-NCAR reanalysis variables in statistical downscaling of daily precipitation. *Climate Res.*, **28**, 95–107, doi:10.3354/cr028095.

Chandran, A., G. Basha, and T. B. M. J. Ouarda, 2016: Influence of climate oscillations on temperature and precipitation over the United Arab Emirates. *Int. J. Climatol.*, **36**, 225–235, doi:10.1002/joc.4339.

Chavez-Demoulin, V., and A. Davison, 2005: Generalized additive modeling of sample extremes. *J. Roy. Stat. Soc.*, **54C**, 207–222, doi:10.1111/j.1467-9876.2005.00479.x.

Coles, G. S., 2001: *An Introduction to Statistical Modeling of Extreme Values*. Springer, 209 pp.

Conte, M., A. Giuffrida, and S. Tedesco, 1989: The Mediterranean oscillation: Impact on precipitation and hydrology in Italy. *Proc. Conf. on Climate, Water*, Helsinki, Finland, Academy of Finland, 121–137.

de Boor, C., 2001: *A Practical Guide to Splines*. Springer, 348 pp.

Di Baldassarre, G., A. Montanari, H. Lins, D. Koutsoyiannis, L. Brandimarte, and G. Blöschl, 2010: Flood fatalities in Africa: From diagnosis to mitigation. *Geophys. Res. Lett.*, **37**, L22402, doi:10.1029/2010GL045467.

Donat, M. G., and Coauthors, 2014: Changes in extreme temperature and precipitation in the Arab region: Long-term trends and variability related to ENSO and NAO. *Int. J. Climatol.*, **34**, 581–592, doi:10.1002/joc.3707.

El Adlouni, S., and T. B. M. J. Ouarda, 2009: Joint Bayesian model selection and parameter estimation of the GEV model with covariates using birth–death Markov chain Monte Carlo. *Water Resour. Res.*, **45**, W06403, doi:10.1029/2005WR004545.

—, A. C. Favre, and B. Bobée, 2006: Comparison of methodologies to assess the convergence of Markov chain Monte Carlo methods. *Comput. Stat. Data Anal.*, **50**, 2685–2701, doi:10.1016/j.csda.2005.04.018.

—, T. B. M. J. Ouarda, X. Zhang, R. Roy, and B. Bobée, 2007: Generalized maximum likelihood estimators for the non-stationary generalized extreme value model. *Water Resour. Res.*, **43**, W03410, doi:10.1029/2005WR004545.

Fisher, R. A., and L. H. C. Tippett, 1928: Limiting forms of the frequency distribution of the largest or smallest member of a sample. *Math. Proc. Cambridge Philos. Soc.*, **24**, 180–190, doi:10.1017/S0305004100015681.

Fowler, H. J., S. Blenkinsop, and C. Tebaldi, 2007: Linking climate change modelling to impact studies: Recent advances in downscaling techniques for hydrological modelling. *Int. J. Climatol.*, **27**, 1547–1578, doi:10.1002/joc.1556.

Friederichs, P., 2010: Statistical downscaling of extreme precipitation events using extreme value theory. *Extremes*, **13**, 109–132, doi:10.1007/s10687-010-0107-5.

—, and A. Hense, 2008: A probabilistic forecast approach for daily precipitation totals. *Wea. Forecasting*, **23**, 659–673, doi:10.1175/2007WAF2007051.1.

Gilks, W. R., S. Richardson, and D. J. Spiegelhalter, 1996: *Markov Chain Monte Carlo in Practice*. Chapman and Hall, 512 pp.

Hastings, W. K., 1970: Monte Carlo sampling methods using Markov chains and their applications. *Biometrika*, **57**, 97–109, doi:10.1093/biomet/57.1.97.

Hertig, E., A. Paxian, G. Vogt, S. Seubert, H. Paeth, and J. Jacobeit, 2012: Statistical and dynamical downscaling assessments of precipitation extremes in the Mediterranean area. *Meteor. Z.*, **21**, 61–77, doi:10.1127/0941-2948/2012/0271.

—, S. Seubert, A. Paxian, G. Vogt, H. Paeth, and J. Jacobeit, 2013: Changes of total versus extreme precipitation and dry periods until the end of the twenty-first century: Statistical assessments for the Mediterranean area. *Theor. Appl. Climatol.*, **111**, 1–20, doi:10.1007/s00704-012-0639-5.

Hofer, M., B. Marzeion, and T. Mölg, 2012: Comparing the skill of different reanalyses and their ensembles as predictors for daily air temperature on a glaciated mountain (Peru). *Climate Dyn.*, **39**, 1969–1980, doi:10.1007/s00382-012-1501-2.
Hundecha, Y., A. St-Hilaire, T. B. M. J. Ouarda, S. El Adlouni, and P. Gachon: 2008: A nonstationary extreme value analysis for the assessment of changes in extreme annual wind speed over the Gulf of St. Lawrence, Canada. *J. Appl. Meteor. Climatol.*, **47**, 2745–2759, doi:10.1175/2008JAMC1665.1.

Jagger, T. H., and B. E. James, 2006: Climatology models for extreme hurricane winds near the United States. *J. Climate*, **19**, 3220–3236, doi:10.1175/JCLI3913.1.

Kaiser, H. F., 1960: The application of electronic computers to factor analysis. *Educ. Psychol. Meas.*, **20**, 141–151, doi:10.1177/001316446002000116.

Kallache, M., M. Vrac, P. Naveau, and P. A. Michelangeli, 2011: Nonstationary probabilistic downscaling of extreme precipitation. *J. Geophys. Res.*, **116**, D05113, doi:10.1029/2010JD014892.

Kalnay, E., and Coauthors, 1996: The NCEP/NCAR 40-Year Reanalysis Project. *Bull. Amer. Meteor. Soc.*, **77**, 437–470, doi:10.1175/1520-0477(1996)077<0437:TNYRP>2.0.CO;2.

Khaliq, M. N., and T. B. M. J. Ouarda, 2007: On the critical values of the standard normal homogeneity test (SNHT). *Int. J. Climatol.*, **27**, 681–687, doi:10.1002/joc.1438.

—, —, J. C. Ondo, P. Gachon, and B. Bobée, 2006: Frequency analysis of a sequence of dependent and/or non-stationary hydro-meteorological observations: A review. *J. Hydrol.*, **329**, 534–552, doi:10.1016/j.jhydrol.2006.03.004.

Kistler, R., and Coauthors, 2001: The NCEP–NCAR 50-Year Reanalysis: Monthly means CD-ROM and documentation. *Bull. Amer. Meteor. Soc.*, **82**, 247–267, doi:10.1175/1520-0477(2001)082<0247:TNNYRM>2.3.CO;2.

Knippertz, P., M. Christoph, and P. Speth, 2003: Long-term precipitation variability in Morocco and the link to the large-scale circulation in recent and future climates. *Meteor. Atmos. Phys.*, **83**, 67–88, doi:10.1007/s00703-002-0561-y.

Mann, H. B., 1945: Nonparametric tests against trend. *Econometrica*, **13**, 245–259, doi:10.2307/1907187.

Maraun, D., and Coauthors, 2010a: Precipitation downscaling under climate change: Recent developments to bridge the gap between dynamical models and the end user. *Rev. Geophys.*, **48**, RG3003, doi:10.1029/2009RG000314.

—, H. W. Rust, and T. Osborn, 2010b: Synoptic airflow and UK daily precipitation extremes. *Extremes*, **13**, 133–153, doi:10.1007/s10687-010-0102-x.

Martín-Vide, J., and J. A. Lopez-Bustins, 2006: The western Mediterranean oscillation and rainfall in the Iberian Peninsula. *Int. J. Climatol.*, **26**, 1455–1475, doi:10.1002/joc.1388.

Metropolis, N., A. Rosenbluth, M. Rosenbluth, A. Teller, and E. Teller, 1953: Equation of state calculations by fast computing machines. *J. Chem. Phys.*, **21**, 1087–1092, doi:10.1063/1.1699114.

Nasri, B., S. El Adlouni, and B. M. J. T. Ouarda, 2013: Bayesian estimation for GEV-B-spline model. *Open J. Stat.*, **3**, 118–128, doi:10.4236/ojs.2013.32013.

Olsen, J. R., J. R. Stedinger, N. C. Matalas, and E. Z. Stakhiv, 1999: Climate variability and flood frequency estimation for the upper Mississippi and lower Missouri Rivers. *J. Amer. Water Resour. Assoc.*, **35**, 1509–1524, doi:10.1111/j.1752-1688.1999.tb04234.x.

Ouachani, R., Z. Bargaoui, and T. B. M. J. Ouarda, 2013: Power of teleconnection patterns on precipitation and streamflow variability of upper Medjerda basin. *Int. J. Climatol.*, **33**, 58–76, doi:10.1002/joc.3407.

Ouarda, T. B. M. J., and S. El Adlouni, 2011: Bayesian non-stationary frequency analysis of hydrological variables. *J. Amer. Water Resour. Assoc.*, **47**, 496–505, doi:10.1111/j.1752-1688.2011.00544.x.

Pettitt, A. N., 1979: A non-parametric approach to the change-point problem. *J. Roy. Stat. Soc.*, **28C**, 126–135, doi:10.2307/2346729.

Poncet, J., 1970: La "catastrophe" climatique de l'automne 1969 en Tunisie. *Ann. Géogr.*, **79**, 581–595, doi:10.3406/geo.1970.15175.

Preisendorfer, R. W., 1988: *Principal Component Analysis in Meteorology and Oceanography*. Developments in Atmospheric Sciences, Vol. 17, Elsevier, 444 pp.

Raftery, A. E., and S. M. Lewis, 1992: Comment: One long run with diagnostics: Implementation strategies for Markov chain Monte Carlo. *Stat. Sci.*, **7**, 493–497, doi:10.1214/ss/1177011143.

—, —, 1995: The number of iterations, convergence diagnostics and generic Metropolis algorithms. *Practical Markov Chain Monte Carlo*, W. R. Gilks, D. J. Spiegelhalter, and S. Richardson, Eds., Chapman and Hall, 115–130.

Reiser, H., and H. Kutiel, 2011: Rainfall uncertainty in the Mediterranean: Time series, uncertainty, and extreme events. *Theor. Appl. Climatol.*, **104**, 357–375, doi:10.1007/s00704-010-0345-0.

Saidi, M. E. M., L. Daoudi, M. E. H. Aresmouk, and A. Blali, 2003: Rôle du milieu physique dans l'amplification des crues en milieu montagne montagnard: Exemple de la crue du 17 août 1995 dans la vallée de l'Ourika (Haut-Atlas, Maroc). *Sécheresse*, **14**, 107–114.

Sillmann, J., V. V. Kharin, X. Zhang, F. W. Zwiers, and D. Bronaugh, 2013: Climate extremes indices in the CMIP5 multimodel ensemble: Part 1. Model evaluation in the present climate. *J. Geophys. Res.*, **118**, 1716–1733, doi:10.1002/jgrd.50203.

Sterl, A., 2004: On the (in)homogeneity of reanalysis products. *J. Climate*, **17**, 3866–3873, doi:10.1175/1520-0442(2004)017<3866:OTIORP>2.0.CO;2.

Toreti, A., E. Xoplaki, D. Maraun, F. G. Kuglitsch, H. Wanner, and J. Luterbacher, 2010: Characterisation of extreme winter precipitation in Mediterranean coastal sites and associated anomalous atmospheric circulation patterns. *Nat. Hazards Earth Syst. Sci.*, **10**, 1037–1050, doi:10.5194/nhess-10-1037-2010.

Tramblay, Y., L. Neppel, and J. Carreau, 2011: Climatic covariates for the frequency analysis of heavy rainfall events in the Mediterranean region. *Nat. Hazards Earth Syst. Sci.*, **11**, 2463–2468, doi:10.5194/nhess-11-2463-2011.

—, W. Badi, F. Driouech, S. El Adlouni, L. Neppel, and E. Servat, 2012a: Climate change impacts on extreme precipitation in Morocco. *Global Planet. Change*, **82–83**, 104–114, doi:10.1016/j.gloplacha.2011.12.002.

—, L. Neppel, J. Carreau, and E. Sanchez-Gomez, 2012b: Extreme value modelling of daily areal rainfall over Mediterranean catchments in a changing climate. *Hydrol. Processes*, **26**, 3934–3944, doi:10.1002/hyp.8417.

—, S. El Adlouni, and E. Servat, 2013: Trends and variability in extreme precipitation indices over North Africa. *Nat. Hazards Earth Syst. Sci.*, **13**, 3235–3248, doi:10.5194/nhess-13-3235-2013.

Vicente-Serrano, S. M., S. Beguería, J. I. López-Moreno, A. M. El Kenawy, and M. Angulo-Martínez, 2009: Daily atmospheric circulation events and extreme precipitation risk in
northeast Spain. Role of the North Atlantic Oscillation, the western Mediterranean oscillation, and the Mediterranean oscillation. *J. Geophys. Res.*, **114**, D08106, doi:10.1029/2008JD011492.

Vrac, M., and P. Naveau, 2007: Stochastic downscaling of precipitation: From dry events to heavy rainfalls. *Water Resour. Res.*, **43**, W07402, doi:10.1029/2006WR005308.

Wald, A., and J. Wolfowitz, 1940: On a test whether two samples are from the same population. *Ann. Math. Stat.*, **11**, 147–162, doi:10.1214/aoms/1177731909.

Wanner, H., S. Brönnimann, C. Casty, D. Gyalistras, J. Luterbacher, C. Schmutz, D. B. Stephenson, and E. Xoplaki,

2001: North Atlantic Oscillation—Concepts and studies. *Surv. Geophys.*, **22**, 321–382, doi:10.1023/A:1014217317898.

Wetterhall, F., S. Halldin, and C. Y. Xu, 2005: Statistical precipitation downscaling in central Sweden with the analogue method. *J. Hydrol.*, **306**, 174–190, doi:10.1016/j.jhydrol.2004.09.008.

Wilcoxon, F., 1945: Individual comparisons by ranking methods. *Biom. Bull.*, **1**, 80–83, doi:10.2307/3001968.

Xoplaki, E., F. J. González-Rouco, J. Luterbacher, and H. Wanner, 2004: Wet season Mediterranean precipitation variability: Influence of large-scale dynamics and trends. *Climate Dyn.*, **23**, 63–78, doi:10.1007/s00382-004-0422-0.