An objective cross-validation framework for mapping rainfall hazard based on rain gauge data

We propose an objective framework for estimating rainfall cumulative distribution function within a region when data are only available at rain gauges. Our methodology is based on the evaluation of several goodness-of-fit scores in a crossvalidation framework, allowing to assess goodness-of-fit of the full distribution but with a particular focus on its tail. Crossvalidation is applied both to select the most appropriate statistical distribution at station locations and to validate the mapping of these distributions. Our methodology is applied to daily rainfall in the Ardèche catchment in South of France, a 2260 km 5 catchment with strong disparities in rainfall distribution. Results show preference for a mixture of Gamma distribution over seasons and weather patterns, with parameters interpolated with thin plate spline across this region. However the framework presented in this paper is general and could be likewise applied in any region, with possibly different conclusion depending on the subsequent rainfall processes.


Introduction
In recent years, Mediterranean storms involving various spatial and temporal scales have hit many locations in southern Europe, causing casualties and damages (Ramos et al., 2005;Ruin et al., 2008;Ceresetti et al., 2012a).Assessing the frequency of occurrence of extreme rainfall in a region is usually done by the computation of return level maps.This requires relating any (large) amount of rainfall at a given location to its return level, i.e. to the frequency such an amount is expected to occur on average at this location.In other words, it requires knowing the cumulative distribution function (CDF) of extreme rainfall at any grid point of the map.However there are other situations when not only the largest rainfalls are of interest, but also smaller and even zero rainfall values.This is for example the case in rainfall simulation frameworks, e.g. when rainfalls are input of spatially distributed hydrological models.In such a case one needs to be able to simulate any possible rainfall fields.This implies knowing both the local occurrence of any rainfall value with the right frequency, and not only the largest ones, and their spatial co-occurrence.Other domains include the evaluation of numerical weather simulations (e.g.Froidurot et al., 2016) or the investigation of the climatology of rainfall events in a region.
A difficulty in producing rainfall return level maps is that knowing the CDF at any grid point ideally requires observing rainfall on a grid scale.However long-enough gridded data with good-enough quality is often lacking.Radar and satellite estimations are usually available for about 10 years at best, and only for selected regions.In addition rainfall estimation in complex topography is particularly tricky, e.g.due to the mountain ranges shielding the radar beam (Germann et al., 2006), or to the complex relationship between satellite-measured radiances and rainfall reaching the ground (Tian and Peters-Lidard, 2010).On the other hand, rain gauge networks are usually operational for 50 to 100 years in the main part of the world, at least at daily scale, but they only provide point observations.Thus, two main methods are usually adopted for learning the CDF of rainfall at any location when observations are only available at selected locations.The first one resorts to the spatial interpolation of point data supplied by rain gauges.This allows transforming point observations into gridded ones, and so to estimate gridded CDFs of rainfall.Among the most performing methods for spatial interpolation of daily rainfall are kriging, Inverse Distance Weighting and spline-surface fitting (e.g.Camera et al., 2014;Creutin and Obled, 1982;Goovaerts, 2000;Ly et al., 2011;Rogelis and Werner, 2013).In complex topography, there may be some gain in applying these methods locally, e.g.considering local precipitation altitude gradients (Frei and Schär, 1998;Gottardi et al., 2012;Lloyd, 2010).However none of the above statistical methods are able to fully account for the statistical properties of rainfall fields.A first difficulty is due to the presence of zeros, which complicates interpolation and can lead to negative interpolated rainfalls -although this could be partially overcome by using analytical transformation of the raw variable.A second difficulty is that rainfall distribution is usually heavy tailed and interpolation methods, by smoothing values, lack quality for representing the most extreme events (Delrieu et al., 2014).
A second way of mapping rainfall hazard is, rather than interpolating the point observations, to map the parameters of CDFs fitted on rain gauge series.In addition to the choice of interpolation models, comes now the choice of the marginal model of rainfall amounts on wet days (referred as nonzero rainfalls).The most commonly used CDFs at daily scale include the exponential, Gamma, lognormal, Pareto, Weibull and Kappa models (Papalexiou et al., 2013).Noting that these distributions tend to underestimate extreme rainfall amounts (Katz et al., 2002), a recent flurry of research developed hybrid models based on mixtures of distributions for low and heavy amounts (Vrac and Naveau, 2007;Furrer and Katz, 2008;Li et al., 2012).
More recently Naveau et al. (2016) proposed a family of distributions that is able to model the full spectrum of rainfall, while avoiding the use of mixtures of distributions.Several studies compared marginal models for rainfall (e.g.Mielke and Johnson, 1974;Swift and Schreuder, 1981;Cho et al., 2004;Husak et al., 2007;Papalexiou et al., 2013), but focusing usually on a couple of CDFs.Other studies compared methods for mapping rainfall hazard, and particularly extreme rainfall, assuming a given CDF (Beguería and Vicente-Serrano, 2006;Beguería et al., 2009;Szolgay et al., 2009;Blanchet and Lehning, 2010;Ceresetti et al., 2012b).However there is, to the best of our knowledge, no study assessing goodness-of-fit of the full procedure of rainfall hazard mapping, i.e. from marginal modeling to the production of hazard maps.
Our study aims at filling this gap by proposing an objective cross-validation framework that is able to validate the full procedure of rainfall hazard mapping starting from point observations.Our framework features three characteristics: i) it selects both the marginal and mapping models, ii) it validates the full spectrum of rainfall, from small to long-term extrapolated amounts, iii) it applies on a regional scale.The framework is illustrated on the Ardèche catchment in South of France.Despite its relatively small size, this test case is particularly challenging as showing extraordinarily strong disparities in rainfall statistics in very short distance.Following previous studies in the region (Evin et al., 2016;Garavaglia et al., 2010Garavaglia et al., , 2011;;Gottardi et al., 2012), the compared marginal distributions involve seasonal and weather pattern subsampling, considering different models for the subclass-dependent distributions.However the proposed cross-validation framework is general, as involving objective criteria, and could likewise be used to select among any other distribution.Section 2 presents the data.Section 3.1 and 3.2 describe the marginal distributions and mapping models considered in this study and present the cross-validation scores of model selection.Section 3.3 detail the procedure of model selection from marginal modeling to hazard mapping.Section 4 gives extensive results for the Ardèche catchment.Section 5 concludes.

