A physically-based parsimonious hydrological model for flash floods in Mediterranean catchments

A spatially distributed hydrological model, dedicated to flood simulation, is developed on the basis of physical process representation (infiltration, overland flow, channel routing). Estimation of model parameters requires data concerning topography, soil properties, vegetation and land use. Four parameters are calibrated for the entire catchment using one flood event. Model sensitivity to individual parameters is assessed using Monte-Carlo simulations. Results of this sensitivity analysis with a criterion based on the Nash efficiency coefficient and the error of peak time and runoff are used to calibrate the model. This procedure is tested on the Gardon d’Anduze catchment, located in the Mediterranean zone of southern France. A first validation is conducted using three flood events with different hydrometeorological characteristics. This sensitivity analysis along with validation tests illustrates the predictive capability of the model and points out the possible improvements on the model’s structure and parameterization for flash flood forecasting, especially in ungauged basins. Concerning the model structure, results show that water transfer through the subsurface zone also contributes to the hydrograph response to an extreme event, especially during the recession period. Maps of soil saturation emphasize the impact of rainfall and soil properties variability on these dynamics. Adding a subsurface flow component in the simulation also greatly impacts the spatial distribution of soil saturation and shows the importance of the drainage network. Measures of such distributed variables would help discriminating between different possible model structures. Correspondence to: H. Roux (helene.roux@imft.fr)


Introduction
Flash floods are defined as sudden floods with high peak discharges produced by severe thunderstorms that are generally of limited area extent (IAHS- UNESCO-WMO, 1974). They represent one of the most destructive hydrological hazards in the Mediterranean region and have caused billions of euros of damages in France over the last two decades (Gaume et al., 2004). Flash flood prediction and risk assessment still lack efficient procedures, mainly because these events are poorly monitored and understood (Marchi et al., 2010). Effectively, this kind of hydrological event is often concerned with poorly gauged small catchments. Moreover, hydrometric stations are vulnerable to destruction or damage in case of flooding and the data flow continuity requirement may not be satisfied. Gaume et al. (2009) and Gaume and Borga (2008) also point out that flash floods are poorly documented phenomena and propose a methodology to increase the existing knowledge of such events using post-flood surveys.
In the Mediterranean climatic zone, precipitation is highly variable, both in time and space, and this variability increases with elevation in mountainous regions (Moussa et al., 2009). Norbiato et al. (2007) studied a flash flood generating storm occurring in the eastern Italian Alps in August 2003 and found extreme spatial gradients up to 80 mm km −1 in precipitation accumulations on 12 h. As a consequence, hydrological processes are also highly variable and therefore difficult to predict (Moussa et al., 2009). Intense, short-duration precipitation rates usually cause flash flood events primarily dominated by overland flow (Gaume et al., 2003).
Modeling the hydrological response of Mediterranean catchments has been explored by several authors who emphasized the importance of the topic. Piñol et al. (1997) Published by Copernicus Publications on behalf of the European Geosciences Union. applied TOPMODEL to simulate the hydrological response of two catchments located in northeastern Spain. The authors emphasized that the major difficulties encountered for such catchments consists in the spatial soil depth heterogeneity as well as the characteristics and the localized nature of downslope flows of water in the soil. They also pointed out that models with very large numbers of parameters would not be easy to calibrate. Several other applications of TOPMODEL in the Mediterranean environment have already been made (Datin, 1998;Durand et al., 1992;Saulnier, 1996;Saulnier and Le Lay, 2009), and lead to promising results even though improvements are needed in modeling wetting-up periods or extreme storm events. Moussa et al. (2009) proposed a spatially distributed model based on simplified physical process representations of the water cycle and flood genesis. This model is well adapted to event-based simulations of floods when surface runoff is the main hydrological process. However, their results show that when the calibrated model parameters were constrained to simulate intense flood events, then the model's performance decreased during years with low rainfall rates. Thus the question arises of the domain and limit of application of the model. It seems difficult to well represent hydrological processes during both drought and flood periods. Blöschl et al. (2008) developed a distributed flood forecasting system that they evaluated over the 622 km 2 Kamp catchment (Austria). They adopted a modeling strategy based on (i) a model structure defined at the model element scale and (ii) multi-source model identification and verification. Within the proposed model structure, the adjustment of 21 parameters for each pixel, 6 routing parameters for each catchment and 7 stream routing parameters for each river reach is required. On the studied catchment, there was a total of 1550 pixels, 12 sub-catchments and 10 stream reaches. The authors relied on the interpretability of model parameters and on their multi-source model identification procedure to facilitate the estimation of model parameters in a realistic way. In their study, they used runoff data along with spatial data including piezometric heads or inundation patterns, both from satellite and ground-based data. However, these are not easily available data and they might not be accessible in the case of poorly-gauged catchments. The study of Reed et al. (2007) was flash flood forecasting dedicated, too: their model produces high-resolution grids of peak flow forecast frequencies during rainfall events. Forecasters can therefore compare these grids to locally derived threshold frequency grids to aid in warning decisions. Their model can be implemented using any continuous simulation distributed hydrologic model. However, if this approach can be helpful for real time forecasting, it doesn't lead to any improvement on the understanding of flash flood hydrological processes.
In this study, rather than adopting directly an existing model and simply going through a calibration/validation exercise, we propose an approach based on building the structure of the model using the understanding of the Mediterranean catchments' hydrological processes. The hydrological rainfall-runoff model MARINE (Modélisation de l'Anticipation du Ruissellement et des Inondations pour deś evéNements Extrêmes) introduced in this paper aims at (i) exploiting the potential of distributed models (ii) using physically meaningful parameters while (iii) maintaining a simple and parsimonious parameterization. The present paper explains the overall structure and methodology of the MARINE model, the infiltration, the subsurface and surface runoff calculation, the data requirements and the model calibration procedure. The MARINE model is applied to the Gardon d'Anduze basin (southern France) to clarify the model data and calibration requirements together with its flood forecasting capabilities. In the present study, water transfer through the subsurface zone was first assumed too slow to contribute to the hydrograph response to an extreme rainfall event and therefore the MARINE model was run without activation of the subsurface lateral flow component. In a second step, the impact of including it in the modeling is discussed.

