M Rousseau
email: m.rousseau@brgm.fr

O Cerdan
email: o.cerdan@brgm.fr

A Ern
email: ern@cermics.enpc.fr

O Le Maître

P Sochala
email: p.sochala@brgm.fr

Study of overland ow with uncertain inltration using stochastic tools

The saturated hydraulic conductivity is one of the key parameters in the modelling of overland ow water uxes. In this study, this parameter is dened as a stochastic parameter, idealized as a piecewise constant random eld with uniform distribution. This paper aims at investigating the eects of the spatial and temporal scales in uncertainty propagation within overland ow models, and at identifying the localization of the most inuential saturated hydraulic conductivity using sensitivity analysis. The results show that the inuence of saturated hydraulic conductivity depends on the soil saturation and its spatial localization.

For instance, in case of low saturated soils, the most inuent parameter is the one located downslope, whereas in case of high saturated soils, the most inuent one is either the most inltrating or the intermediate one. The results indicate where eorts should be concentrate when collecting input parameters to reduce modelling uncertainties.

Introduction

.

Indeed, the performance of models directly depends on the validity of the input parameters.

.

. Dierent methodologies have been elaborated to measure directly saturated hydraulic

conductivities. However, the obtained values for K s depend on the chosen methodology and most importantly, the spatial representativity of these measurements remains rather limited [START_REF] Wösten | Pedotransfer functions: bridging the gap between available basic soil data and missing soil hydraulic characteristics[END_REF]. In most model applications, values for K s are estimated through the application of pedo-transfer functions (PTF) using basic soil properties [START_REF] Bouma | Using soil survey data for quantitative land evaluation[END_REF][START_REF] Bouma | Transfer functions and threshold values: from soil characteristics to land qualities[END_REF][START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF]. By testing dierent PTF's to predict K s , it was concluded [START_REF] Tietje | Accuracy of the saturated hydraulic conductivity prediction by pedo-transfer functions compared to the variability within fao textural classes[END_REF] that predicting K s using a PTF is not always accurate owing to the inherent variability of K s . Furthermore, using a set of data to compare dierent measurement and estimation methods, a high variability of K s (more than 79%) has been observed [START_REF] Baroni | Uncertainty in the determination of soil hydraulic parameters and its inuence on the performance of two hydrological models of dierent comple[END_REF]. To overcome this lack of accuracy, a possible approach consists in calibrating parameters, but the resulting values are often valid only for the used conguration and moderate variations thereof.

An alternative approach already suggested in [START_REF] Tietje | Stochastic modeling of the unsaturated water ow using auto-correlated spatially variable hydraulic parameters[END_REF] consists in considering K s as a stochastic parameter instead of being estimated by deterministic approaches. It is today well admitted that probabilistic modelling provides ecient means to quantify parameter uncertainty.

Uncertainty Analysis (UA) considers the uncertain parameters of a model as random objects, and the objective is to compute or characterize the induced variability in the model solution or in quantities of interest. For highly uncertain data within non-linear models, as in hydrology applications, so-called global UA methods, which study the eects of all the input parameters simultaneously, are needed. One essential step in UA is the denition of a random model for the uncertain parameters. Random models with dierent levels of complexity can be considered. For instance, a relatively simple approach is the Generalized Likelihood Uncertainty Estimation (GLUE) procedure [START_REF] Beven | The future of distributed models: Model calibration and uncertainty prediction[END_REF] which is a Monte Carlo (MC) method generating a high number of parameter sets to compare the predicted model responses with observed responses and to accept or not some simulations through some chosen likelihood measure. Being a Bayesian approach, this likelihood measure can be updated for each new set of observed responses. Numerous studies are based on a Bayesian framework [START_REF] Huard | Calibration of hydrological model gr2m using bayesian uncertainty analysis[END_REF][START_REF] Kuczera | Monte carlo assessment of parameter uncertainty in conceptual catchment models: the metropolis algorithm[END_REF][START_REF] Laloy | Modelling intercrop management impact on runo and erosion in a continuous maize cropping system: Part i. model description, global sensitivity analysis and bayesian estimation of parameter identiability[END_REF][START_REF] Thiemann | Bayesian recursive parameter estimation for hydrologic models[END_REF]. Bayesian statistics mean that input parameters are considered as probabilistic variables having a joint posterior probability density function (pdf ). Dierent methods exist for sampling posterior pdf 's. The Monte Carlo Markov Chain (MCMC) sampler is often used in hydrology models, the earliest general (and most popular) method being the MetropolisHastings algorithm [START_REF] Hastings | Monte carlo samping methods using markov chains and their applications[END_REF][START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF]. In hydrology, various recent studies have aimed at improving MCMC samplers: the Shued Complex Evolution (SCE) Metropolis algorithm (SCEM) [START_REF] Vrugt | A shued complex evolution metropolis algorithm for optimization and uncertainty assessment of hydrologic model parameters[END_REF], which is a modied version of the SCE global optimization algorithm [START_REF] Duan | Eective and ecient global optimization for conceptual rainfall-runo models[END_REF], or the Dierential Evolution Adaptive Metropolis algorithm (DREAM) [START_REF] Vrugt | Treatment of input uncertainty in hydrologic modeling: Doing hydrology backward with markov chain monte carlo simulation[END_REF]. In contrast with Bayesian statistics, an alternative approach providing a complete probabilistic description considers the unknown parameters as random variables described by a xed pdf. This approach, which is more adapted to problems where little data is available, is undertaken in the present work. Once the probabilistic framework is settled, it remains to characterize the model output variability in terms of input variability. To this purpose, we rely on MC methods since such methods provide an eective and robust methodology to generate a sample set of model solutions by sampling input parameters. We mention another recent methodology based on stochastic spectral methods [START_REF] Ghanem | Stochastic nite elements: a spectral approach[END_REF]. The advantage is that a more complete probabilistic description of model output is achieved, but the methodology becomes computationally demanding when the input parameters are described by a large number of random variables.