Data
We illustrate our framework on the Ardèche catchment (2260 km 2 ) located in South of France (see Figure 1).The region includes part of the south-eastern edge of the Massif Central, where the highest peaks of the region are located (more than 1500 m.a.s.l), and the lower elevated Rhône valley (down to 10 m.a.s.l).The south-eastern slope of the Massif Central is known to experience most of the extreme storms and resulting flash floods (Figure 2 of Nuissier et al., 2008).These so-called "Cévenol" events are produced by quasi-stationary mesoscale convective systems that stabilize over the region during several tens of hours.The positioning and stationarity of these systems are largely influence by the topography of the surrounding mountain massifs (Nuissier et al., 2008).We use two daily rain gauge networks maintained respectively by Electricité de France and Météo-France.We consider the 15 rain gauges inside the catchment, together with the 27 stations located less than 15km outside.This gives a total of 42 stations with 20 to 64 years of data between January 1, 1948 and December 31, 2013.
In both databases, daily values are recorded every day at 6AM UTC, corresponding to rainfall accumulation between 6AM of the previous day and 6AM of the present day.
The Ardèche catchment is chosen for illustration purpose and because, despite its relatively small size, it shows very large disparities in rainfall distribution.To illustrate these disparities, we show in Figure 2 the averages of annual totals and annual maximum daily rainfalls for each station.Computing the ratios between the largest and lowest values in Figure 2 gives a ratio of 2.6 for the annual totals and 3.2 for the annual maxima.For comparison the latter ratio is barely lower than the ratio found over the whole of France, which amounts 4. For both annual totals and annual maxima, the strongest values in the region concentrate along the Massif Central ridge, while much smaller values are found a few km apart in the Massif Central plateau or in the Piémont.Concentration of daily rainfall and particularly of extreme daily rainfall along the Massif Central ridge has already been documented in many studies, see e.g. Figure 10 of Blanchet et al. (2016a).We assume in this study temporal stationarity of rainfall.Case of potential nonstationarity due to climate change will be discussed in Section 5.

Considered marginal models
Let R be the random variable of daily rainfal amount at a given station.R is zero with probability p 0 and, for any r > 0, we have the following decomposition: where G is the CDF of nonzero rainfall at the considered station.Choice of G is an issue.One of the difficulties is that we wish to model adequately both the bulk of the distribution of nonzero rainfall and its tail, i.e. the probability of extreme rainfall to

Distribution
CDF G(r) or density g(r), for r > 0 Parameters Gamma 1. Considered models for the marginal distributions of nonzero rainfall.Γ in the Gamma density is the complete Gamma function occur.The most common models for nonzero rainfall include the Gamma, Weibull and lognormal models (Papalexiou et al., 2013), whose CDF G(r) or densities g(r) = ∂G(r)/∂r, r > 0 are given in Table 1.Although less common, another family of models for nonzero rainfall relies on univariate extreme value theory, which tells that probabilities of the form pr(R ≤ r|r > q), with q large, can be approximated by either an exponential or a Generalized Pareto tail (Coles, 2001, chapter 4).This led Naveau et al. (2016) to propose the extended exponential and extended Generalized Pareto distributions, whose CDF are given in Table 1.Note that less parsimonious models are given in Naveau et al. (2016) but they are not be considered in the present study.The extended exponential and extended Generalized Pareto distributions of Table 1 insure that the occurrence probability of small (but nonzero) rainfall amounts is driven by κ while the upper tail of nonzero rainfall is equivalent to a Generalized Pareto tail.
The extended exponential model is also called "Generalized exponential" and it has been used previously for extreme rainfall in Madi and Raqab (2007); Kaźmierczak and Kotowski (2015).
In the models of Table 1, rainfall is implicitely assumed to come from a single distribution.This assumption may be questioned.Indeed, different climatological processes trigger precipitation, leading to the occurrence of rainfall of different natures and intensities (e.g.convective vs. stratiform precipitation).Furthermore, rainfall occurrence and intensities often vary with season, reflecting both variations in temperature and in storm tracks, for example.For this reason, Garavaglia et al. (2010) proposed the use of subsampling based on seasons and weather patterns (WP).Each day of the record period is assigned to a WP.If S seasons and K WP are considered, then days are classified into S × K subclasses.The law of total probability gives, for all r > 0, where p s,k is the probability that a given day is in season s and in WP k (thus s k p s,k = 1).Following (1), R in season s and WP k is zero with probability p 0 s,k and, for any r > 0, we have the decomposition where G s,k is the CDF of nonzero rainfall at the considered station for a day in season s and WP k.This gives in (2), for all r > 0, where p 0 = S s=1 K k=1 p s,k p 0 s,k is the probability of any day to be dry.Nonzero precipitation amounts defined by (3) have CDF: where p s,k = p s,k (1−p 0 s,k )/(1−p 0 ).(4) defines a mixture of S ×K distributions, e.g. a mixture of S ×K Gamma distributions.
Analogously, the CDF of nonzero precipitation amounts in a given season s writes where A similar idea is used in Wilks (1998) for example, but considering a mixture of 2 (exponential) distributions is an unsupervised way, i.e. without relying on a priori subsampling.It shows the advantage of not requiring prior knowledge on the classification but it is at the same time more difficult to estimate, in particular if the models for different seasons and WP do not differ much.
In this article, we will consider the supervised case (3), with S = 2 seasons and K = 3 WP, considering the five models of Table 1 for the distribution of precipitation amounts G s,k (see Section 3.3).This implies that estimation of G s,k can be made independently of each others, by considering only the days of the record belonging to season s and WP k.A variety of inference methods exists.For rainfall analysis, two options are popular: maximum likelihood (ML) estimation and a method of moments based on probability weighted moments (PWM).However, as noted in Naveau et al. (2016), ML estimation may fail for rainfall because the discretization due to instrumental precision strongly affects low values, which biaises ML estimation if not accounted for.One way to circumvent this issue is to resort to censored likelihood but choice of the censoring threshold is in itself an issue.Results on our data (not shown) reveal that the threshold has to be no smaller than 5mm.PWM, on the other side, is much more robust against discretization since it is based on summary statistics, rather than on the exact values of observations (Naveau et al., 2016).For this reason, we estimate in this study the distributions of precipitation amounts G s,k by PWM, while p 0 s,k at a given station is estimated empirically as p0 s,k = d 0 s,k /d where d is the number of observations and d 0 s,k is the number of zero values in season s and WP k.Combining estimations in (3) gives an estimation of the rainfall CDF at the considered station, and in (4) an estimation of the CDF of nonzero rainfall.
Estimates of return levels are then obtained as follows.The T -year return level r T is the level expected to be exceeded on average once every T years.It satisfies the relationship pr(R ≤ r T |R > 0) = 1 − 1/(T δ) where δ is the mean number of nonzero rainfall per year at the considered station.When subsampling (4) is considered, there is not an explicit formulation and estimation of r T is obtained numerically by solving pr(R ≤ r T |R > 0) = 1 − 1/(T δ) in (4).