Structure and methodology of the MARINE model
The modeling approach followed herein consists in building a distributed model for flash flood forecasting. The predominant factor determining the formation of runoff is represented by the topography: slope and downhill directions. MARINE runs on a fixed time step and is structured into three main modules (Fig. 1). The first module allows separating the precipitation into surface runoff and infiltration; the second module represents subsurface downhill flow, and the third one the overland and channel flows: the transfer function component allows routing the rainfall excess to the catchment outlet using different approximations of Saint-Venant equations. Both infiltration excess and saturation excess are represented within MARINE. The spatial discretization of the catchment is performed using the Digital Elevation Model grid resolution, a regular grid of squared cells. Evapotranspiration is not represented since the model's purpose Nat. Hazards Earth Syst. Sci., 11, 2567-2582, 2011 www.nat-hazards-earth-syst-sci.net/11/2567/2011/ was to simulate individual flood events during which such process is negligible. The model simulates the flood hydrograph at any point of the drainage network. It is also possible to follow the evolution of the distributed variables such as soil moisture or overland flow velocities all over the catchment. A description of each procedure is detailed hereafter.

Infiltration
Local scale infiltration is described assuming onedimensional flow to occur into independent, vertically homogeneous soil columns using the Green and Ampt model. Infiltration rate is equal to rainfall intensity as long as rainfall intensity doesn't exceed potential infiltration rate. When rainfall rate exceeds the infiltration rate, ponding occurs (Gourley and Vieux, 2006).
The infiltration rate i (m s −1 ) is given at the local scale by: where r is rainfall rate (m s −1 ), t p is time to ponding (s), K is the saturated hydraulic conductivity (m s −1 ), ψ is the soil suction at wetting front (m), θ s and θ i are saturated and initial water contents (m 3 m −3 ) respectively and I is the cumulative infiltration (m). I is calculated from rainfall rates accumulated over time and i(t) and I (t) are related by: i (t) = dI (t) dt . Infiltration stops when soil water content θ (m 3 m −3 ) exceeds saturated water content θ s .