In the present paper, we focus on the impact of the variability in the saturated hydraulic conductivity K s on overland ows produced by various rainfall events in various geometries.

A general probabilistic description of the saturated hydraulic conductivity is to model it as a random eld. Although very rich, this type of model needs a substantial amount of information for its description, and is, therefore, not well adapted to the present setting. Moreover, extracting simple information in view of practical hydrology purposes from complex probabilistic models is an intricate issue. For these reasons, we rely on simpler probabilistic models where the saturated hydraulic conductivity is idealized as a piecewise constant random eld over distinct portions of the simulation domain which are specied a priori. Each spatial subdomain is characterized by a unique Random Variable (RV) yielding the corresponding saturated hydraulic conductivity in the subdomain. Specically, the random conductivity eld K s (x, θ), where θ is a random event, is written as

K s (x, θ) = i 1 Ωi (x)K s,i (θ), (1)
where 1 Ωi is the indicator function of the i-th spatial subdomain and {K s,i (θ)} is a set of (positive) real-valued random variables. In addition, the saturated hydraulic conductivity within a subdomain is assumed to be statistically independent from the others. As a result, the uncertain hydraulic conductivity eld is modelled using a nite set of independent random variables, whose cardinality is equal to the number of subdomains considered in the simulation, and the joint probability density function (pdf ) of the random variables K s,i factorizes in the form

p({K s,i }) = i p i (K s,i ). (2)
This idealization of the actual conductivity eld is motivated by the physical reality. If one thinks of subdomains as homogeneous agricultural plots (or parcels of vegetation), the variability of K s inside the parcel is usually negligible compared to the variability from a parcel to another because of the homogenization created by agricultural practices. Moreover, the present model can be subsequently rened by introducing inner variability within the parcels if additional information on soil properties within parcels is available.

The objective of this work is twofold. Firstly we consider test cases involving dierent spatial and temporal scales so as to investigate the eect of these scales in uncertainty propagation within overland ow models. Indeed, one salient result discussed hereafter is that the most inuential input parameter on model output variability depends on the spatial and temporal scales of the processes of interest. Secondly, with an eye toward practical hydrology issues, we use sensitivity analysis tools as a means to identify the most inuential input parameters. This information can be, for instance, valuable to decide on where to concentrate additional measurement eorts to improve eld knowledge. Moreover, within a given test case, we consider various possibilities for the spatial organization of the parcels so as to study the eect of this organization on model output variability. One interesting result in this direction, relevant to landscape management issues, is the comparison presented hereafter between three grass strips distributed evenly within a eld and a single grass strip located near the eld outlet. Finally, we observe that the present work is not dedicated to a specic eld study, but the purpose is instead rather general and aims at assessing the possibilities oered by the present methodology to improve our understanding of variability in overland ows due to inltration. For this reason, we focus on a simple two-dimensional setting where the ow is described by a one-dimensional shallow-water ow model including friction and inltration, the latter being described by the GreenAmpt model. More elaborate ow models can be considered. We also mention a dierent approach [START_REF] Wang | Probabilistic predictions of inltration into heterogeneous media with uncertain hydraulic parameters[END_REF][START_REF] Wang | Reduced complexity models for probabilistic forecasting of inltration rates[END_REF] to compute pdf 's of inltration rates and inltration depths. This paper is organized as follows. Section 2 briey describes the rainfall-runo model and the numerical method used in the deterministic overland ow simulations once values for the random input parameters are specied. Section 3 introduces the stochastic approach and the statistic tools used to propagate and analyze the uncertainties in model output.

Section 4 presents the two test cases designed to evaluate the impact of uncertainties in K s and of the spatial localization of these uncertainties on overland ow. Results are discussed in Section 5. Finally, conclusions are drawn in Section 6.

The setting

In this section, we present the physical model and its numerical resolution.

Physical model

We are concerned with overland ows where the water depth is much smaller than the characteristic horizontal size of the eld of study (see Figure 1). Such ows can be de- scribed by the 2D shallow water (SW) equations which are obtained from the 3D incompressible free-surface NavierStokes equations by averaging on the vertical direction under some simplifying assumptions, in particular hydrostatic pressure and negligible vertical velocity [START_REF] Gerbeau | Derivation of Viscous Saint-Venant System for Laminar Shallow Water; Numerical Validation[END_REF][START_REF] Stoker | Water Waves: The Mathematical Theory with Applications[END_REF][START_REF] Viollet | Mécanique des uides appliquée[END_REF]. Neglecting also the ow transverse to the main slope direction, we obtain the 1D SW equations which express mass and momentum conservation as follows:

z x R(x,t) u(x) I(x) h(x) z(x)
∂ t h + ∂ x (hu) = R -I, (3)
∂ t (hu) + ∂ x hu 2 + 1 2 gh 2 = -gh (∂ x z + S f ) , (4)
where h [L] is the water depth, u [L/T] the depth-averaged velocity, z [L] the ground surface elevation, and g [L/T 2 ] the gravitational constant (where L and T denote length and time units, respectively). The source term R -I [L/T] corresponds to the dierence between the rainfall rate R and the inltration rate I. The quantity S f [L/L] accounts for friction eects. The value of S f depends on the properties of the soil surface and can be estimated from calibration or published values. DarcyWeisbach's formula is often used [START_REF] Chow | Open Channel Hydraulics[END_REF][START_REF] Esteves | Overland ow and inltration modelling for small plots during unsteady rain: numerical results versus observed values[END_REF][START_REF] Fiedler | A numerical method for simulating discontinuous shallow ow over an inltrating surface[END_REF][START_REF] Viollet | Mécanique des uides appliquée[END_REF]:

S f = f |u| 8gh u, ( 5
)
where f is the possibly time and space-dependent DarcyWeisbach's roughness coecient and |u| the module of the velocity u.

We use the GreenAmpt model [START_REF] Green | Studies on soil physics: 1, ow of air and water through soils[END_REF] to predict cumulative inltration through dry or wet soils. Herein, we consider the formulation developed in [START_REF] Mein | Modeling inltration during a steady rain[END_REF] for rainfall-runo predictive models. This formulation postulates, at any point x in space, a well-dened wetting front propagating vertically and separating a fully saturated zone from a zone at the initial soil moisture. At any point on the surface, the inltration capacity I c (t) [L/T] at time t is calculated as follows (the dependency on the space variable is omitted for simplicity):

I c (t) = K s 1 + (θ s -θ i ) h f + h(t) I * (t) , (6)
where I * (t)

Numerical resolution

A well-balanced nite volume method is used to discretize the SW equations, which we rewrite in the general form ∂ t U + ∂ x F (U ) = S(U ), where U is the vector of conservative variables, F the ux vector, and S the source term. Specically,

U = h hu , F (U ) = hu hu 2 + 1 2 gh 2 , S(U ) = R -I -gh(∂ x z + S f )
.

The domain is divided into cells (indexed by i) of the form C i = [x i-1/2 , x i+1/2 ] and of length ∆x > 0 taken constant for simplicity. The GreenAmpt model is applied locally in each mesh cell. To obtain a second-order scheme, the variables need to be reconstructed at cell interfaces. We denote by U i+1/2 ± the conservative variables computed at either side of the interface x i+1/2 using an ENO-type (Essentially Non Oscillatory) reconstruction [START_REF] Harten | Uniformly high-order accurate nonoscillatory schemes[END_REF]. Moreover, the ground surface elevation z is described as a piecewise constant function, and cell-interface values z i+1/2 ± are also reconstructed. This yields the following scheme written here in space semi-discrete form:

∆x d dt U i (t) + F i+1/2 -F i-1/2 = ∆x R i -I i -gh i S f,i + ∆x 0 S s,i ,
where the subscript i refers to the value in the cell C i and the uxes F i±1/2 are computed using the HLL ux (see e.g., [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws and Well-Balanced Schemes for Sources[END_REF]). The source term S s,i accounts for the term -gh∂ x z in the source term S. To evaluate the uxes F i±1/2 and the source term S s,i , an hydrostatic reconstruction scheme is applied, as described in [1, 2, 7, 27]. Specically, we set

           h i+1/2 L = max 0, h i+1/2 -+ z i+1/2 --max(z i+1/2 -, z i+1/2 + ) , h i-1/2 R = max 0, h i-1/2 + + z i-1/2 + -max(z i-1/2 -, z i-1/2 + ) , U i+1/2 L = h i+1/2 L , h i+1/2 L u i+1/2 - t , U i-1/2 R = h i-1/2 R , h i-1/2 R u i-1/2 + t .
Then, the HLL ux is evaluated using (U i+1/2 L , U i+1/2 R ), and the source term S s,i is evaluated as

S s,i = 1 ∆x g 2 h 2 i+1/2 L -h 2 i+1/2 -+ h 2 i-1/2 + -h 2 i-1/2 R + h i-1/2 + -h i+1/2 -z i+1/2 --z i-1/2 + .
Finally, for time discretization, we use a second-order explicit RungeKutta method based on the Heun scheme, except for the friction term S f which is treated semi-implicitly at each stage of the Heun scheme [START_REF] Bristeau | Boundary Conditions for the Shallow Water Equations solved by Kinetic Schemes[END_REF]. This leads to a second-order accurate overland ow model with inltration that we now use to study uncertainty propagation.

Stochastic model and statistic tools

In this section, we describe the stochastic model and the statistic tools used to analyze the results.

Stochastic model

We are interested in uncertainty propagation stemming from the uncertain parameter K s in the SW equations. As discussed in the introduction, we consider a simple approach where we subdivide the physical domain into p (with typically p = 2 or 3 in our numerical results) subdomains and assign to each subdomain a single uncertain parameter K s which is a RV with known pdf. The uncertain input parameters of the model are then regrouped into a single vector-valued RV such that X = (X 1 , . . . , X p ) = (K s,1 , . . . , K s,p ). We assume that the RV's X i are independent, but can have dierent pdf 's denoted by p Xi . We consider a uniform distribution for all the subdomains because of the relatively low range of values taken by K s within each subdomain (however, high contrasts are considered between subdomains).

In the present setting, the pdf p Xi depends on its corresponding subdomain i only through the minimal and maximal bounds on K s,i . Moreover, since the RV's X i are assumed to be independent, the pdf of X factorizes into the form

p X (x 1 , • • • , x p ) = p i=1 p Xi (x i ). (7)
The uncertain output quantities of the model are the peak runo rate and the runo coecient for a specic rainfall event. Let Y denote any of these output quantities. Once a realization of X, say x, is known, a realization of Y , say y(x), is obtained by solving numerically the corresponding deterministic problem described in Section 2.

Moments and pdf

Assuming that Y has nite second-moment, the expectation and the variance of Y are dened as

E[Y ] = y(x)p X (x)dx and V (Y ) = y(x) -E[Y ] 2 p X (x)dx, so that V (Y ) = E[Y 2 ] -E[Y ] 2 .
We are interested in evaluating various statistical quantities related to the model output Y . To this purpose, we use Monte Carlo (MC) simulations. Let X = {x (1) , . . . , x (M ) } be a sample set of the input stochastic parameters, where M is the sample set dimension and

x (m) , 1 ≤ m ≤ M , are realizations of X. Let Y = {y (1) , . . . , y (M ) } be the corresponding sample set of the model output such that, for each 1 ≤ m ≤ M , y (m) = y(x (m)
) is the model response to the vector of input parameters x (m) . The empirical estimators for the expectation and the variance are

Ê[Y ] = 1 M M m=1 y (m) and V (Y ) = 1 M M m=1 y (m) -Ê[Y ] 2 . ( 8
)
To estimate the pdf of a random variable, we use the kernel density estimator, also called ParzenRozenblatt method, which is a generalization of the histogram method. The pdf of Y is estimated as

pη (y) = 1 M η M m=1 G y -y (m) η ,
where G is a specic pdf used as kernel and η is a smoothing parameter. The most commonly used kernel is the Gaussian function G(x) = (2π) -1/2 exp (-1 2 x 2 ). Thus, the pdf at a point y is estimated by the number of observations close to y and counterbalanced by the distance of these observations to y. The kernel distribution function allocates more important weights to observations near the point y and weaker weights to distant observations. The parameter η xes the kernel function width and, therefore, controls the smoothness of the estimated pdf pη . The smaller the parameter, the more accurate the estimation of the pdf; however, too small values for η can generate spurious data artifacts if the sample set is not suciently large. An illustration is presented in Section 4.2.

Sensitivity analysis

Sensitivity Analysis (SA) allows one to assess the relative contribution of each uncertain input parameter to model output variability and, in particular, to identify key parameters by establishing a hierarchy within the input parameters according to their inuence on the output variability. Variance-based global SA methods based on Sobol indices [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF] determine which part of the response variance results from the variance of each input or group of inputs. The sensitivity of the response Y to the input parameter X i can be quantied by the rst-order sensitivity index S i dened as

S i = V i V (Y ) , V i = E E[Y |X i ] 2 -E[Y ] 2 ,
where E [Y |X i ] is the conditional expectation of Y given the value of X i (see [START_REF] Bouma | Transfer functions and threshold values: from soil characteristics to land qualities[END_REF] below for its denition). More generally, higher-order sensitivity indices quantify the sensitivity of the model response to interactions among input parameters. Let i denote a non-empty subset of indices such that i ⊆ {1, . . . , p} and let ∼i = {1, . . . , p} \ i. The sensitivity index S i is dened as

S i = V i V (Y ) , V i = V (E [Y |X i ]) - ∅ =j i V j , where V (E [Y |X i ]
) is the variance of the conditional expectation of Y given the value of X i . This conditional expectation is dened as

E [Y |X i ] = y(x)p X ∼i (x ∼i )dx ∼i , (9)
where p X ∼i and dx ∼i are, respectively, the density and the probability measure of x ∼i (con-

ventionally, E [Y |X i ] = Y if i = {1, . . . , p} and ∼i is empty). Observing that E [E[Y |X i ]] = E[Y ], we obtain V (E[Y |X i ]) = E E[Y |X i ] 2 -E[Y ] 2 .
Furthermore, the law of total variance states that

∅ =i⊆{1,...,p} V i = V (Y ), so that ∅ =i⊆{1,...,p} S i = 1.
Following Homma and Saltelli [START_REF] Homma | Importance measures in global sensitivity analysis of nonlinear models[END_REF], it is convenient to consider for a single index i ∈ {1, . . . , p}, the total sensitivity index S T,i which evaluates the total sensitivity of the model response Y to the input parameter X i , including X i alone and all interactions with the other input parameters X j , j = i. Computing this index instead of the high-order sensitivity indices allows one to reduce computational costs by avoiding tedious calculations. The total sensitivity index S T,i is evaluated as follows:

S T,i = 1 - V ∼i V (Y ) , V ∼i = E E[Y |X ∼i ] 2 -E[Y ] 2 ,
where V ∼i is the variance of the conditional expectation of Y given all the parameters except X i . The interpretation of the indices S i and S T,i is the following: X i is an inuential parameter if S i is important, whereas X i is not an inuential parameter if S T,i is small. Moreover, S i close to S T,i means that interactions between X i and the other parameters are negligible.

MC simulations are used to estimate the rst-order sensitivity indices S i and the total sensitivity indices S T,i . To this purpose, we estimate the quantities V i and V ∼i . To save computational costs when evaluating these variances [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]

, the expectations E E[Y |X i ] 2
and E E[Y |X ∼i ] 2 are computed as a unique integral by making use of two independent M-samples of input variables, X and X * , in such a way that

Vi = 1 M M m=1 Y x (m) × Y x (m) i , x * (m) ∼i -Ê[Y ] 2 , V∼i = 1 M M m=1 Y x (m) × Y x * (m) i , x (m) ∼i -Ê[Y ] 2 .
Finally, the rst-order sensitivity index S i and the total sensitivity index S T,i are estimated

as Ŝi = Vi V (Y ) , ŜT,i = 1 - V∼i V (Y ) .
In practice, the computational procedure requires two samples of input parameters, each of dimension M , and M × (2p + 1) deterministic model evaluations to calculate all the rst-order and total sensitivity indices.

Test cases

The section presents the test cases and a brief performance evaluation of the methodology.

Presentation

To evaluate how uncertainties on the values of K s and its spatial localization can impact the surface runo during various types of rainfall events, we focus on two output quantities:

the peak runo rate at the outlet, Q max , and the runo coecient, RC, dened as the total volume of runo divided by the total volume of rainfall.

A one-dimensional at slope of length L with an inclination gradient ∂ x z = 2% is considered. Uniform friction coecient and inltration parameters (except K s ) are chosen with values The values of the rainfall intensity R, the rainfall duration T R , and the simulation time T are specied in Tables 1 and3 for the two spatial congurations described below.

f = 0.25, θ s -θ i = 0.3, h f = 0.023. A constant rainfall intensity R(t) [L/T] is
For the rst spatial conguration, named Three-eld, we consider a 4.8 m long domain divided into three subdomains, each one with its own saturated hydraulic conductivity K s,i , i ∈ {1, 2, 3}, which is a RV independent of K s,j , j = i. Each K s,i has a uniform distribution K s,i ∼ U[K min s,i , K max s,i ], where K min s,i and K max s,i are the minimal and maximal values which can be taken by K s,i in its corresponding subdomain. To assign these minimal and maximal values, we consider three choices, each representing realistic values for a given soil type. We refer to these choices using an index 2(a). Figure 2(b) presents the hydrographs for the case where K s,+ , K s,o , and K s,-are all equal to their respective mean values (Table 2), and the impermeable conguration. The signicant dierences observed emphasize the importance of inltration processes.

For the second spatial conguration, named Grass strip(s), we consider a 318 m long domain in which grass strips (GS) of total width equal to 6 m are added. Two spatial localizations of the GS are considered, as represented in Figure 3: either three narrow, 2 m wide GS are equally spaced or one large, 6 m wide GS is located at the outlet. We assign a saturated hydraulic conductivity to the GS and another one to the remaining part of the soil surface called the eld. Both K s are independent RV's with uniform distribution. The values taken by K s on the GS are higher than those on the eld (see Table 4).

Performance evaluation

Before discussing our results in the next session, we verify the numerical procedure on the Three-eld test case with SRE and the spatial localization + -o for K s . Figure 4(a) presents the convergence of the rst-order and total sensitivity indices as a function of the sample set dimension M . For both order indices (rst and total), the convergence is obtained quickly around M = 1000. 3 narrow GS 1 l arge GS 2% 2%

L 2% R(t) K s,- K s,o K s,+ (
f i el d f i el d R( t ) R( t )

GS conf i gurat i on 1 GS conf i gurat i on f i el d f i el d

Results and discussion

This study focuses on the eects of the uncertainties in K s and of the spatial localization of the soil types on the peak runo rate at the outlet, Q max , and the runo coecient, RC.

The two test cases Three-eld and Grass-strip(s), detailed in Section 4, are investigated.

In particular, the dierent statistic estimators (moments, pdf 's, and sensitivity indices) are calculated for the two rainfall events SRE and LRE, with a sample dimension equal to 100,000 and a uniform distribution for all RV's.

Three-eld test case

Figures 5 and6 present the 100,000 couples (Q max ,RC ) for the six possible spatial localizations of soil types and for SRE and LRE, respectively. An important point is that there is an important correlation between the two outputs Q max and RC for each spatial localization. Indeed, all the correlation coecients are greater than 0.9. Concerning SRE (Figure 5), the simulations even tend to line up in a curve. We observe that Q max = 0 when K s,+ is downslope due to a complete inltration of the rain and the upslope runo.

The cloud of points for the spatial localization + -o contains the one for -+ o which corresponds to the weakest discharges. The congurations where K s,-is located downslope are similar whatever the positions of K s,+ and K s,o because the values of K s,+ and K s,o are suciently important to inltrate all the rain. Therefore, the clouds of points for SRE depend on the K s located downslope. Concerning LRE (Figure 6), Q max and RC increase since the rainfall duration is longer. As a result, the inuence of K s,+ is more pronounced for LRE and contributes more to the discharge at the outlet, whereas the inuence of K s,- decreases. Therefore, the clouds of points for LRE depend essentially on the position of K s,+ within the domain. Figures 5 and6 stress the importance of the spatial distribution of the soil types since the outputs are inuenced by the K s located downslope for SRE and by the localization of K s,+ for LRE. To better understand why this dierence is observed by changing the rainfall duration, we focus on the inltration process over the domain. Equation [START_REF] Boegh | Incorporating remote sensing data in physically based distributed agro-hydrological modelling[END_REF] implies that the ratio I c /K s tends to 1 when the inltrated water volume tends to innity (corresponding to a saturated soil). To study the eect of increasing the rainfall duration on soil saturation, Figure 7 presents the condence interval (i.e. minimal and maximal values for the 100,000 model responses) of the ratio I c /K s at nal time, as a function of spatial position, in grey for SRE and in black for LRE. As expected, the soil is more saturated for LRE and the ratio is closer to 1. Besides, the variability of the condence interval is in general the highest for the subdomain where K s,-is located and the weakest for the subdomain where K s,+ is located. The variability decreases as a function of soil saturation since the more saturated the soil, the smaller the variability, except for some limit cases where there is no runo on the concerned subdomain.

Table 5 presents the mean µ, the standard deviation σ, the coecient of variation cov = σ/µ, the median P 50 , and the 90th percentile P 90 related to Q max (white lines) and RC (grey lines). On the whole, there is more dispersion on the estimated values for SRE. For instance, cov varies between 10% when K s,-is located downslope and reaches 42% (for + -o) and 217% (for -+ o) when K s,o is located downslope. Besides, the values are higher for SRE than for LRE. The increase of the rainfall duration leads to a decrease in the dispersion values, therefore dispersion depends on the state of soil saturation (as the variability of the ratio I c /K s observed previously in Figure 7). Furthermore, for SRE, the distribution is not uniform when K s,o is downslope since the median and the mean are dierent. Finally, similar conclusions can be drawn from the statistical values associated with RC. Figures 8(a) and 8(b) present the pdf 's of Q max estimated by the ParzenRozenblatt method and plotted by groups in function of the K s which inuences the most the discharge at the outlet. A rst important point is that the pdf shape shows that the process studied is not Gaussian. Concerning SRE, the spatial localizations of soil types yielding clouds of points that are correlated and uniformly distributed in Figure 5 (the ones where K s,-is downslope), generate a spread pdf looking like a rectangular function. The pdf resulting from the conguration -+ o has a marked peak owing to the numerous null discharges observed. This marked peak does not have the expected form on the left part because it is dicult to approximate accurately such a pdf (resembling a Dirac function) by a Gaussian kernel. Concerning LRE, the six curves on Figure 8(b) are very close and have the form of a at bell on top and are almost symmetrical with respect to zero (i.e, with respect to the mean because the samples are standardized). So, contrary to SRE, the spatial distribution of K s does not inuence the distribution of Q max for LRE. The estimated pdf 's for RC lead to the same conclusions. We can conclude that the dispersions calculated in Table 5 for SRE are conrmed by the non-uniform distribution obtained in the pdf curves. the highest index corresponds to the parameter K s located downslope, thus corroborating the previous conclusions on the most inuent K s . For instance, in Figure 5, for the spatial localizations where K s,-is located downslope, the clouds of points are similar. Switching K s,o and K s,+ does not impact the outlet discharge, meaning that only K s,-inuences this quantity, and indeed the indices S -are equal to 1. Concerning LRE, since the rainfall duration is longer, more runo is generated in the most inltrating subdomain because of the decreasing of the inltration capacity. In Figure 9(b), the most inuent parameter is either K s,+ (with S + ≈ 72%) or K s,o (with S o ≈ 70%) when K s,+ is located upslope. Moreover, K s,-is not very inuent, and contrary to Figure 5, the three parameters K s are not negligible in the sensitivity analysis. Moreover, the total sensitivity indices are equal to the corresponding rst-order index, that is, S T,i ≈ S i . These equalities mean that there is no signicant interaction between the input parameters. Concerning the runo coecient RC, the sensitivity analysis leads to the same conclusions. In practice, in case of soils with low levels of saturation (for SRE), it is important to focus the measurements on the parameter closer to the outlet. For more saturated soils (e.g., for LRE), the measurements should focus on the most inltrating parts of the domain.

To study the eect of the length L of the domain, we have also tested the case where L = 48 m with LRE. It is interesting to notice that the length of the domain does not aect the results. The clouds of points and the pdf 's have the same shape, and the most inuent sensitivity index is the same, i.e, K s,o in cases where K s,+ is located upslope or K s,+ in other cases. The only signicant dierence is that the sensitivity index for K s,+ vanishes when K s,+ is located upslope. This result can be explained by the fact that longer domains lead to an augmentation of the distance between the upslope subdomain and the outlet.

Grass strip(s) test case

Figure 10 presents the couples (Q max ,RC ) corresponding to the 100,000 model responses for the four congurations (1 or 3 GS; SRE or LRE). In each conguration, the clouds of points are well correlated and, as previously, the values of Q max and RC are higher for LRE than for SRE. Concerning Q max , for both SRE and LRE, the values are contained approximately in the same intervals whatever the spatial localization of the GS. Concerning RC, its values are slightly higher for the spatial conguration with 1 GS, and this eect is more signicant for SRE. We conclude that the spatial localization of the GS has very little inuence on the variability of the runo, and almost none on that of the ow at the outlet.

Figure 11 presents the condence interval (for the 100,000 model responses) of the ratio I c /K s at nal time, as a function of spatial position. As expected, because of the duration of the rainfall events, the ratios are closer to 1 for LRE. Besides, compared to Figure 7, the values taken by the ratio are very close to 1, meaning that the soil is almost saturated.

For each rainfall event, the values taken by I c /K s for the two spatial congurations (1 GS and 3 GS) are very close. Furthermore, we observe that for SRE, the variability of the ratio I c /K s is between 2 and 4 times higher for the GS than for the eld. Conversely, for LRE, the variability is approximately 3 times more important for the eld than for the GS. An interesting result concerning SRE is that the eect of having 3 GS instead of 1 GS downslope, is to somehow homogenize the level of variability of I c /K s along the eld.

Statistical values (µ, σ, cov, P 50 , and P 90 ), not presented here, conrm that the spatial congurations with 1 GS and 3 GS are similar regarding Q max for LRE, and very close for SRE. Concerning LRE, in agreement with the almost essentially at shape of the pdf 's (Figure 12), we report the same values for the model outputs with the mean values of the parameters, the mean estimation, and the median. Concerning SRE, highly marked peaks are observed with signicantly dierent values (4.2 for 1 GS versus 3.1 for 3 GS). These peaks explain the dierence between the mean and the median. Moreover, the mean values of the model outputs dier from the model outputs with the mean parameters. This underlines the importance of non-linear processes. The statistical values and the estimated pdf 's for RC lead to the same conclusions.

Figures 13 and14 present the three statistic estimators µ, P 50 , and P 90 for the peak runo rate max t Q(x, t) as a function of spatial position, and the two deterministic values of this quantity (taking K s = K min s and K max s

). The curves for P 50 almost coincide with those for max t Q(x, t) calculated with the value K s = K mean s . Contrary to LRE where equality is obtained, the median is inferior to the mean for SRE. Both for the 1 GS and 3 GS congurations, the distribution is not uniform in space. Moreover, for both SRE and LRE, RC is slightly higher with the 3 GS conguration. Although the runo volumes are comparable for 1 GS and 3 GS, the spatial distribution of maximal discharges varies. Indeed, both in Figures 13 and14, the discharges along the spatial domain are weaker for 3 GS, owing to the presence of the three GS which slow down the ow. Moreover, this eect is more signicant for the SRE because of the saturation of the soil. Therefore, for processes like soil erosion, which are inuenced by the maximal discharge, the main result of Figure 13 is that the 3 GS conguration reduces (especially for SRE) the occurrence of high values for max t Q. Moreover a relevant information obtained with the stochastic approach is that, for SRE, (resp. LRE) the 90th percentile is 33% (resp. 11%) lower with the 3 GS conguration Concerning the sensitivity analysis, for the four congurations (1 GS or 3 GS; SRE or LRE), the rst-order sensitivity indices related to the eld (in the range 92% to 96%) are much higher than those related to the GS. This shows that only the K s of the eld is an inuent parameter, owing to the very important inltration capacity of the GS.

To study the eect of the values for K s , we have also tested the Grass strip(s) test case with less inltrating GS. The obtained results corroborate the previous observations. There is no signicant dierence in terms of runo and discharge at the outlet, but the presence of three GS slows down the ow and diminishes the occurrence of extreme values for the ow rates.

Conclusion

In this work, we have studied the impact of the variability in soil properties on overland ows caused by rainfall events. We have considered the soil saturated hydraulic conductivity K s as the most uncertain input parameter in the framework of the GreenAmpt inltration model. To model uncertainties, the ow domain has been split into subdomains reecting the spatial organization of the landscape (e.g., agricultural elds, grass strips), and the saturated hydraulic conductivity has been described by statistically independent and uniformly distributed random variables, with one random variable assigned to each subdomain. Concerning output quantities, we have focused on the discharges at the outlet (peak runo rate and runo coecient). Two test cases, named Three-eld and Grass strip(s), have been investigated.

The Three-eld test case aims at understanding the role of spatial organization in uncertainty propagation. The conclusions depend on the level of soil saturation. For long rainfall events leading to highly saturated soils, the variability of model outputs remains moderate. Moreover, the most inuent input parameter is the K s taking the highest values, except when the most inltrating subdomain is located upslope, in which case the most inuent input parameter is the K s taking intermediate values. For short rainfall events with moderately saturated soils, the most inuent input parameter, regardless of its value, is the K s located downslope, that is, the closest to the outlet. Moreover, since the variability of model responses is high, a practical conclusion in this situation is that additional measurements near the outlet should be the most eective to help reduce output uncertainties.

The Grass strip(s) test case aims at comparing, in view of possible land management issues, runo uncertainties obtained with two possible spatial localizations of grass strips within a eld, namely three narrow, equally-spaced grass strips versus one large grass strip located at the eld outlet. The rst conclusion is that the duration of the rainfall event substantially impacts the shape of the probability density function (pdf ) of the model outputs. Specically, highly peaked pdf 's are obtained for short rainfall events (and moderately saturated soils), while relatively at pdf 's are obtained for long rainfall events (and highly saturated soils). The second conclusion is that the localization of the grass strips does not impact the variability of the model outputs. However, one important dierence concerns the spatial distribution of maximal discharges since the conguration with three grass strips leads to less probable extreme values, as reected by the lower values taking by the 90th percentile for the conguration with three grass strips. This observation is relevant in view of assessing erosion risks, since the detachment of soil particles is very sensitive to the peak discharge.

To sum up, the present methodology can be used as a tool to determine where additional eorts should be concentrated when collecting input parameters to reduce output uncertainties, and also as a tool to assess risks due to extreme values, keeping in mind that quantitative results depend on the probabilistic model selected for the input parameters.

Possible further applications can target soil erosion based on sediment transport models and groundwater contamination based on pollutant transport models.

[1] E. Audusse. Modélisation hyperbolique et analyse numérique pour les écoulements en eaux peu profondes. PhD thesis, Université Pierre et Marie Curie -Paris VI, 2004.

Figure 1 :

 1 Figure 1: Geometric conguration and basic notation.

  imposed during a time T R [T] and stopped afterwards. The simulation time is denoted by T [T]. Two congurations are studied for the spatial distribution of the uncertainty parameter K s . For each spatial conguration, two rainfall events are simulated, a short rainfall event (SRE) and a long rainfall event (LRE).

Figure 2 :

 2 Figure 2: Three-eld test case with the spatial localization + -o: (a) initial conguration ; (b) rainfall hydrograph for the impermeable conguration and the case where K s,+ , Ks,o, and K s,-all take their respective mean value.

Figure 4 (

 4 b) presents the MC estimate of the expectation and standard deviation of Q max with ±3 bootstrap standard error bounds plotted against the sample set dimension M . A sample set dimension equal to 100,000 appears to be sucient to achieve convergence, and in agreement with the previous convergence of the sensitivity indices. This value for M is used in what follows. Figure4(c) illustrates the inuence of the bandwidth η on the pdf estimation. An under-smoothed pdf is obtained with a small value (η = 0.01) whereas an over-smoothed pdf is obtained with a large value (η = 0.5). The value η = 0.05 yields a suciently smoothed pdf without spurious oscillations. This value for η is used in what follows.

Figure 3 :

 3 Figure 3: Grass strip(s) test case: initial conguration.

Figure 4 :

 4 Figure 4: Numerical verication for the Three-eld test case with the spatial localization + -o: (a) convergence of the rst-order and total sensitivity indices for Qmax as a function of sample dimension M ; (b) convergence of the MC estimate of Qmax expectation and standard deviation with ±3 standard error as a function of sample dimension M ; (c) probability density estimation of Qmax, using the kernel density estimator for dierent bandwidth values η.

Figure 5 :

 5 Figure 5: Three-eld test case and Short Rainfall Event: peak runo rates Qmax and runo coecients RC for the six possible spatial localizations of soil types.

Figure 6 :

 6 Figure 6: Three-eld test case and Long Rainfall Event: peak runo rates Qmax and runo coecients RC for the six possible spatial localizations of soil types.

Figure 7 :

 7 Figure 7: Three-eld test case: condence interval of the ratio Ic/Ks at nal time as a function of spatial position and for the six possible localizations of soil types; Short Rainal Event (in grey) and Long Rainfall Event (in black).

- 3 Table 5 :Figure 8 :

 358 Figure 8: Three-eld test case: probability density function of the peak runo rate Qmax estimated with a bandwidth η = 0.05; linear scale (top) and logarithmic scale (bottom); (a) Short Rainfall Event (SRE); (b) Long Rainfall Event (LRE).

Figures 9 (Figure 9 :

 99 Figures 9(a) and 9(b) present the rst-order sensitivity indices estimated for the 100,000 standardized samples and for the two rainfall events. The sensitivity indices related to K s,+ , K s,o , and K s,-are respectively denoted by S + , S o , and S -. Concerning SRE (Figure 9(a)),

Figure 10 :

 10 Figure 10: Grass strip(s) test case with Short and Long Rainfall Events (resp. SRE and LRE): peak runo rates Qmax and runo coecients RC for the two congurations (one large grass strip downslope (1 GS), or three narrow grass strips (3 GS)).

Figure 11 :

 11 Figure 11: Grass strip(s) test case: condence interval of the ratio Ic/Ks at nal time as a function of spatial position for the two congurations (one large grass strip downslope (1 GS), or three narrow grass strips (3 GS)); Short Rainfall Event (in grey) and Long Rainfall Event (in black).

Figure 12 :

 12 Figure 12: Grass strip(s) test case: probability density function of the peak runo rate Qmax for the Short and Long Rainfall Events (SRE and LRE), estimated with a bandwidth η = 0.05; linear scale (top) and logarithmic scale (bottom).

Figure 14 :

 14 Figure 14: Grass strip(s) test case: statistical estimations of the peak runo rate maxt Q(x, t) as a function of spatial position (mean µ = Ê[maxt Q(x, t)]), median P 50 , and 90th percentile P 90 ), and some deterministic values of this quantity (taking Ks = K min s or K max s ) for the Long Rainfall Event.

Table 1 :

 1 Three-eld test case, data for the two rainfall events: rainfall intensity R, rainfall duration T R , and total simulation time T .

		-	o	+
	K min s	2.78•10 -7 2.78•10 -6 1.10•10 -5
	K max s	1.10•10 -6 5.50•10 -6 1.66•10 -5

Table 2 :

 2

Three-eld test case: minimal and maximal values of Ks (m.s -1 ) for the soil types.

Table 3 :

 3 Grass strip(s) test case, data for the two rainfall events: rainfall intensity R, rainfall duration T R , and total simulation time T .

		eld	GS
	K min s	3.57•10 -6 2.22•10 -5
	K max s	6.35•10 -6 3.33•10 -5

Table 4 :

 4 Grass strip(s) test case: minimal and maximal values of Ks (m.s -1 ).

  Moreover, we observe that for both SRE and LRE, cov takes higher values for RC than for Q max .

			Short Rainfall Event				Long Rainfall Event
		Ks,-downslope	Ks,o downslope	Ks,+ upslope	Ks,+ midslope	Ks,+ downslope
		+ o -	o + -	+ -o	-+ o	+ o -	+ -o	o + -	-+ o	o -+	-o +
	µ	1.8•10 -5	1.8•10 -5	1.5•10 -5	2.6•10 -7	4.0•10 -5	4.0•10 -5	3.8•10 -5	3.9•10 -5	3.9•10 -5	3.9•10 -5
		2.8•10 -4	2.8•10 -4	8.7•10 -5	5.9•10 -7	8.2•10 -3	8.1•10 -3	6.9•10 -3	6.2•10 -3	6.2•10 -3	6.2•10 -3
	σ	1.8•10 -6	1.8•10 -6	6.1•10 -6	5.6•10 -7	1.8•10 -6	1.8•10 -6	3.4•10 -6	3.3•10 -6	3.3•10 -6	3.3•10 -6
		6.3•10 -5	6.3•10 -5	6.1•10 -5	1.5•10 -6	5.2•10 -4	5.4•10 -4	6.6•10 -4	8.2•10 -4	8.3•10 -4	8.3•10 -4
	cov	10%	10%	42%	217%	5%	5%	9%	9%	9%	8%
		22%	22%	70%	250%	6%	7%	10%	13%	13%	13%
	P50	1.8•10 -5	1.8•10 -5	1.5•10 -5