Evaluation at regional scale in a cross-validation framework
The goal of this evaluation is to assess which marginal model performs better at the regional scale, i.e. for a set of n stations taken as a whole, rather than individually.We follow the split sample evaluation proposed in Garavaglia et al. (2011) and Renard et al. (2013).We divide the data for each station i into two subsamples, C (1) i and C (2) i , and consider nonzero rainfall for these two subsamples.We fit a given competing model on each of the subsamples, giving two estimated distributions of G in (4): Ĝ(1) i , estimated on C (1) i , and Ĝ(2) i , estimated on C (2) i .Our goal is to test the consistency between validation data and predictions of the estimates, both for the core and tail of the distributions, and the stability of the estimates when calibration data changes, focusing particularly on the tail which is usually less stable.
As shown in Table 2, three families of scores are computed, assessing respectively i) accuracy of the estimations along the full range of observations (MEAN(NRMSE)), ii) reliability of the tail of the estimated distribution, checking in particular systematic over or under-estimation of the observations (AREA(FF ) and AREA(N T )), and iii) stability of the tail at extrapolation (MEAN(SPAN)).The scores relating the tail of the distribution have been proposed and used in Garavaglia et al. (2011); Renard et al. (2013); Blanchet et al. (2015).In the split sample evaluation framework, four scores can be derived of a given score Sc: Sc (12) is the regional score when G (2) i is validated on the nonzero rainfall subsample C (1) i .Sc (21) , Sc (11) and Sc (22)  are obtained symmetrically.Sc (11) and Sc (22) are calibration scores, while Sc (12) and Sc (21) are cross-validation scores.For the sake of conciseness, we detail below the case of Sc (12) for the different scores.
The NRMSE (Normalized Root Mean Squared Error) evaluates reliability of the fits in the whole observed range of nonzero rainfall, by comparing observed and predicted return levels of daily rainfall.For a given station i ∈ {1, . . ., Q}, where n (1) i is the number of nonzero rainfall in C (1) i for station i, T k ranges the observed return periods of nonzero rainfall in i,T k is the observed daily rainfall associated to the return period T k for the subsample C (1) Without loss of generality we assume T 1 , . . ., T n (1) i to be sorted in descending order (so T 1 is associated to the maximum over C (1) i ).If station i has δ i nonzero rainfall per year on average, usual practice is to consider the kth largest return period as i , and to estimate r (1) i .For the set of Q stations, we obtain a vector of NRMSE (12) of length Q which should remain reasonably close to zero.A regional score is obtained by computing the mean of the Q values: For competing models, the closer the mean is to 0, the better the goodness-of-fit.
NRMSE assesses goodness-of fit of the whole distribution in the observed range.Now let have a closer look at the tail of the distribution, and in particular at the maximum over C (1) i,T1 in (6), that for shortness we denote m (1) i .If G i is the true distribution of nonzero rainfall, then the corresponding random variable M (1) i has distribution G i to the power n (1) i , whose variance is large.Thus computing error based on the single realization m (1) i would be very uncertain.For this reason, Renard et al. (2013) proposed to make evaluation by pulling together the maxima of the Q stations, after transformation to make them on the same scale.It is based on the idea that if X has CDF F , then F (X) follows the uniform distribution on (0, 1).Taking i should be a realization of the uniform distribution.For the set of Q stations, this gives a uniform sample ff (12) of size Q.
Hypothesis testing for assessing the validity of the uniform assumption is challenging because the ff (12) i are not independent from site to site, due to the spatial dependence between data.Thus Blanchet et al. (2015) proposed to base comparison on the divergence of the density of the ff (12) to the uniform density.A reasonable estimate of the latter is obtained by computing the empirical histogram of the ff (12) with 10 equal bins between 0 and 1.As illustrated in Figure 3 12) should be reasonably uniform on (0, 1).
i ) tends to overestimate the true G i ), or in other words the return period of the maximum over C (1) i tends to be underestimated (case of over-estimated risk).If the histogram is right-skewed, the return period of the maximum over C (1) i tends to be over-estimated (case of under-estimated risk).Although any scenario of misfitting could theoretically be possible, in practice the histograms of ff (12) show mainly the three above alternatives: either a good fit (flat histogram), or a tendency towards a systematic under-or over-estimation (left-or right-skewed histograms).By focusing of maximum values, the histogram of ff (12) can be seen as a way of assessing systematic bias in the very tail of the distribution.For a more quantitative assessment, we compute the area between the density of the ff (12) and the uniform density as follows: AREA(FF (12) ) = 1 18 where # is the number of elements of the set.The term inside the absolute value in ( 8) is the difference between densities in the cth bin.The division by 18 forces the score to lie in the range (0, 1) with lower values indicating better fits (the worst case being all values lying in the same bin).Figure 3 shows that, when Q = 42 stations are considered, a value of AREA(FF (12) ) around 0.2 corresponds to no systematic bias in the very tail of the distribution at regional scale, whereas a value around 0.5 corresponds to a strong over-or under-estimation.In the latter case, only looking at the histogram can inform on whether overor under-estimation applies.
The N T criterion is an alternative to FF assessing reliability of the fit of the tail but focusing on prescribed (large) quantiles rather than on the overall maximum.It applies the same principle as FF , involving a transformation of X to F (X), but considering X as K (1) i,T , the random variable of the number of exceedances over C (1) i of the T -year return level, i.e.K (1) should be a realization of the discrete uniform distribution.Randomisation to transform n ).Taking T as e.g.half to one-quarter the length of the observations allows to assess reliability of the close tail of the distribution.As such, it is a good complement to FF that focuses on the far tail (i.e. on the maximum).
Last but not least, the SPAN criterion evaluates the stability of the return level estimation, when using data for each of the two subsamples.More precisely, for a given return period T and station i, where r(1) i,T , e.g., is the T -year return level for the distribution G estimated on subsample C (1) i of station i, i.e. such that Ĝ(1) i {r (1) is the relative absolute difference in T -year return levels estimated on the two subsamples.
It ranges between 0 and 2; the closer to 0, the more stable the estimations for station i.For the set of Q stations, we obtain a vector of SPAN T of length Q with a distribution which should remain reasonably close to zero.A rough summary of this information is obtained by computing the mean of the Q values of SPAN i,T , i = 1, . . ., Q: For competing models, the closer the mean is to 0, the more stable is the model.When T is larger than the observed range of return periods, MEAN(SPAN T ) evaluates the stability of the return levels in extrapolation.Note that it is by definition 0 in calibration and thus it is only useful in cross-validation.
For the sake of concision, in the rest of this article the scores MEAN(NRMSE), AREA(FF ), AREA(N T ) and MEAN(SPAN T ) will be referred to as the NRMSE, FF , N T and SPAN T scores.

Considered mapping models
Let R i be the random variable of daily rainfall at station i, i = 1, . . ., Q. Applying Section 3.1 at station i gives an estimate Ĝi (r) of the CDF G i (r) = pr(R i ≤ r|R i > 0).Our goal is to derive an estimate of the CDF of nonzero daily rainfall at any location l of the region, pr(R(l) ≤ r|R(l) > 0), based on the Q estimated CDFs Ĝi (r).Location l refers here to the three coordinates of ground projection coordinates and altitude, that we write l = (x, y, z).Let θi be the set of estimated parameters for station i and θi,j its jth element.θi is composed of the S × K probability of zero rainfall p 0 s,k and the 2 × S × K or 3 × S × K parameters of the distributions G s,k , depending on the marginal distribution (see Table 1).We assume the θ i,j ordered so that the first S × K elements are the p 0 s,k .We aim at estimating the surface response θ j (l) at any l of the region, knowing θ j (l i ) = θi,j .In this study we consider three of the most popular method: kriging interpolation, linear regression methods and thin plate spline regressions.However the parameters θ j s are constrained whereas these models apply the unbounded variables: the probabilities p 0 s,k lie in (0, 1), while the parameters of Table 1 are all positive.Therefore we apply the mapping models to transformations of θ j , i.e. to ψ j = tr(θ j ) where tr maps the range of values of θ j to (−∞, +∞).In this study we consider ψ j (l) = Φ{θ j (l)} if j ≤ S × K (i.e. if θ j is any p 0 s,k ), where Φ is the standard Gaussian CDF, and to ψ j (l) = log{θ j (l)} otherwise.Other transformations would be possible, in particular p 0 s,k may be transformed with the logit function, but will not be considered here for the sake of concision.Thus we aim at estimating ψj (l) given values ψ j (l i ) = ψi,j at station locations, with obvious notations.If l ≤ S × K, estimates of θ j (l) are then obtained as θj (l) = Φ −1 ( ψj (l)).Otherwise surface response estimates are obtained as θj (l) = exp( ψj (l)).For the sake of clarity, we omit below the index j, considering a surface ψ(l) to be estimated for all l in the region, given values ψ(l i ) = ψi .
The considered mapping models are listed in Table 3.Three families of method are considered: kriging, linear regression and thin plate spline.Additionally to how they map values, there is a fundamental difference between these models: kriging is an exact interpolation, i.e. ψ(l i ) = ψi at any station location l i used to estimate the model.On the contrary the linear regression models and thin plate splines provide inexact interpolations: in the great majority of the time, ψ(l i ) = ψi (the goal being obviously to minimize the overall error).
For the kriging interpolation, cases with and without external drift are tested (chapter 3.6 of Diggle and Ribeiro, 2007).
The external drift, if any, is modeled as a linear function of altitude (i.e. of the form a 0 + a 1 ζ), considering ζ as either the altitude of the station (z), or, following Hutchinson (1998), as a smoothed altitude (Z) derived by smoothing a 1km Digital Elavation Model (DEM) with 5km moving windows (i.e.taking Z as the average altitude of 25 DEM grid points).The results that will be presented below correspond to the case of an exponential covariance function of the form ρ(h) = e −h/β , with β > 0. We also considered the case of a powered exponential covariance function ρ(h) = e −(h/β) ν with 0 < ν ≤ 2, but this resulted in a slight loss of stability due to the additional degree of freedom, without improving the accuracy at regional scale.
For the sake of concision, these results are not presented here.Combining alternatives for the drift part gives a total of three kriging interpolation models with 2 to 3 unknown parameters each, for each ψ.Estimation of the kriging models are made by maximizing the likelihood associated to the ψi , assuming that ψ(l) is a Gaussian process (see chapter 5.4 of Diggle and Ribeiro, 2007).Alternatives are to estimate drifts and variograms by least squares in different steps, with the risk of biaising estimates (chapter 5.1 to 5.3 of Diggle and Ribeiro, 2007).Both estimation methods are available in the R package 'geoR' (e.g. functions 'likfit' and 'variofit').In the case without drift, prediction at any site l of the region is obtained as where the weights w i (h i ) are derived from the kriging equations and satisfy The weights depend on the estimated covariance function and on the distance h i between l and station location l i in the (x, y) space (i.e. ).In the case with external drift, prediction at any location l of the region is then obtained as where ζ is the altitude at location l (i.e.either z or Z).Predictions ( 11) and ( 12) are exact: ψ(l i ) = ψi , and consequently θ(l i ) = θi .
For the linear regression models, we start from regressions of the form ψ(l) = a 0 +a 1 x+a 2 y +a 3 ζ +a 4 x 2 +a 5 y 2 +a 6 xy + a 7 xζ + a 8 yζ + (l), where (l) ∼ N (0, σ 2 ) and ζ is, as before, either the altitude of the station (z) or the smoothed altitude (Z).We consider Akaike Information Criteria (AIC), defined as AIC= 2η − 2 log L, where η is the number of parameters (10 at the start) and L is the maximum likelihood value of the regression model.The lower AIC, the better the model.Then we repeatedly drop the variable that increases most the AIC -if any-, and add the variable that decreases most the AIC -if any.This stepwise method is implemented in the 'step' function of the R package 'stats'.At algorithm stop, the model may contain 1 to 10 parameters, for each ψ.Predictions θ(l) are then obtained as the back-transformation the estimated regressions.
These methods are implemented in the function 'Tps' of the R package 'fields'.In the bivariate case, ψ(l) is modeled as where λ is the so-called smoothing parameter, which controls the trade-off between smoothness of the estimated function and its fidelity to the observations.It can be estimated by generalized cross-validation.Then predictions are obtained as where h i is the Euclidean distance in the (x, y) space between l and station location l i .The partial trivariate case assumes that ψ(l) − a 3 ζ is a bivariate thin plate spline, where a 3 is fixed and ζ is either z or Z.To make the connection with kriging, ψ(l) can thus also be seen as a bivariate thin plate spline with (fixed) drift in ζ.The coefficient a 3 is estimated in a preliminary step by regressing ψi against ζ i .Estimation of the bivariate thin plate spline for ψ(l) − a 3 ζ is made as described above given the values of ψi − a 3 ζ i .Predictions are obtained as Finally in the trivariate case, we have ψ(l) = u(x, y, ζ) + (x, y, ζ).The minimization problem is similar to (13) with a penalization enlarged by several terms (Wahba and Wendelberger, 1980).Predictions are then obtained as where h i is the Euclidean distance in the (x, y, ζ) space between l and station location l i , scaling the altitude by a factor of 10 following Boer et al. (2001); Hutchinson (1998 ). Coefficients a i and b i in ( 14) to ( 16) are estimated by solving a linear system of order Q involving the smoothing parameter λ.Note that the trivariate case (16) differs from the bivariate case with drift (15) in two ways.First, (16) considers distance in the (x, y, ζ) space whereas (15) considers distance in the (x, y) space.Second, in (16) the weights associated to the stations are linear functions of the distance, unlike in (15) (see the term h 2 i log(h i ) vs. h i ).

Evaluation at regional scale in a cross-validation framework
Evaluation is performed in two ways.The first one is a leave-one-out cross-validation scheme aiming to test at regional scale how the interpolated distributions are able to fit the data of the stations when these data are left out for estimating the mapping model.The second step assesses spatial stability by comparing the interpolated distributions obtained at a given station whether the data of this station are used or not in the mapping estimation.In other words, it is a comparison between leave-one-out and leave-zero-out mappings.So the two evaluations differ in that the first one compares an interpolated distribution to data, while the second step compares two interpolated distributions.
First, let consider a given parameter estimate θ(1) i,j obtained at station i based on the subsample C (1) i (for a given marginal model).We apply a leave-one-out cross-validation scheme: for each station i 0 alternatively, we use the set of θ(1) i,j , i = i 0 to estimate the surface response θ(1) j (l).We store the value of this estimate at station location i 0 , i.e. θ(1) j (l i0 ) = θ(1) i0,j Repeating this for every parameter θ j gives estimation of the full set of parameters at station i 0 , i.e. estimation G(1) i0 of G i0 .Iterating over the stations, we obtain new estimates Q of G 1 , . . ., G Q .Applying similarly for the subsample C (2) Q .We can assess reliability of these estimates at regional scale by computing the scores Sc (11) , Sc (22) , Sc (12) and Sc (21) of Section 3.1.2,where Ĝi is replaced by Gi .Note all these scores are cross-validation scores since the estimates G(1) 1 and G(2) 1 are computed without using any data of station i.Second, we consider the set of all θ(1) i,j , i = 1, . . ., Q to estimate the surface response θ(1) j (l).We store the value of this function at every station location, giving new estimates G * (1) 1 , . . ., G * (1) Q .Note that in the particular case of kriging, G * (1) i is exactly Ĝ(1) i since it is an exact interpolation method, so every θ(1) i,j equals θ * (1) i,j .We can assess the stability of the interpolated distributions at a given location l i when observations are available or not at this location by comparing G * (1) i (r) and G(1) i (r) for all r.For this we discretize r between 0 and 450mm (which is the overall maximum rainfall) with 1mm step and we compute the total variation distance (TVD) between G * (1) i , which are given by: TVD where e.g.

g(1)
i is the density function associated to G(1) i .Note that the KLD is not symmetric.Written as such, it can be interpreted as the amount of information lost when G(1) i as the "true" distribution of data.TVD and KLD differ in that TVD focuses on the largest deviation between the two CDFs, whereas KLD somewhat integrate deviations along rainfalls.Obviously, one would like the interpolation to be as stable as possible when data are available or not at station i, i.e. that the lower TVD i and KLD i , the more stable the interpolation at station i.

Regional scores MEAN(TVD
(1) i ) and MEAN(KLD i ) are obtained similarly for the subsample C (2) .For competing models, the closer the means are to 0, the more spatially stable is the interpolation.For shortness, we will refer to MEAN(TVD) and MEAN(KLD) as the TVD and KLD scores, respectively (Table 2).

Procedure of model selection at regional scale
We wish to evaluate and compare performance of both marginal and mapping models for estimating rainfall frequency across the region.We consider models both with and without season/WPs.For the cases involving the use of WPs, we use the WP classification described in Garavaglia et al. (2010), which is obtained by clustering synoptic situations (geopotential heights) for France and surrounding areas into eight classes.However a grouping of the eight WPs into three is made in order to improve the robustness of the method while conserving the diversity of the rainy synoptic situations.The choice of the grouped WPs is made in a separate analysis based on the spatial correlation of rainfall in the different WPs.The range of spatial correlation is twice as big in WP1 than in WP2, and three times as big in WP1 than in WP3.The occurrence statistics of the three WPs for the period 1948-2013 are presented in Figure 4.The yearly occurrence of the three WPs is roughly similar (27% for WP1, 36% for WP2, 37% for WP3).However the WPs show very different seasonality.In particular WP1 is more frequent in spring and autumn, which correspond to wetter periods, particularly in autumn (see the monthly averages of nonzero rainfall in Figure 4).
WP3 is more frequent in summer, which is the dryest season, while WP2 features almost a reversed seasonality compared to WP3.This shows that, although being based on the spatial dependence, the WPs are linked to the seasonality of rainfall in the region.
In cases where subsampling is also undertaken by season, we impose a restriction of S being two seasons, representing the season-at-risk during which most of the annual maxima are observed, and the season-not-at-risk.Furthermore, we impose the season-at-risk to be the same for all the stations due to the little extent of the region.Based on Figure 4, we define the season-at-risk as the three months of September, October and November, as in Garavaglia et al. (2010); Evin et al. (2016) for example.Alternative for bigger regions would be to select the months composing the season-at-risk following the procedure described in Blanchet et al. (2015).

Marginal selection procedure.
The full cross-validation procedure for selecting both the marginal and mapping models is summarized in Figure .First we consider the marginal distributions of Table 1 and select the best of them at regional scale, as described in Section 3.1.2: 1. We divide the days of 1948-2013 into two subsamples of equal size, denoted C (1) and C (2) .Given the weak temporal dependence of rainfall in the region (80% of the wet periods have length lower than 3), division is made by randomly choosing blocks of five consecutive days to compose C (1) , the remaining blocks composing C (2) .2. For every station i, we consider the set of observed days in C (1) and C (2) , giving C (1) i and C (2) i .4. We compute the scores of Section 3.1.2,getting two calibration scores (( 11) and ( 22)) of NRMSE, FF and N T and two cross-validation scores (( 12) and ( 21)) of NRMSE, FF , N T and SPAN T .For N T , we consider T = 5 years, which is lower than the minimum length of the calibration data and allows focusing on the tail but still having several exceedances of the T -year return level at every station.So FF , by focusing on the maximum of roughly 10 to 30 years of data, can be seen as an evaluation score for the far tail, while N 5 can be seen as a evaluation score for the close tail.For SPAN T , we consider T = 100 and T = 1000 years in order to test extrapolation far in the tail but at a scale still commonly used for engineering purposes (dam building, protections etc, Paquet et al., 2013).

We fit each distribution of
5. We repeat 50 times the steps 1-4.
We obtain 100 values of each calibration score and 100 values of each cross-validation score.We apply this procedure to the four distributions of Table 1, considering the four alternatives: no season nor WP (S = 1, K = 1), two seasons but no WP (S = 2, K = 1), no season but three WPs (S = 1, K = 3), two seasons and three WPs (S = 2, K = 3).Comparing the distributions of the scores of the 16 models allows us the select the marginal distribution yielding to the best fit at regional scale.We select this marginal model for further consideration.

Mapping selection procedure.
Second we consider the mapping models of Section 3.2.1 for interpolating the selected marginal model, and we select the best of them in two ways, as described in Section 3.2.2.
2. We estimate the mapping models of Section 3.2.1 following the leave-one-out cross-validation framework of Section 3.2.2.We obtain new estimates G(1,t) i for each station i and each mapping model.Each i and G (2) i since the computation of G(1,t) i did not use any data of station i.
3. We compute the scores of Section 3.1.2associated to G(1,t) i , i = 1, . . ., Q.We obtain for each score two values (e.g. (11)and FF (21) when considering G(1,t) i and the maximum value over either C (1) or C (2) ).All these score are cross-validation scores.

FF
4. We estimate the mapping models of Section 3.2.1, using all the stations to make interpolation.We obtain new estimates G * (1,t) i for each station i and each mapping model.NRMSE scores) than for its tail (FF and N 5 ).Reliability scores are even better when both seasons and WPs are considered (case (2,3)), particularly for the tail of the distribution.
Obviously, there is a loss of stability when considering seasons and/or WPs due to the increased number of parameters.
However the score of SPAN 100 ranges 0.08-0.14, which means that the two estimates of the 100-year return levels over C (1)   and C (2) differ by 8 to 14%, which seems acceptable.