Subsurface flow
The subsurface model is based on Darcy's law. Using frequently invoked assumptions, (i) the slope of the water table in the saturated zone is assumed to coincide with local topographic slope, (ii) the local transmissivity is an exponential function of the local storage deficit (Original TOP-MODEL assumption, Beven and Kirkby, 1979), the flow per unit width q is expressed as: where T 0 is the local transmissivity of fully saturated soil (m 2 s −1 ), θ s and θ are saturated and local water contents (m 3 m −3 ), m is transmissivity decay parameter (-), and β is local slope angle (rad). When soil water reaches the drainage network, it is assumed that the flow into the drainage network occurs with a velocity calculated using Eq.
(2). This subsurface model represents the horizontal flow in the unsaturated zone which usually occurs in a layer of limited thickness and with high hydraulic conductivity due to the preferential flow paths and macroporosity (Ciarapica and Todini, 2002).

Surface water
The surface runoff calculation is divided into two parts: the overland flow and the flow along the drainage network. Both are simulated using the 1-D kinematic wave approximation of the Saint-Venant equations with the Manning friction law.

Overland flow
To represent the overland flow, it is assumed that conditions hold for the application of the kinematic model. The onedimensional water mass balance equation for the overland flow along a uniform slope is described as follows: where h is water depth (m), t is time (s), u is overland flow velocity (m s −1 ), x is space variable (m), r is rainfall rate (m s −1 ), and i is infiltration rate (m s −1 ). In the kinematic wave analogy, the momentum equation reduces to S 0 = S f , where S 0 stands for bed slope (m m −1 ) and S f for energy gradient line (m m −1 ). The Manning friction law provides a relationship between flow depth h and flow velocity u: where n o is the Manning roughness coefficient (m −1/3 s). Along with the mass conservation Eq. (3), this law allows simulating the overland flow: The forcing function on the right hand side of Eq. (5) expresses the rainfall excess that is the difference between the rainfall rate and the soil infiltration rate. Soil infiltration is treated by the Green and Ampt equation as explained above.

Drainage network
When the drainage area becomes greater than 1 km 2 , the overland flow is structured in a drainage network. Flow in this drainage network is simulated using the kinematic wave  (Rawls et al., 1992). approximation by taking into account the width of the network reach. A hypothesis is made on the network reach cross-section geometry (Fig. 2). This leads to a new transport equation in the drainage network. Characteristics of the network reach are calculated using geomorphological considerations (Liu and Todini, 2002): network reach width W Di and depth H Di are taken to increase as a function of the area drained by the ith cell, such that : where W Dmax is the maximum width at the basin outlet, W Dmin is the minimum width, corresponding to the threshold area estimation of Giannoni et al. (2000) a th =1 km 2 , which is the minimum upstream drainage area required to initiate a channel, a tot is the total area and a di is the area drained by the ith cell.

Parameter description and data requirements
Parameterization of the MARINE model and the required input data to run the MARINE model without activation of the subsurface lateral flow component are described below.

Parameter description
On each cell, the model needs the following parameters: (i) 5 parameters for the estimation of the surface runoff and infiltration, namely the saturated hydraulic conductivity K, the saturated and initial water contents θ s and θ i , the soil suction at wetting front ψ, and the soil thickness Z, (ii) 2 parameters for the calculation of overland flow if the cell is not in the drainage network, namely the local slope S 0 and the surface roughness n o , (iii) 6 parameters for the transfer function in the drainage network, namely the hillslope S 0 , the depth H D , the width W D and the cross-sectional slope S D of the network reach, and two roughness coefficients n D1 and n D2 (river bed and flood plain Manning coefficient). Coefficient n o varies from 0.03 up to 0.1 m −1/3 s and n D1 varies from 0.025 up to 0.05 m −1/3 s. Hence the model needs 7 parameters for a hillslope cell and 11 parameters for a cell in the drainage network. Most of these parameters can be estimated using information on topography, soil and land cover as explained later in Sects. 3.2 and 4. When a separate modeling of the channel hydrodynamic is required, that is to say when the kinematic wave assumption doesn't hold in the river (Moussa and Bocquillon, 1996), the mass transfer model requires more knowledge about the geometric characteristics of the river (cross sections geometry, roughness, hydraulic structures).