5
We illustrate the quality of the fit for the station Antraigues, located in the very foothills of the Massif Central slope (see Figure 1).We focus on the tail of the distribution by looking at the return level plot (here beyond the yearly return period).
Antraigues is chosen to illustrate the case of an heavy-tailed distribution of rainfall: its 20-year return level is about 3 times larger than its yearly return level.Of course, some variability is found in the return level estimations depending on the subsample used for estimation.The values of (S, K) are indicated in the title.The dotted lines show the 95%-envelope of return level estimates over the 100 subsamples.
The plain line shows the median estimates.The gray points show the full sample (35 years).Each estimation is based on half of these points.intervals (that could be obtained by bootstrapping for example) but variability when only half the data in used from calibration.
Figure 7 clearly shows that considering seasons or WPs allows to get heavier-tailed distributions.The median estimate follows quite closely the empirical points.Variability is relatively low, altough it is slightly larger than with the other marginal models involving less parameters.Coefficient of variation of the 100-year return level is less than 7%, in coherence with the SPAN 100 of Figure 6 at regional scale.
Due to its better fit for the Gamma model (Figures 6 and 7) as for the other distributions (not shown), we select the mixture model with S = 2 seasons and K = 3 WPs for further investigation.Figure 8 shows the scores of cross-validation when the parent distribution G s,k is either the extended exponential, extended Generalized Pareto, Gamma, lognormal or Weibull distribution.The reliability scores NRMSE, FF and N 5 in the lognormal case are missing because they lie far above the upper range of the depicted values (e.g. the median NRMSE is about 0.7), which clearly rules out the use of the lognormal model for this region.The reliability criteria of the four other distributions all show the same pattern: a better performance of the Gamma model, closely followed by the extended exponential case.Then comes the extended Generalized Pareto, itself closely followed by the Weibull model.A closer look at the values of ff i and n i,5 for all stations and samples reveals that the weaker reliability of the Weibull and extended Generalized Pareto models is due to their tendency to systematically overestimate the probability of occurrence of large values (i.e. to underestimate their return period), with ff i and n i,5 tending to be too frequently small (remind case G 1 of Figure 3).
The stability score SPAN 100 in Figure 8 shows that the most stable model is the lognormal case but this is because the lognormal distribution gives unreasonably huge estimates of large return values (as illustrated in Figure 9 for the station Antraigues, for example), giving very large normalization terms in the SPAN criteria (see ( 9)).The fact that the lognormal model has by far the worst reliability scores but the best stability score preaches for the conjoint use of these two family of scores not to misinterpret results.Stability of the Gamma and extended exponential distributions are very similar and fairly less good than the lognormal case.Then comes the Weibull distribution, and finally the Generalized Pareto distribution, which is clearly the least stable.Figure 9 illustrates the quality and spread of the fits depending on the distribution for the station Antraigues, when estimation is made on either subsample.Compared to Figure 7, it confirms that the Gamma and the extended exponential models perform almost likely.Median estimations differ by about 5% for the 100-year return level (303mm for the Gamma vs. 287mm for the extended exponential model) and by about 7% for the 1000-year return level (414mm vs. 386mm), with very similar widths of the 95% envelopes (e.g. ± 40mm for the 100-year return level).The lognormal model clearly fails to reproduce return periods larger than one year, giving much too heavy tails despite a resonably good fit of the bulk.Actually the skewness -which informs somehow on the "asymetry of the bulk"-is reasonably well estimated, whereas the kurtosis -which informs on the heaviness of the tail-is much overestimated.This is in line with Fig. 2 of Hanson and Vogel (2008), which shows that when the skewness of daily rainfall across the US is well estimated by the lognormal distribution, then the kurtosis is much overestimated.Note that Papalexiou et al. (2013) did not find such ill-fitted tails with the lognormal distribution but in their case fitting is made on the tail (i.e. on the largest values), whereas the lognormal model seems to fail when adjusting both the bulk and the tail of rainfall distribution.The Weibull and extended Generalized Pareto models give very similar fits up to the 50-year return period but the return level plot of the extended Generalized Pareto model is more convex than for the Weibull model, giving median estimation 8% larger for the 100-year return level (390mm vs. 358mm) and 35% larger for the 1000-year return level (799mm vs. 522mm).Width of the 95%-envelope is also larger both in absolute and relative values, in coherence with the SPAN 100 of Figure 8 at regional scale.Finally both the Weibull and extended Generalized Pareto models overestimate the return levels associated to 1-5 years, unlike the Gamma and extended exponential models.This tendency towards overestimation of the tail is actually a quite general feature observed for most the stations, giving to frequently low values of ff i and n i,T , as already stated.
The results of Figures 8 and 9 lead us to conclude that the best performance for the region is achieved by the Gamma and extended exponential models, which actually perform very similarly for Antraigues station.Note that exactly the same conclusions hold when focusing on the season-at-risk rather than considering the whole year, i.e. when computing the crossvalidation scores for the estimated seasonal distribution G s in (5) rather than for the year-round distribution G in (4.Due to its slightly better performance at regional scale for adjusting the tail of the distribution (FF and N 5 in Figure 8), we select the Gamma model (with two season and three WPs) for further consideration.

Selection of the mapping model
Figure 10 shows the 10 scores of evaluation of the mapping models of Table 3.The first comment is that, compared to Figure 8, the FF scores of leave-one-out cross-validation are for any mapping method of the same order as for the the local fits, while the SPAN scores are even slightly better.This means that i) no mapping method gives systematic over-or under-estimation of the very tail, and ii) mapping gives more stable estimations by smoothing out the sampling effect.However NMRSE scores are all larger, meaning that any mapping gives less accurate estimations of the full distributions than the local fits.Loss in accuracy is equivalent and relatively small for all kriging interpolations and the bivariate thin plate splines (with or without drift), while the trivariate thin plate spline and even more the linear model are less accurate.A closer look at the fits of all stations reveals that the strong loss in NRMSE for these two methods is actually due to few stations that are systematically very badly fitted -among which the station of Mayres of Figure 11-, which strongly deteriorates the spatial mean of the scores.Their less good performance is due to a lack of flexibility, which prevents them from adapting to local features.However at the same time, the lack of flexibility of these methods allow them for slightly increased stability in the tail, as shown by the SPAN scores in Figure 10.
Back to the kriging methods, the three tested alternatives give very similar fits, with slightly less stability when considering a drift in station altitude z, while considering the smoothed altitude Z is useless because a 1 in ( 12) is almost always zero.
The best kriging method for our region study in thus the simple kriging interpolation.This method is only slightly beaten in accuracy by the bivariate thin plate spline (with or without drift), but which is slightly less stable.However the TVD and KLD scores comparing the spatial stability of the mappings show that the bivariate thin plate splines are clearly more stable in space than all kriging methods.The linear model is even more stable but, as already said, it is much less accurate.Finally, comparing the five cases of thin plate spline shows that the three bivariate cases clearly outperform the trivariate case, both in terms of accuracy and stability.Comparing the bivariate case with drift (15) to the the trivariate case ( 16) shows the usefulness of considering non-linear weights of the distance (through the term h 2 i log(h i ) vs. h i ).Last but not least, whatever the method but particularly for the thin plate spline, better accuracy and stability is achieved when the smoothed altitude Z is considered rather than the station altitude z, as also found in Hutchinson (1998) for interpolating rainfall data.
We illustrate the results in Figure 11 for the Antraigues station, adding to that the case of the worst fit of the thin plate spline, which is for the station of Mayres.Mayres lies at about 500 m.a.s.l, as Antraigues, but it is located at the end of a funnel shaped valley surrounded by steep slopes (see Figure 1).This creates favorable conditions to the orographic intensification of rainfall, with the consequence that Mayres receives more rainfall than expected at this altitude, as also confirmed by Figure 2.For this reason, although the local fit of the Gamma model is reasonably good, the interpolated distributions underestimate the empirical values, even the small ones.This can be seen in Figure 11 by comparing the black curves, which were obtained independently of the data of Mayres, to the red curves of the kriging case, which are equal to the local fits since kriging is an exact interpolation.Although return levels are underestimated with all models, kriging and the bivariate thin plate spline manage however better the fit the data of Mayres in the leave-one-out framework, in coherence with the NMRSE values of Figure 10 at regional scale.For the station Antraigues, underestimation is also found for all methods due to smoothing but with much smaller extent than for Mayres.For both stations, comparing the red and black curves shows that kriging and the trivariate thin plate spline are too dependent on the data used for fitting since large differences are obtained whether the station is included or not in the estimation, in coherence with the TVD and KLD values of Figure 10    than the kriging method.The trivarite thin plate spline and the linear model should be avoided for our data due to their lack of flexibility.

Final regional model
Figure 12 illustrates the final regional models when both the Gamma parameters and the mapping models are estimated using all the available data.The map of the probability of daily rainfall to exceed 1mm is obtained from (3) with r = 1mm.The maps of the mean nonzero rainfall in the WPs of the season-at-risk (S2) are obtained as the product λ 2,k κ 2,k , k = 1, . . ., 3, with the notations of Table 1.The four maps of Figure 12 reveal the double effect of the Massif Central ridge, which both creates a climatological border and enhances orographic precipitation.The map of rainfall probability is conform to the climatology of the region (as shown by the colored points), with a smaller probability of rainfall in the Rhône valley and increased probability when approaching the relief due to orographic effect.Being more exposed to the west fluxes -which are the most common in showing a spatial correlation of rainfall twice as big as in WP2 and three times as big than in WP3.Remarkably, roughly the same factors are found when comparing the range of values of the means (resp. 5-36mm, 3-11mm, 2-10mm).There is thus a strong link between the spatial correlation of rainfall and the mean amounts.However the WPs do not only differ in the range of values of the mean amounts but also, and maybe even more, in the way these amounts are usually distributed over the region.This emphasizes once again the usefulness of considering subsampling over WPs in order to distinguish contrasted spatial pattern.The map of WP1 shows a strong intensification of rainfall along the Massif Central slope, while a clear decrease in the mean rainfall is found when passing the Massif Central ridge both towards the Massif Central plateau with means divided by three in 20km, and towards the Rhône valley with means divided by two in 20km.In WP2 the topography builds somehow a mask effect.The larger means are found along the Massif Central slope with a fast break when passing the Massif Central ridge.Daily means in the Massif Central plateau are half the values of the slope, while daily means in the Rhône valley are just slightly lower than in the slope.Finally the map of the mean nonzero rainfall in WP3 shows an inverse pattern to that of the probability of rainfall.The mean almost linearly decreases from the Rhône valley to the Massif Central plateau while the probability of rainfall almost linearly increases.The largest rainfall events in this WP are usually convective events of small extent occurring mainly in the Rhône valley, reason why the mean values are larger in this area, although the probability of rainfall is relatively low.
Last but not least, Figure 13 shows the map of the probability of daily rainfall to exceed 100mm.It reveals a clear concentration of higher probabilities of exceedance along the Massif Central ridge, with actually quite similar patterns as the averages of annual totals and annual maxima in Figure 2, with however much more pronounced disparities.It is up to 10 times less likely to exceed 100mm rainfall in the Rhône valley than along the ridge, and up to 1000 times less likely in the Massif Central plateau.Actually 100mm is expected to be exceeded several times a year along the ridge, about every year on the slope, and on average every 100 to 1000 year in the Massif Central plateau.