Data requirements
The MARINE model requires field data, usually from Digital Elevation Models (DEM), soil survey and vegetation or land-use, as well as precipitation measurements. The DEM application consists in identifying connections between cells, thereby giving the catchment extent and the flow pathways, calculating the hillslope and cumulating the drainage area for detecting the drainage network and calculating the geometric characteristics of the network reaches. In the MARINE model, as in many rainfall-runoff models (Liu and Todini, 2002), drainage is only possible along the cardinal directions for the four adjacent cells at each edges. The overbank crosssectional slopes S D of each network reach are also derived from the DEM. The S D parameter is approximated by the slope between the cell in the network reach and the neighbouring cells in the catchment. Soil surveys, provided by INRA and BRGM, allow deriving soil texture and thickness. Soils are assumed to be vertically homogeneous. According to the soil texture, a Rawls and Brakensiek (1983) soil class is assigned to each cell (Fig. 3). For each soil class, the estimated values for the Green and Ampt front suction, saturated hydraulic conductivity and porosity are determined according to the results of Rawls and Brakensiek (1983). The vegetation and land-use map (2000 Corine Land Cover: Service de l'Observation et des Statistiques) is used to derive distributed surface roughness (Chow, 1959).

Model calibration and formulation of calibration criteria
The chosen approach for model calibration and validation followed by the sensitivity analysis procedure is described in the following section.

Estimated parameters and calibration procedure
Parameterization and calibration is a crucial issue for hydrologic models. As it is physically based, the MARINE model can take advantage of the information on topography, soil characteristics and land cover for parameter estimation. In order to avoid a model over-parameterization, the number of parameters to estimate was kept as low as possible. According to Refsgaard's (1997) recommendations, this can be done by fixing a spatial pattern of a parameter and allowing its absolute value to be modified by calibration. This approach has been chosen for two parameters, namely the distributed saturated hydraulic conductivities K and soil thicknesses Z.
The spatial patterns of these parameters are derived from soil surveys and a unique coefficient of correction is then applied to each parameter map. The calibration procedure consisted in estimating (i) these two coefficients of correction: one for the saturated hydraulic conductivities, named C K and the other one for the soil thicknesses, named C Z , (ii) the overbank roughness of the drainage network n D2 , (iii) the initial soil saturation condition θ i . The choice of these calibration parameters followed observations made during a first calibration process carried out manually using a trialand-error procedure. As a matter of fact, the model was not very sensitive to some parameters, especially to both other roughness coefficients n o (surface roughness) and n D1 (main channel roughness of the drainage network), therefore they were chosen according to vegetation and land-use information and n D1 was kept constant all over the catchment. The estimation of C K , C Z and n D2 has been implemented for one flood event and then, as a validation procedure, the estimated values were used to simulate other events that occurred in the same catchment. The initial soil saturation condition θ i had to be set for each event. Eventually, only four parameters needed to be calibrated for the whole catchment.
Calibration of the MARINE model has more to do with an adjustment than with a conventional calibration and could have been carried out by a simple trial-and-error method. However, in order to be able to qualify model outputs, the chosen procedure was to achieve Monte-Carlo simulations to derive the sensitivity of the model to individual parameters and to determine a calibrated parameter set defined as the set giving the best simulated hydrograph according to a chosen criterion. Calibration criterion and sensitivity analysis procedure are introduced in the following paragraphs.

Model performance criterion and sensitivity analysis
The first step of the sensitivity analysis consists in the definition of a likelihood measure intended as an evaluation of how well the model conforms to the observed system behavior. The possibility of including various types of criteria into the likelihood measure makes the concept attractive for evaluating reliability in flood extent modeling, as demonstrated in several examples (Aronica et al., 1998(Aronica et al., , 2002Romanowicz and Beven, 2003;Werner, 2004). Well known performance criterion functions have been used to build a criterion evaluating the performance of the MARINE model: it consists in a linear combination of the efficiency coefficient (Nash and Sutcliffe, 1970) and the error of peak time and runoff (Lee and Singh, 1998): where N is the number of observation data, Q s and Q o are respectively the simulated and the observed runoff, Q s P and Q o P are respectively the simulated and observed peak runoff, T s P and T o P are respectively the simulated and observed time to peak, T o C is the time of concentration of the catchment. The error of peak time and runoff is designed to aid in warning decisions in emphasizing peak flow characteristics. That is one of the attractive aspects of this method: in formulating the likelihood measure, explicit thought must be given to how model performance is assessed in the light of model application. This L NP criterion is therefore an attempt to conciliate real time flood forecasting requirements with a better understanding of the physical phenomena involved in flood event generation. It has been calculated only for specific observed discharges Q o greater than 0.3 m 3 s −1 km −2 at the Anduze station since the aim was to focus on reproducing extreme events and very high flows. Therefore, smaller discharges were neglected in the evaluation of the goodness of simulated hydrographs.
In order to derive the sensitivity of the model to individual parameters, Monte-Carlo simulations were achieved by running the model with different randomly chosen sets of parameter values. Initial range of parameter values to be considered is selected with the intent of preserving physically realistic parameter values. Uniform parameter distributions within their range of variation are mainly used in lack of prior information. Each set of parameter values is then assigned a likelihood of being a simulator of the system, on the basis of the chosen likelihood measure (Eq. 7).

Model implementation on the Gardon d'Anduze river basin
Since its advent in 2000, the MARINE model has been applied to several catchments for uses such as flood forecasting or extreme flood analysis. For instance, it has been implemented on the Thoré river (Goutorbe et al., 2002) (Tanguy et al., 2005).

Catchment characteristics
The Gardon d'Anduze river is located in southern France, 70 km northeast of the city of Montpellier. The catchment drains an area of 545 km 2 . The river flows in a southeast direction to the confluence with the Rhône river. Over its course, the Gardon d'Anduze river is joined by tributaries including the Gardon de Sainte Croix, Gardon de Mialet and Gardon de Saint Jean (Fig. 4). The river water course has a total length of approximately 50 km. Local climatic tendencies produce the highest flooding risk in autumn with the maximum rainfall rate in this period. Summers are hot and dry; however summer storms can also present a nonnegligible flooding risk. The W Dmin parameter (Eq. 6) is equal to the nearest value of drain width corresponding to intermittent flow process in the region and is set to one meter according to field observation. Field observations also led us to propose that W D varies from 1 m to 30 m, and H D varies from 0.1 m up to 2 m.

Topography
The catchment has a highly marked topography consisting of mountain peaks, narrow valleys and steep hillslopes. The highest areas are found in the Cévennes, where the elevation rises till 1200 m a.s.l. near the mount Aigoual. The river basin elevation at Anduze is approximately 130 m. A DEM data file of the study site with a grid scale of 50 m was available from the National Geographic Institute (IGN -BD TOPO®) (Fig. 4). The mean slope of the whole river basin is approximately 20 %.

Soils
About 64 % of the catchment area develops on metamorphic terrain. The substrate is made of shale and crystalline rocks overlain by silty clay loams -on 83 % of the area -and sandy loam top soil (Moussa and Chahinian, 2009

Vegetation and land use
Vegetation is dense and mainly composed of chestnut trees, pasture, Holm oaks, conifers, waste land and garrigue. Chestnut trees are located in the upstream area and on the south-facing slopes (sunny sides or adret) while forested garrigues and Holm oaks are located in the downstream area and on the north-facing slopes (shady sides or ubac). A vegetation and land-use map (2000 Corine Land Cover provided by the Service de l'Observation et des Statistiques (SOeS) of the French Ministry of Environment, www.ifen.fr) was used to derive distributed surface roughnesses.

Hydrometeorological data availability
Radar rainfall measurements combined with rain gauge data have been available since 1994, with a 6 min time step from 1994 to 2001, and 5 min since 2002. Measurement grid spatial resolution is 1 km by 1 km. Several floods have been experienced in this catchment since 1994. Four events, occurring in 1994, 1995, 2000 and 2002, were retained according to their representativeness of the various hydrological behaviors observed in the basin. The flood of 2002 is an exceptional event with a return period of more than 50 yr. This maximal water level at Anduze has been exceeded in 1907 and 1958. The events of 1994, 1995 and 2000 are medium ones but the flood of 1994 is the consequence of two distinct rainfall events separated by a 30-h interval while the flood of 1995 occurred in October. October often presents high soil moisture, but it also occurs that the first rain after summer takes place in November. The total rainfalls range between 187 mm for the flood of 1995 and 297 mm for the flood of 2002, the runoff coefficients -that is to say the ratio of total streamflow volume to the total precipitation over the catchment area for the considered event -between 24 % for the flood of 2000 and 48 % for the flood of 1995 and the maximum discharges between around 800 m 3 s −1 for the flood of 1994 and around 3600 m 3 s −1 for the flood of 2002. Characteristics of the studied flood events are summarized in Table 1 and total cumulated rainfalls are shown in Fig. 6.

Sensitivity analysis and model calibration
Sensitivity analysis was achieved on two flood events: the intermediate flood of 1994 and the exceptional flood of 2002.
The aim was to compare the model sensitivity to individual parameters for different kinds of flood. Sensitivity to the estimated 4 parameters was tested: the two correction coefficients -C K for the saturated hydraulic conductivities and Manning roughness coefficient of the overbanks (m −1/3 s) 0.06 1 0.2 C Z for the soil thicknesses -the overbank roughness of the drainage network n D2 and the initial soil saturation condition θ i . Parameter variation ranges are listed in Table 2. Prior parameter distributions have been chosen uniform. The scatter plots of Fig. 7 correspond to the results obtained for the flood of 2002 but the results obtained for the flood of 1994 are quite similar. The values of simulated discharges seem to be very sensitive to the overbank roughness of the drainage network n D2 . In Fig. 7, it can be seen that good and poor simulations are available throughout the same range for the 4 parameters. It suggests that the parameter response surface is very complex and confirms that the value of one single parameter has little meaning when taken outside the context of the other parameter values. Parameters C K , C Z and θ i show a large range of equifinality. However, low values of the correction coefficients C K and C Z and high values of the initial soil moisture content θ i always correspond to negative likelihood values. There seems to be a threshold in those three parameter values: above this thresholdor below this threshold for θ i -the model is less sensitive to the chosen parameter. These parameters govern the infiltration mechanism, they determine how much water is infiltrated and at which rate: they therefore control the value of the simulated runoff coefficient.
The model was calibrated on the basis of these Monte-Carlo simulations using the medium flood of 1994. The parameter set giving the best simulated hydrograph for the likelihood measure L NP is described in Table 2. With these parameter values, the soil depths range from 0 m to 5.3 m with an average of 1.5 m and the hydraulic conductivities range from 9 mm h −1 to 101 mm h −1 with an average of 61 mm h −1 . The next step was to test the calibration by simulating other flood events with the parameters estimated from the flood of 1994.

Model validation
Using the parameter values of Table 2, simulations were carried out for the floods occurring in 1995, 2000 and 2002. Comparison of observed and simulated hydrographs for all floods is shown in Fig. 8. Table 3 summarizes the simulation results: they show a good agreement with the corresponding observations for all three events, except for the simulated runoff coefficient and peak discharge of the 1995 event. Indeed, the floods of 1994, 2000 and 2002 occurred in September and show similar runoff coefficients, below 40 %. The flood of 1995 occurred in October and has therefore a greater runoff coefficient. If it seems to be relevant to use the same initial soil moisture content θ i for the floods of 1994, 2000 and 2002, this is not the case for the flood of 1995. This may explain the lower value of the L NP criterion for this event: 0.57 against 0.94 for the 2000 flood and 0.91 for the 2002 flood. Indeed, using the same parameters but a higher value of the initial soil moisture contentθ i = 78 % -the goodness criterion value increases to L NP = 0.82 as it can be seen in Fig. 9. This emphasizes the need for a model initialization related to the time of occurrence of each event.

Simulated hydrographs at upstream locations
Output hydrographs are available at any point of the drainage network, as it can be seen in Fig. 10. It is consequently possible to follow the evolution of the discharge along the drainage network. Therefore, the hydrographs simulated using the parameters estimated from the observed discharges Nat. Hazards Earth Syst. Sci., 11, 2567-2582, 2011 www.nat-hazards-earth-syst-sci.net/11/2567/2011/  of the 1994 event at Anduze have also been compared with the observations made at Saumane and Mialet, two discharge gauging stations located upstream of the Anduze station (Fig. 4). Results are shown in Tables 4, 5 and Fig. 11.
The hydrographs at the Mialet station present poor L NP values as discharges are systematically overestimated. Tributaries on this area of the catchment are indeed significantly affected by karst processes: part of the discharge probably flows in fissure-karstic paths and is therefore not measured by the gauging station at Mialet. Simulated hydrographs at the Saumane station show a better agreement with observed discharges, especially for the 2000 event with L NP = 0.57 against L NP = −0.28 at Mialet and L NP = 0.94 at Anduze. However, even at the Saumane station, simulated discharges are less satisfactory than at the Anduze station. Indeed, the estimated parameters enabling the model to reproduce the integrated response at Anduze and the Saumane subcatchment characteristics differ from the ones of the entire catchment: for instance, the mean soil depth is approximately 1m for Saumane catchment against 1.5 m for the Gardon catchment. This emphasizes the need for regionalization methods and the assessment of the relationship between local parameter identifiability and catchment characteristics (Wagener and Wheater, 2006). This shows the impact of rainfall spatial distribution on saturation dynamics as the major part of the 2002 rainfall event occurring on the downstream part of the catchment as it can be seen in Fig. 6, whereas it occurred on the upstream part in 1994. Moreover, in the upstream part of the catchment, soil depths are lower than in the downstream part as it can be seen in Fig. 5a).

Soil saturation dynamics
To emphasize the importance of soil properties spatial distribution on the evolution of soil saturation dynamics, the mean saturation state of 1994 event has been calculated for 7 soil categories combining Rawls and Brakensiek's (1983) Nat. Hazards Earth Syst. Sci., 11, 2567-2582, 2011 www.nat-hazards-earth-syst-sci.net/11/2567/2011/ soil classes and soil depths existing in the catchment (Fig. 13). It can be seen that soil properties' spatial variability has a great impact on this dynamic. Soils with high hydraulic conductivities (classes 3, sandy loam, and 4, loam) and low depths are rapidly saturated: the mean saturation state for C4, depth 0-1 m increases from 58 % at the beginning of the event till 93 % in only 12 h to reach 100 % after 36 h. Soils with low hydraulic conductivity and high depth (class 7, silt and depths ranging between 4 m and 5 m) exhibit little dynamic: mean saturation state at the end of the event is of 66 % against 58 % at the beginning.

Impact of subsurface flow
The influence of including subsurface flow in the simulations has been tested using the same parameters of Table 2 and adding a subsurface flow component with an horizontal hydraulic conductivity presenting the same spatial distribution as the vertical one but 1000 times greater. As it can be seen in Fig. 14a for the 1994 event, the resulting hydrograph modifications are most important at the beginning of the rainfall event, between the two peaks and during the recession period. Assuming that water transfer through the subsurface zone is too slow to contribute to the hydrograph response to an extreme event may not be true for periods with a lower rainfall rate. Indeed, simulations including subsurface flow show better agreement with observed discharges for these periods (Fig. 14). The likelihood criterion L NP is The inclusion of subsurface flow also greatly modifies the soil saturation dynamics. Figure 15 clearly shows the importance of the drainage network on saturation dynamics when subsurface flow is activated. This is due to the fact that exfiltration can occur in the drainage network. Measures describing the spatial distribution of saturation state would be helpful to choose which physical phenomenon should be included in the model.

Discussion and conclusions
The MARINE model is structured around the understanding of Mediterranean catchment hydrological response in order to be dedicated to flash flood prediction and analysis. It was tested on the Gardon d'Anduze river basin as part of a hydrological rainfall runoff model intercomparison project lead by the French central hydrometeorological service for flood forecasting (Ministère de l'écologie du développement et de l 'aménagement durables, 2003). Model construction and assignments of prior values to model parameters were based on easily available spatial data. Only four parameters needed to be estimated in the model. Two flood events were chosen to implement a sensitivity analysis of the model prediction to these parameters using Monte-Carlo simulations. The model was then calibrated on one flood event and tested on three others. Simulation results were compared on the basis of a model performance criterion representing both efficiency and the error of peak time and runoff.  In the following section, the overall quality of the results is discussed taking into account the sensitivity analysis, the modeling hypothesis and the model structure.

Sensitivity analysis
The MARINE model aims at using physically interpretable parameters in order to facilitate their estimation. Results of the sensitivity analysis show that the model is very sensitive to the Manning roughness coefficient of the overbanks n D2 . Indeed, this parameter is related to the transfer time to the outlet and then to the peak position. As shown on the scatter plots, the model is also sensitive to the three other parameters C K , C Z and θ i as they directly affect the runoff rate by fixing the infiltration rate and the soil capacity. However it is likely that there are interactions between these parameters and this may explain that the corresponding scatter plots are not as meaningful as the one of the roughness coefficient (Fig. 7).

Calibration, validation and model structure
Consistent results were found between calibration and verification events at the Anduze station. Results concerning the 1995 event show that a possible improvement of the model concerns the initialization: the soil moisture at the beginning of each flood strongly depends on the date of occurrence of this event as it has been shown in Sect. 6.2. An initial soil moisture specification based on the SAFRAN-ISBA-MODCOU model provided by Meteo-France (Habets et al., 2008) is currently being tested (Braud et al., 2010  of initial moisture conditions prior to a flood event where no measurements are available (Tramblay et al., 2010). The hydrographs simulated at stations located upstream the station used for calibration are less satisfactory. Results emphasize the need for implementation of regionalization methods.
The MARINE model has been built for flash flood prediction and analysis on ungauged catchments. It is therefore compatible with raster-based Geographic Information Systems and may be used with spatial data sets. The choice of a distributed model was borne out by the importance of the spatial variability of rainfall and topography in the flash flood generation. When performing tests for assessing whether distributed model simulations are different from lumped model simulations under parametric and input uncertainties representative of present-day operational flow-forecasting conditions, Carpenter and Georgakakos (2006) found that a distributed model showed better performance with respect to peak flow magnitude in approximately 60 % of the events for the two study catchments, whereas the lumped model showed better performance in less than 25 % of the events. Their main conclusion is that even on the scales of current lumped operational forecasting models, distributed models offer clear performance advantages under present day parametric and input uncertainties, when used to produce ensemble streamflow simulations. On the contrary, Saulnier and Le Lay (2009) found that spatial extent of the rainfall patterns is not always of major importance. However their results show that for the 8-9 September 2002 event, the accurate geographical localization of the storm cells was needed to significantly improve the discharges simulations. Examining the impact of spatial aggregation of rainfall and soil properties on extreme flood modelling, Sangati et al. (2009) confirm that a correct rainfall volume is not enough for an accurate reproduction of flash flood events characterised by large rainfall variability. Moreover they found that the soil properties' aggregation length exerts a similar effect on peak discharge errors as increasing the rainfall aggregation length. The impact of the spatial variability of soil and rainfall description is also supported by the results presented here, especially those concerning soil saturation dynamics (Figs. 12 and 13).
Concerning the model structure, the results show that assuming that water transfer through the subsurface zone was too slow to contribute to the hydrograph response to an extreme event may not be true for the recession period in particular (Fig. 14). Adding a subsurface flow component in the simulation also greatly impacts maps of soil saturation and emphasizes the importance of the drainage network in this dynamic. Measures of such distributed variables would help discriminating between different possible model structures.
Concerning the values of C K and C Z , the multiplicative constants of soil maps' properties contain underlying physical properties. A correction coefficient C K greater than 1 may be interpreted as the accountancy of vertical macropores existence which accelerates the wetting front vertical displacement. A correction coefficient C Z greater than 1 appears necessary to simulate the total outflow volume during the flash Nat. Hazards Earth Syst. Sci., 11, 2567-2582, 2011 www.nat-hazards-earth-syst-sci.net/11/2567/2011/