Conclusion and discussion
In this article we have presented an objective framework for estimating rainfall cumulative distribution function within a region when data are only available at rain gauges.For this we have proposed an objective procedure involving split sampling crossvalidation and using several evaluation scores allowing to validate the accuracy of both the bulk and tail of the distribution and the stability of these estimates when calibration data changes.We have applied this procedure to daily rainfall in the Ardèche catchment in southern France, a particularly challenging test case showing very strong disparities in rainfall in very short distance.Our results show that for this region, the best marginal model (among those tested) is a mixture of Gamma distributions over seasons and weather patterns, and that the best mapping model (among those tested) is the bivariate thin plate spline.
Although our procedure of selection is general and could be applied to any region of the world -and possibly to other variables (temperature etc)-, we stress that our conclusions are in themselves not universal.In particular, other marginal distributions may be more suitable than the Gamma in other regions of the world showing less or more heavy tails.Although the mixture of distributions over weather patterns has revealed efficient in other countries (e.g. in Norway and Canada, Brigode et al., 2014;Blanchet et al., 2015), it might be less relevant in, e.g., monsoon climate regions where the consideration of seasons seems essential and might be enough.We also stress that the goal of this article was neither to be exhaustive nor to find the best model for the region, but rather to present an objective cross-validation framework that can be applied to other variables and/or models of interest.Possible direction for improvement regards the choice of the spatial covariates to be used in the interpolation method.
There might be more relevant covariates than the geographical coordinates used in this study, e.g.considering atmospheric and terrain characteristics (Carreau et al., 2013;Kyriakidis et al., 2001).Finding good gridded covariates (and good regression models) is a subject of research in itself, and it is particularly tricky in areas with complex topography (Prudhomme and Reed, 1999;Weisse and Bois, 2001;Drogue et al., 2002;Beguería et al., 2009;Rogelis and Werner, 2013).The geographical distance itself might also be improved, e.g. by better accounting for the terrain characteristics (Gottardi et al., 2012;Evin et al., 2016) or by considering statistical distance (Ahrens, 2006).Also more robust estimates of the marginal parameters at station locations (i.e. of the θs) might be obtained by gathering observations of neighbor stations in order to increase the sample size, following the concept of regions-of-influence proposed by Burn (1990).Such idea has been quite widely used in the context of rainfall extremes (e.g.Carreau et al., 2013;Evin et al., 2016, for the studied region).However we anticipate the gain to be much less pregnant when interest is in modeling any rainfall -as in this study-, and not only the extreme ones since parameter estimation is already based on many data (several thousands).
Despite the above reservations of prudence, some other results seem to us to be generalizable.Among these is the fact that distributions showing sub-exponential tails (EGP for example) give usually unrealistic return levels when extrapolating far out the observations.Such distributions should not be used if return periods of several hundreds of years are of interest, as it is in many cases related to civil protection (dykes, dams etc, Paquet et al., 2013).Also the kriging method gives usually too patchy maps of rainfall hazard by sticking the observations, unless nugget effects are considered (which was not the case in this study).
Finally the linear model with spatial covariates usually fails to map rainfall hazard because it is highly unlikely to be ruled by simple enough physics for the parameters to be well represented as linear functions of the covariates, in particular in such complex topography (Carreau et al., 2013).
Last for not least, we put this study in a framework of temporal stationarity and we addressed the question of the spatial nonstationarity of the margins.Yet several studies showed temporal trend in the rainfall distribution in the region, particularly since the 80s and particularly along the Massif Central slope where daily rainfall is usually more intense (Blanchet et al., 2016b;Tramblay et al., 2011Tramblay et al., , 2013;;Vautard et al., 2015).Extending the proposed procedure to the case of nonstationary rainfall would be possible by considering the marginal parameters as parametric functions of time, e.g.linear models.This would increase the number of parameters but the split sample framework would still be valid.However the scores would have to be adapted to account for changing distributions.One way of doing this would be to transform the rainfall at time t to some variate independent of t.For example considering R t = exp{− exp(−G t (R t ))} would transform R t with CDF G t to a stationary Gumbel variate, to which the scores presented in this article could be applied for model evaluation and selection.A drawback however would be that the value of the scores would depend upon the chosen transformation.Also the SPAN score might have to be thought over because return levels in changing climates are not meaningful for quantifying risk (Katz, 2013).

Figure 1 .
Figure 1.Region of analysis.The blue polygon is the Ardèche catchement.The red points show the location of the stations.The upper triangle is the station Antraigues and the lower triangle the station Mayres (both lie at about 500 m.a.s.l.).The background shows the altitude in gray scale (1km raster cells).The top left insert shows a map of France with the studied region in red.The black lines are the 400 and 800 m.a.s.l.isolines.
k as the kth largest observed rainfall over C (1) i .Estimate r(2) i,T k is obtained numerically from Ĝ(2) i as described in Section 3.1.The normalization by the mean rainfall of C (1) i in (6) allows comparison of NRMSE over stations with different pluviometry.The smaller NRMSE
to a continuous uniform variate on (0, 1) is proposed inRenard et al. (2013) and extensively described inBlanchet et al. (2015).For i ranging over the set of Q stations, we thus obtain a sample of Q uniform variates.Scores are calculated as for FF by comparing the empirical densities of N (12) T to the theoretical uniform density, giving the scores AREA(N (12) T y) where u is an unknown smooth function and (x, y) are uncorrelated errors with zero means and equal variances.The function u is estimated by minimizing

i
) are then obtained by computing the mean of the Q values.MEAN(TVD

Figure 4 .
Figure 4. Right: Monthly percentage of occurrence of the three WPs.Left: Boxplot of the monthly averages of daily nonzero rainfall.Each boxplot contains 42 points (one point per station).

Figure 5 .
Figure5.Schematic summary of the full cross-validation procedure for selecting both the marginal and mapping models.

Figure 6 .
Figure 6.Scores of cross-validation when G s,k are Gamma distributions and the number of seasons and WT varies: S ∈ {1, 2} and K ∈ {1, 3}.The values of (S, K) are indicated in the x-labels.Each boxplot contains 100 points.

Figure 7 .
Figure 7. Case of Antraigues when G s,k are Gamma distributions and the number of seasons and WT varies: S ∈ {1, 2} and K ∈ {1, 3}.

Figure 8 .
Figure 8. Scores of cross-validation when G s,k is either the extended exponential (eexp), extended Generalized Pareto (egp), Gamma (gamma), lognormal (lnorm) or Weibull (wei) distribution, with S = 2 and K = 3.Each boxplot contains 100 points.The boxplots of reliability scores in the lognormal case are missing because they lie far above the upper range of depicted values.

Figure 9 .
Figure 9. Case of Antraigues when G s,k is either the extended exponential (eexp), extended Generalized Pareto (egp), lognormal (lnorm) or Weibull (wei) distribution, with S = 2 and K = 3.The dotted lines show the 95%-envelope of return level estimates over the 100 subsamples.The plain line shows the median estimates.The gray points show the full sample (35 years).Each estimation is based on half these points.Case of the Gamma distribution is shown in the right panel of Figure 7.
at regional scale.Finally, comparing the envelope widths in red and black in Figure 11 confirms that interpolation increases stability of the estimates, as also revealed by the SPAN score of Figure 11.

Figure 10 .Figure 11 .
Figure 10.Scores of mapping when G s,k are Gamma distributions with S = 2 and K = 3 whose parameters are interpolated with the mapping models of Table 3.The two first rows show leave-one-out cross-validation scores.Each boxplot contains 200 points.The third row compares interpolations at a given station whether the data of this station are used or not in the interpolation.Each boxplot contains 100 points.

Figure 13 .
Figure 13.Map of the probability of daily rainfall to exceed 100mm.The points show the location of the stations.

Table 2 .
Summary of the considered scores for evaluating marginal and mapping models.

Table 3 .
Mapping models considered in this study, with involved coordinates.Kriging method provides exact interpolation, unlike the linear regression and thin plate spline.
Table 1 to the two subsamples, getting estimates Ĝ(1) i of each distribution and for every station.
We conclude following the results of Figures 10 and 11 that the best interpolation method (among those tested) is the bivariate thin plate spline with drift in smoothed altitude, which is slightly more accurate but much more stable spatially stable