A modeling approach to assess the hydrological response of small mediterranean catchments to the variability of soil characteristics in a context of extreme events

This paper presents a modeling study aiming at quantifying the possible impact of soil characteristics on the hydrological response of small ungauged catchments in a context of extreme events. The study focuses on the September 2002 event in the Gard region (South-Eastern France), which led to catastrophic flash-floods. The proposed modeling approach is able to take into account rainfall variability and soil profiles variability. Its spatial discretization is determined using Digital Elevation Model (DEM) and a soil map. The model computes infiltration, ponding and vertical soil water distribution, as well as river discharge. In order to be applicable to ungauged catchments, the model is set up without any calibration and the soil parameter specification is based on an existing soil database. The model verification is based on a regional evaluation using 17 estimated discharges obtained from an extensive post-flood investigation. Thus, this approach provides a spatial view of the hydrological response across a large range of scales. To perform the simulations, radar rainfall estimations are used at a 1 km2 and 5 min resolution. To specify the soil hydraulic properties, two types of pedotransfer function (PTF) are compared. It is shown that the PTF including information about soil structure reflects better the spatial variability that can be encountered in the field. The study is focused on four small ungauged catchments of less than 10 km 2, which experienced casualties. Simulated specific peak discharges are found to be in agreement with estimations from a post-event Correspondence to: S. Anquetin (sandrine.anquetin@hmg.inpg.fr) in situ investigation. Examining the dynamics of simulated infiltration and saturation degrees, two different behaviors are shown which correspond to different runoff production mechanisms that could be encountered within catchments of less than 10 km2. They produce simulated runoff coefficients that evolve in time and highlight the variability of the infiltration capacity of the various soil types. Therefore, we propose a cartography distinguishing between areas prone to saturation excess and areas prone only to infiltration excess mechanisms. The questions raised by this modeling study will be useful to improve field observations, aiming at better understanding runoff generation for these extreme events and examine the possibility for early warning, even in very small ungauged catchments.


Introduction
Flash floods represent the most destructive natural hazard in the Mediterranean region causing around a billion Euros of damage in France over the last two decades (Gaume et al., 2004).Flash floods are rare events that usually occur in ungauged river basins.Amongst them, small-ungauged catchments are recognized to be the most vulnerable to storms driven flash flood (Ruin et al., 2008).
Several methods for predicting flash floods in ungauged river basins are now accepted.The flash flood guidance (Georgakakos, 2006) and the discharge threshold exceedance approach (Reed et al., 2007;Younis et al., 2008) are built to give an early flash flood warning suitable to organize the civil protection.These methods rely either on conceptual or Published by Copernicus Publications on behalf of the European Geosciences Union.
C. Manus et al.: Hydrological responses associated to soil characteristics physically based hydrological models and need to be complemented with sensitivity studies to improve the understanding of the major hydrological factors associated to the flood event.
There is no unique and simple theory about the runoff production on watersheds during flash flood events.The main reason is that various processes can be involved which are usually grouped in two types of overland flow: infiltration excess (Horton runoff) or saturation excess (Dunne process).Dunne and Black (1970) suggest that saturated areas result from saturation of a surface layer (in case of a less impermeable layer below) or the whole soil profile by incident rainfall whereas others invoke the rising of groundwater tables.Due to the high heterogeneity and space variability of the watershed characteristics (land use, soil type and depth, sub-soil, local slope, upstream contributing area) and to antecedent moisture conditioning, these processes are likely to be active at the same time in various combinations (Smith and Goodrich, 2005).Latron andGallart (2007, 2008) showed the existence of both mechanisms on a small Mediterranean catchment using field survey and piezometer and tensiometer data analysis.
But flash flood events are poorly understood due to the lack of experimental sites and long-term hydrometeorological data with adequate space-time resolution (Foody et al., 2004;Delrieu et al., 2005).In particular, few data are in general available regarding discharges during flash-floods.When gauges exist, they can be seriously damaged and turn out of work during the event.If data are provided, discharge estimation is prone to large errors as stage discharge relationships are generally extrapolated far beyond the range of gauged values.Traditional gauging in such conditions is all the more impossible because it is too dangerous for operators.Furthermore, given the scarcity of these events, it is quite difficult to have a long series of events at a given location.In this context, Gaume et al. (2004) propose the use of post flood field survey after a major event.The idea is to provide a spatial view of the hydrological response across a large range of catchment scales by reconstructing peak discharge and timing of the flood using flood marks combined with simple hydraulics consideration and witness interviews.Although restricted to maximum flow, these information are in general the only one available and should be complemented by a modelling approach, using re-analysis of radar based precipitation estimates, in order to get a regional and consistent view of the event (Borga et al., 2008).
In order to gain insight into flash floods, several hydrometeorological observatories are being set up.One of them, the long-term observatory OHMCV1 , covers an area of about 160×200 km 2 in the Cévennes-Vivarais region (Fig. 1).Various numerical models, distributed or not, have been tested in the Cévennes region.The goal of these studies was the derivation of efficient prediction tools for flash flood forecasting more than towards process understanding.They were generally set up on gauged catchments.That condition allows the calibration of model parameters on part of the data and the validation of the corresponding models with independent data and/or on sub-catchments (Moussa et al., 2007).Examples can be found in Sempere-Torres et al. (1992) and Nalbantis et al. (1995) on the Gard River, Saulnier et al. (1997) and Saulnier and Datin (2004) with Topmodel (Beven and Kirky, 1979) on the Ardèche River and Ayral et al. (2005) with Althaïr model (based on Horton infiltration capacity excess) on the Gardon.The comparison between all the models shows that, once calibrated, Horton and Dune schemes reproduce equally well all the flood events.
As regards to prediction in ungauged catchments, few studies have been conducted at regional scale.Ayral et al. (2006), with the Althaïr model, and Le Lay and Saulnier (2007), with the event-based Topmodel approach, tested various levels of sophistication in the inputs and model parameters regionalization.Ayral et al. (2006) obtained a systematic overestimation of peak discharge and a satisfactory simulation of the time of the peak, when the model was used with spatially homogeneous parameters.Le Lay and Saulnier (2007) showed that the model efficiency significantly increased when the spatial variability was taken into account.Nevertheless, for some of the catchments, the mis-performances remained unexplained.The authors raised some hypotheses about the impact of the spatial structure of soil hydraulic characteristics.
Field studies have been undertaken to locally analyze the mechanisms of runoff generation, but they only cover a small part (545 km 2 ) of the whole OHM-CV observatory (3200 km 2 ) and are located on steep areas (Ayral, 2005;Marchandise, 2007).Based on hillslope infiltration experiments with either controlled rainfall inputs or real events, the experimental results showed that the soil infiltration capacity is generally very high (more than 100 mm h −1 ) and that a large proportion of the flow is generated by subsurface lateral paths.However, the obtained velocities were too small to explain the quick response of the catchment.Ayral et al. (2005) suggested that the transit time of the water within the soil might be very short and the exfiltrated water quickly mobilized in temporary surface channels.
In order to gain insight and propose hypothesis about active process controlling runoff generation, numerical simulation can provide a useful complement to observations (Piñol et al., 1997;Vivoni et al., 2007).Using a regional distributed modeling approach we propose to assess the possible role of soil variability on runoff generation and catchment response during the flash flood Gard event in September 2002 in France (Delrieu et al., 2005).A critical outcome is a map highlighting the most probable major active processes (infiltration excess or saturation excess) and their spatial variability.The results presented here will also contribute to improve To study the sensitivity of the model responses to soil characteristics, the following modeling options are chosen given the scarcity of the available information at this scale.The proposed distributed modeling approach is carried out without calibration phase.The approach is thus applicable at the scale of the whole Cévennes -Vivarais region and gives a regional view of the problem.Following the recommendation of Borga et al. (2008), radar rainfall data are used.The soil properties are extracted from the available Languedoc-Roussillon soil database.The model accounts for different soil profiles and hydraulic characteristics, as well as a routing scheme using a kinematic wave approximation of the St-Venant equation (Moussa and Bocquillon, 1996).In the following, Sect. 2 first describes the available data and the case study.Section 3 describes the model, its set up and the methodology retained for model verification.Then, results are presented in Sect. 4 and conclusions and perspectives are finally drawn, giving some recommendations about required field studies for the future.

Region of interest and case studied
The Cévennes -Vivarais region (Fig. 1a) is the Southeast part of the Massif Central (85 000 km2 ; i.e. one sixth of the France area).The relief of the southeasterly facing slope starts from the Mediterranean shore and the Rhône Valley.The altitude ranges from sea level up to 1700 m (Mount Lozère) over roughly 70 km.The main Cévennes rivers, Ardèche, Cèze, Gard, right bank tributaries of the Rhône River, and the Vidourle have a typical almost intermittent hydrological regime: very low water levels in the summer, floods occurring mainly during autumn.The above-mentioned catchments are medium size catchments (almost 2300 km 2 for the largest) with travel times of less than 12 h.The modeling approach is set up at the scale of the whole Cévennes-Vivarais region, but the paper focuses on four small ungauged catchments shown in Fig. 1b: Domazan (6 km 2 ), Rousson (12 km 2 ), Saint Quentin la Poterie (12 km 2 ) and Quissac (2 km 2 ).They are spread within the region and located both in mountainous and plain areas.This study is motivated by the analysis of Ruin et al. (2008) on the hydrometeorological circumstances that led to accidental casualties in these small ungauged catchments, that are much vulnerable to intense storms.Ruin et al. (2008) showed that these catchments do not systematically react to the first heavy rainfall.To rate the impact of soil variability on the hydrological response of small catchments, soils characteristics are extracted from the Languedoc-Roussillon soil data base (later referred as BDSol-LR), provided by the INRA 3 from the IGCS 4 program.This database gives soil information (i.e.texture, depth of horizons etc.) on pedological landscape units called Cartographic Soil Unit (CSU).Those units are established at the 1/250 000 resolution and are georefer- 3 The French National Institute of Agronomical Research 4 http://gissol.orleans.inra.fr/enced.They are composed of Typological Soil Unit (TSU) whose vertical heterogeneity is described by stratified homogeneous layers of soil.The proportion of TSUs is given within a particular CSU, but not their precise location.Each entity is described through tables providing both quantitative and qualitative information.The latter can be used to check the consistency of quantitative data.As an example, the soil depth in the whole region is given in Fig. 2a.The average depth in the studied region does not exceed 55 cm and more than 50% of the soils are shallow (depth below 50 cm).Figure 2b shows the texture of the soils encountered within the data base.They span over a large part of the textural triangle with a mean texture of 30% silt, 20% clay and 50% sand.Figure 2b also shows the large diversity of the soils encountered in the region, which should result in a large variability of soil hydraulic properties.
The meteorological and hydrological data of the September 2002 event used in this study were collected and analyzed in the framework of the Cévennes -Vivarais Mediterranean Hydrometeorological Observatory (OHM-CV).Rain intensity was derived from the radar observation at Bollène (see Fig. 1a for its location).Data are available at 1×1 km 2 resolution every 5 min (Boudevillain et al., 2006).
The meteorological event started early in the morning of the 8th of September 2002 with the formation of first convective cells over the Mediterranean Sea.Then, convection progressed northward to form inland a mesoscale convective system (MCS) over the Gard river watershed.The quasi-stationary MCS stayed over the same region until approximately the next morning, and then moved eastward together with the surface front.As emphasized by Delrieu et al. (2005), three distinct rainfall phases were observed during the episode.In phase 1, from 08:00 UTC to 22:00 UTC, 8 September 2002, the highest rainfall developed and became stationary over the plain region just north of the city of Nîmes.Maximum rainfall amounts were about 300 mm during this phase.In phase 2, between 22:00 UTC and 04:00 UTC, 9 September 2002, the rain system slowly moved northwards and northwestwards.The MCS produced rainfall amounts greater than 100 mm over this 6 h-period.Lastly, in phase 3, after 04:00 UTC, rainfall amounts exceeding 100 mm during an 8 h-period occurred during the eastward evacuation of the storm.
For the whole event, the raingauge network locally recorded 24 h cumulated rainfall greater than 600 mm, which was confirmed by radar observation (Fig. 3).
As mentioned before, river discharge information is more critical to obtain.The available water level stations cover watersheds of more than several hundreds km 2 whereas the basins of interest are around 10 km 2 .As proposed by Gaume et al. (2004) for a previous similar event, an extensive postflood investigation was carried out during the months following the event.The objective was to collect information concerning the flood and to analyze the hydrological behavior of watersheds with an area of 2 to 300 km 2 .Basically, the procedure aims at estimating maximum discharges by means of water level marks and also interviewing witnesses to document the chronology of the flood.During the post-event investigation, 93 river cross sections were surveyed and 143 witnesses were interviewed.These estimations were used to produce a maximum specific peak discharge map (Delrieu et al., 2005).In Table 1, the estimated peak discharges of 17 of these catchments, used for model verification, show that most of the tributaries of the Gard River have a peak specific discharge greater than 5 m 3 s −1 km −2 .In the first uphill of the mountains, in the region of Alès (Fig. 1), peak of specific discharges are even estimated over 20 m 3 s −1 km −2 .These specific discharges are among the highest ever reported for this range of watershed sizes (Gaume and Bouvier, 2004).
To specify the intensity of this event, it is important to notice that the 10 years return period discharge for such catchments is about 2 m 3 s −1 km −2 in this region (Delrieu et al., 2005).The 4 small ungauged catchments under study were not explicitly surveyed but Delrieu et al. (2005).Gaume and Bouvier (2004) present in more details the methodology of the flood peak estimation for several small nearby catchments.

Model presentation
The hydrological model, used in this study, is developed within the LIQUID hydrological modeling platform (Viallet et al., 2006).This platform allows the elaboration of independent process modules that can exchange variables and fluxes with other modules.They can be assembled to build "à la carte" models, according to the modeling objectives and the available data.Each module is run with its own characteristic time and space scales and the geometry of modeling units does not need being regular.A time sequencer, which is an event-based simulator, allows running the simulation and synchronizing the various modules.
To represent soil vertical heterogeneity and both infiltration excess and saturation excess mechanisms, a 1-D vertical soil water transfer module is used.It is based on the mixed form of the Richards Eq. (1).
In this equation, h is the water pressure head (m), θ is the soil water content (m 3 m −3 ), t is the time (s), z the depth (m) and K the hydraulic conductivity (m s −1 ).
Equation ( 1) is solved through the simplified Ross' algorithm (Ross, 2003), as validated by Varado et al. (2006).It requires the hydraulic characteristics of the soils through the hydraulic conductivity curve K(θ ) (Eq. 2) and the retention curve θ (h) (Eq. 3) based on the Brooks and Corey's model (1964).The hydraulic conductivity is expressed as: where K s and θ s are the hydraulic conductivity and water content at saturation (m s −1 and m 3 m −3 ); η defines the dimensionless shape parameter of the hydraulic conductivity curve.
The retention curve is expressed as follows: where h bc is the pressure head at air entry (m) and λ is the shape parameter of the water retention curve.Equation ( 4) is used in the following to express the relationship between the two shape parameters (Childs and Collis-George, 1950;Haverkamp et al., 1999).
The corresponding 1-D soil water transfer module is run for each hydro-landscape (defined as the elementary computing volume as described in Sect.3.2).It allows the computation of infiltration and soil water distribution within the soil profiles.Both saturated and unsaturated flow can be sim-ulated.Therefore, it is possible to simulate infiltration excess and ponding, perched water tables and/or the full saturation of soils profile.If a water table is present in the hydrolandscape, the model is also able to simulate its rising due to vertical infiltration.On the other hand, no lateral transfer is taken into account in the model at present stage.Thus the water table cannot rise due to lateral transfer.The flow routing in the river is provided by a 1-D kinematic wave approximation of the Barré de St-Venant equations (Moussa and Bocquillon, 1996) used with simplified trapezoidal cross sections.The river network is extracted from DEM analysis and discretized into river reaches, draining the sub-catchments.Each river reach has its own As a first guess and in order to avoid calibration on parameters that surely may vary from one catchment to another and within the catchment itself, the initial wetted section and roughness coefficient are fixed along the whole network with the constant values of 0.06 m 2 and 0.05, respectively.The value of the Manning coefficient has been provided by the post flood field survey (Gaume and Bouvier, 2004).
The cross section is reduced along the network according to the Strahler order of the reaches.The first Strahler reach is chosen to be rectangular and 1-m width.The average slope of each reach is computed using the DEM analysis proposed by Travis et al. (1975).

Hydraulic properties of soils
To run the 1-D soil transfer module, the retention and hydraulic conductivity curves are specified thanks to pedotransfer functions (PTF) which relate the parameters of the Brooks and Corey models K s , θ s , h bc , λ and η (Eqs. 2 and 3) to measured soil data (i.e.soil texture, organic matter content and/or other routinely measured data by soil surveys (Wösten et al., 1999).Measurements of soil infiltration are only available very locally and are therefore difficult to extend to the whole region.In the present work, we evaluate the relevance of two PTFs for the whole region of interest.The BDSol-LR is therefore used to document the soil properties.
The first class of PTF, based on the Cosby et al. (1984) approach (later refer as C84), uses only textural information.Its formulation is presented in Table 2.
Experimental studies (Ayral, 2005;Marchandise, 2007) highlighted the determinant role of soil structure in the Cévennes -Vivarais region flood genesis.To account for soil structure influence, the Rawls and Brakensiek (1985) (later refer as RB85) formulation is used as a second PTF class.Table 3 gives the RB85 formulation.In this case, the soil structure is taken into account through porosity.However, while the soil texture is explicitly given in the BDSol-LR database, the porosity needs being derived.Brakensiek et al. (1981) used a statistical approach to derive an average value of the effective porosity and the associated standard deviation for each textural class according to FAO (1990).The porosity is therefore obtained from the Brakensiek et al. (1981) results, which provide range value as shown in Fig. 4, and from qualitative information available in the BDSol-LR, such as "massive structure" or "gritty structure", which is used to choose between the minimum, maximum or average value within this range.
In addition, the observed rock fragment rate (not shown here but present in the BDSol-LR database) has a 20% average value in the region with some areas where it exceeds 50%.Based on these observations, the values of K s and θ s were reduced to account for the presence of rocks that can strongly affect the infiltration rate (Brakensiek and Rawls, 1994;Morvan et al., 2004).Simply, the values of K s and θ s estimated by the PTFs, are then reduced by the rock fraction (considered impervious and non porous).For instance, in areas where the rock fragment rate is 50%, the estimated values are divided by two.
The consistency of the hydraulic properties estimated using the two PTFs is evaluated by comparison with the available information about soil hydraulic properties (see Sect. 4).

Modeling setup
The model is implemented for the whole Cévennes -Vivarais region and run on the small catchments without any parameter calibration.Within the LIQUID framework, the  Rawls and Brakensiek (1985).The domain of validity is defined as 5%< sand <70% and 5%< clay <60%.model setup requires choosing the simulated processes, their representation and consistent catchments discretization.The followed methodology is the one outlined by Dehotin and Braud (2008) who propose a catchment discretization consistent with the modeling objectives, the process representation and the available data.The discretization (Fig. 5) is based on two embedded levels.
-For the first level, the catchments are sub-divided into sub-catchments or Representative Elementary Watersheds (REWs).They are extracted from the stream network previously derived from a DEM analysis (Digital Elevation Model).The latter was conducted using the SAGA 5 software, using 0.1 km 2 as the surface threshold and the first Strahler order.This first level of discretization is not sufficient to describe the soil heterogeneity within the sub-catchments.A second level of discretization is therefore required.
-The variability inside the sub-catchments can be accounted for through homogeneous zones, called hydrolandscapes, which are derived from GIS layers analysis.
In this study, this second level is based on the Cartographic Soil Units (CSU) presented in Sect. 2. Each 5 http://sourceforge.net/projects/saga-gis/CSU contains several Typological Soil Units (TSU) that are not georeferenced.In the present study, we affect the dominant TSU to each CSU for the sake of simplicity.The soil depth is fixed to the average depth for each TSU.To solve the Richards equation, the vertical soil profile is divided into cells of 1 cm thickness.A different set of hydraulic parameters is assigned to each soil horizon composing the soil profiles.
Because the September 2002 event is the first main rainy event after the summer dry season, the model is initialized as follows.We first consider a uniform soil water pressure profile.The value of this pressure corresponds to a 75% saturation of the first horizon.A simulation without rainfall but with a constant potential evaporation of 2 mm day −1 is then run for two months.The obtained soil moisture profiles are used as initial condition for the simulation of the September 2002 event.
As far as the boundary conditions are concerned, a zero flux (impermeable bedrock) condition is fixed at the bottom boundary.Preliminary tests showed that the use of a gravitational flux as boundary conditions could modify considerably the soils response.This result shows that field studies would be required to better document this boundary and possible percolation within the fractured bedrock.In the following, we will only present results with the no flux bottom boundary condition.The upper boundary condition of the soil is directly calculated by the model using input rainfall series and a zero potential evapotranspiration (evaporation is neglected during the event).Rainfall data at the 1 km 2 resolution are interpolated per hydro-landscapes using weighted averages on the surfaces.The module time step is automatically calculated as a function of the dynamics (Ross, 2003) and the management of all the hydro-landscapes units is ensured by the time sequencer of the LIQUID platform (Viallet et al., 2006).
As mentioned above, ponding water is directly transmitted to the nearest river reach and transmitted within the river network using the kinematic wave routing module.It provides the simulated streamflow at the outlet of the studied catchments and for their sub-catchments.The model also provides, for each hydro-landscape, the infiltration flux and the vertical fluxes between the cells.The saturation degree profiles are also available and are used in the analysis provided in the next section.

Model verification strategy
Since very few discharge data are in general available during flash flood, the use of post flood field survey is very interesting and provides a spatial view of the hydrological response across a large range of scales (Borga et al., 2008).Of course, only information on peak discharges can be retrieved and information on the event timing can be available when witnesses are interviewed.Our hypothesis is that a regional evaluation using peak discharge is a good starting point for the evaluation of our model.It is also a good complement to compare the simulation within gauged catchments.Nevertheless, concerning the rainfall inputs, the availability of radar data is a requisite for a fair evaluation of models at scales that are not properly covered by raingauge data.Given the current state of radar data, it is quite difficult to have a long reliable series of radar data at a given location.This point is especially critical when we need high resolution (space and time) rain inputs for physically based hydrological models where the simulated process has a time step lower  1984) and (b) Rawls and Brakensiek (1985) formulations.These results are obtained after a two months evaporation of initially wet surface soils.
than the hourly time step.This point is also valid for lumped models that need accurate rain inputs as well.The model proposed in our study has been evaluated (Anquetin et al., 2009) on one gauged catchment (Saumane, 99 km 2 ) using the same numerical implementation.The results show very good performances since the radar protocol is adapted to mountainous areas (volume-scanning protocol).
The model verification strategy is, therefore, based on a regional modelling approach based on the simulation of the discharges of 17 watersheds documented during the post flood survey (Gaume and Bouvier, 2004).The size of the catchments and the estimated discharges are given in Table 1.Their locations are indicated in Fig. 1b.
The numerical configurations for each studied catchments are summerized in Table 4.The numbers of sub-catchments and hydro-landscapes used for the simulations, give an idea of the heterogeneity taken into account.

Soils infiltration properties
To illustrate the impact of the two different PTFs on modeled soil hydraulic properties, simple statistical scores on the whole database of the region of interest are summarized in Table 5.The compared values clearly show the impact of the PTF on the estimation of the Brooks and Corey parameters.The variability range is much larger for the RB85 formulation, especially for the λ, K s and h bc parameters.In Table 5, the two PTFs lead to very close maximal water storage capacity, defined as the product of the depth of the soil and the water content at saturation θ s .Figure 6 shows however that the obtained initial available soil water capacity calculated for the September 2002 event presents large differences between the C84 and the RB85 formulations.These are of 35 mm in average and can reach 535 mm for very deep soils (Table 6).The average initial available soil water capacities are consistent with the estimation cited in Delrieu et al. (2005).
The evaluation of K s is difficult as only few point measurements on the Gardon d'Anduze catchment (545 km 2 ) are available.On this area, Ayral et al. (2005) and Marchandise (2007) reported values larger than 100 mm h −1 .For the same area, the RB85 estimates saturated hydraulic conductivities ranging from 1 to 90 mm h −1 , which are lower than these observations.However, in situ data take into account local macroporosity whereas it is not included in the derivation of PTFs.This result highlights the need for further work in order to include macroporosity in the definition of PTFs.

Fig. 7.
Comparison between the estimated (empty squares) and the simulated (black dots) specific peak discharges at the 17 validation outlets.The vertical black line stands for the range of uncertainties given by Gaume and Bouvier (2004).
For the present study, the use of PTF is the only way to get hydraulic properties for the whole region, in a homogeneous way.In the following, we use the RB85 PTF as it provides a more realistic range of variation than C84.Nevertheless, it would be necessary to complement the existing field study in order to document the variety of soils encountered at the regional scale.This could allow establishing a PTF for the region.

Model verification
Figure 7 displays the comparison between the estimated and the simulated peak discharges at the 17 validation outlets.As shown in Fig. 7, most of the simulated peak discharges are within the range of the uncertainties given by the post flood survey.Apart for catchments #7 and #10, the relative error between the simulated and the estimated peak discharges ranges from 3% (#12) to 47% (#15).The mean error value is about 22% whereas the estimation uncertainty provided by the post flood survey reaches 21%.This is a fairly good result for this type of regional validation without any calibration.The simulation of the peak discharge at the outlet #6 fails.As mentioned by Bonnifait et al., 2009, this catchment is located in a region where the rainfall gradient has been estimated to be the strongest.The rainfall input is, therefore, probably at the origin of the overestimation of the peak discharge.The underestimations of the simulated peak discharge for the catchment #10 is difficult to explain since there is no modeling, background on this catchment.This catchment will be a good candidate for a special observation device within the framework of the future HyMEx experiment in order to better document the soil characteristics.

Flood dynamics for the 4 catchments
Figure 8 displays the observed rain field and the simulated specific discharge at the 4 outlets (Rousson, St. Quentin la Poterie, Quissac and Domazan).Due to the size of the catchments and the importance of the event, the hydrographs approximately follow the rain dynamics.In Fig. 9, the simulated peaks for the 4 studied catchments and the 17 validation outlets are compared against the values retrieved from the post-flood field investigation (Gaume and Bouvier, 2004), used as validation data.For the 4 studied catchments, the specific discharges reach 15 to 25 m 3 s −1 km −2 .The compared values are of the same order of magnitude and the hierarchy between catchments is correctly simulated by the model, which was not obvious from the mere examination of the rain amounts collected by the watersheds.
In order to illustrate the influence of soil characteristics on the water available amount for runoff, Fig. 10 displays the time evolution of i) the average infiltration rate for each catchment, ii) the mean rain intensity and iii) the corresponding cumulative runoff coefficient.This latter is computed as the simulated cumulative water amount available for runoff (ponding) divided by the cumulated rainfall up to this time.Runoff generation is obviously strongly linked to the rain intensity but also to the soil storage capacity and the infiltration capacity.In Fig. 10, the soils of the Rousson and St. Quentin la Poterie catchments are fully saturated at 06:00 UTC September the 9th.On the other hand, the soils    (Gaume and Bouvier, 2004).Simulated peak discharges using LIQUID hydrological model for the 17 validation catchments (grey triangles) and the 4 studied catchments (black triangles).   of the Domazan and Quissac catchments do not reach their saturation state as the infiltration rate remains positive during the whole event.
These results highlight the combined contribution of the rain intensity and the soil properties in the runoff generation at the catchment scale.The next section aims at identifying how the soil properties act on the infiltration rate and if any main hydrological process can be identified.

Infiltration dynamics associated to soil properties
In Fig. 11, the average simulated runoff coefficients are evaluated for the different sub-catchments and plotted against the cumulated rainfall during the whole simulation (08:00 UTC the 8th of December to 00:00 UTC the 10th of December).
Despite the usual strong correlation between runoff coefficient and rainfall intensity, some discrepancies are observed that are assumed to be linked to the soil signature.For example, the largest runoff coefficient for the Rousson catchment occurs where the cumulated rainfall is the lowest.For the Quissac and more specifically, the St.-Quentin la Poterie catchments, the total rainfall amount is almost spatially homogeneous but the calculated runoff coefficients can be seen highly variable amongst sub-catchments.
To illustrate the role of soil properties on the runoff genesis, the average infiltration rates are plotted for each subcatchment in Fig. 12 for the Domazan catchment.This figure, combined with the analysis of ponding (not shown), highlights two different tendencies for the given example.The west/southern (W/S) sub-catchments of Domazan present a globally higher infiltration capacity than the east/northern (E/N) ones.This is particularly emphasized during the quasi no-rain phase around 11:00-13:00 UTC (see Fig. 10d).The infiltration process is rapidly stopped within the W/S sub-catchments as no ponding is simulated.On the contrary, the infiltration still goes on within the E/N subcatchments as ponding remains at the surface.At 14:00 UTC, the rain rate has increased.The higher infiltration capacity of the W/S sub-catchments allows infiltration until 22:00 UTC on September the 8th.At this time, the sub-catchments are fully saturated.These sub-catchments probably react as saturated source areas and present lower average runoff coefficients than the E/N ones.
For the other catchments (not presented here), similar tendencies are observed as far as the infiltration dynamics is concerned.
In Fig. 13, the vertical profiles of the soil saturation state are presented at different instants of the event and for two hydro-landscapes representative of the two tendencies.HL1 (Fig. 13a, b, c, d) refers to the W/S sub-catchments of Domazan (larger infiltration capacity) while HL2 (Fig. 13e, f,  g, h) stands for the E/N sub-catchments that have lower infiltration rates.These two hydro-landscapes receive approximately the same average rainfall.Therefore, Fig. 13 shows different behaviors that can be related to soil properties.).The black line stands for the mean evolution of the simulated rates in grey.For each sub-catchment, the average specific discharge (m 3 s −1 km −2 ) (in bold) and the maximum discharge (m 3 s −1 ) (in parenthesis) are indicated.
The time evolution of the infiltration in HL1 (Fig. 13a,  b, c, d) presents a strong correlation with the rain dynamics.At the beginning of the event (Fig. 13a) the first layers of the soil saturate and when rain stops (Fig. 13b), the water infiltrates deeper while the upper layers start drying up.During this phase, there is no ponding at the HL1 surface.When rain starts again (Fig. 13c), the surface saturates and ponding starts increasing.The infiltration front moves on quickly towards the deeper layers; that leads to the complete soil column saturation at 21:47 UTC the 8th of December.Then, runoff generation follows the rainfall dynamics.So, this hydro-landscape experiences both infiltration excess and saturation excess mechanisms.
In the case of HL2, ponding starts immediately (Fig. 13e) and the surface remains saturated (Fig. 13f), contrary to HL1 (Fig. 13b).The infiltration front moves on slowly towards the deeper layers.The first layer of the soil saturates the 9th of September at 10:05 UTC; the second remains very moist, but does not saturate during the event.Due to the low infiltration capacity of the soil, Hortonian infiltration excess is Fig. 13.Time evolution of the vertical profile of the saturation state for the two hydro-landscapes of Domazan recalled in the map.Hydrolandscape 1, referred as HL1 (Fig. 12a, b, c, d), is representative of the W/S sub-catchments that present important infiltration capacity.Hydro-landscape 2, referred as HL2 (Fig. 12e, f, g, h), stands for the E/N sub-catchments that have lower infiltration rates.For both hydrolandscapes, the time run from left to right: the 8 September 2002 1) from 08:50 UTC to 09:40 UTC, 2) from 09:40 UTC to 15:00 UTC, 3) from 15:00 UTC to 21:47 UT, 4) from 21:47 UTC to 9 September 2002 16:00 UTC.
the only active hydrological process on this hydro-landscape.On other soils, we have found soils prone to saturation excess due to a lower infiltration capacity of the deeper layer than the upper one, leading to a perched water table.Given the initial conditions (without pre-existing water table) and the absence of lateral transfer in the present version of the model, we are only able to evidence saturation excess linked to saturation of the topsoil or complete saturation of the soil reservoir (referred to as type-B saturation excess by Latron and Gallart, 2007).Saturation excess due to groundwater rising (type-A saturation excess of Latron and Gallart, 2007) could be simulated in case of a soil with very high hydraulic conductivity for which the infiltration front would reach quickly the column bottom and subsequent infiltration would produce a rising of the water table.
The same analysis procedure is performed for the three other catchments.Figure 14 synthesizes the results as a map distinguishing the areas where the soils fully saturate during this event and the areas where it does not.Note that during such an extreme event, rainfall intensities are so high that infiltration excess occurs over all the soils, but it seems to be preponderant on the white areas of Fig. 14.This map can serve as basis for an experimental set up, aiming at verifying these hypotheses.For instance, areas suspected to be prone to saturation excess could be equipped with small piezometers, whereas the other areas could be instrumented with surface sensors, measuring soil moisture.

Conclusions
This study is a first step towards the set up of a modeling approach for ungauged catchments in a region prone to flashfloods events.A simple model is established using the facilities provided by the LIQUID modeling platform.The soil parameters required for the simulation are derived from the available database BDSol-LR.The model was verified using a regional approach based on estimated peak discharges, provided by a post flood survey.Then, the analysis focuses on the simulated hydrological response of four smallungauged catchments to the variability of soil properties, for the September 2002 catastrophic event.The results highlight the importance of soil hydraulic properties on the hydrological response.They also show, even on small catchments of less than 10 km 2 , the possible co-existence of different runoff mechanisms (infiltration excess and saturation excess).The simulated maximum specific peak discharges are consistent with the post-event investigation for the four sub-catchments.These results must be confirmed at larger scales where observed simulated discharges are available.This work is under way through the extension of the present analysis to larger catchments.
This first study also highlights a lack of knowledge in various domains, which provides guidance for the establishment of field study priorities in the future HyMeX experiment.One can mention: -The necessity to validate and/or calibrate pedotransfer function for characteristic soils of the region and to document soil properties, including structure; -The results are related to the definition of the bottom boundary condition and would require investigations about infiltration capacity of bedrock layers, which are generally considered as impermeable; -The extension of the modeling study presented here will allow mapping areas principally prone to infiltration excess and to saturation excess.It would be interesting to validate, even with qualitative information, these hypotheses; -The validation of a regional modeling approach would also require to enhance stream flow estimation on small upstream catchments, using embedded scales.Image video (Muste et al., 2005) should be encouraged in such a context.Previous sensitivity studies have shown the strong impact of the initial soil moisture (which specifies the initial storage capacity of the soil).Continuous simulations, including evapotranspiration simulation will be included in the modeling approach to assess this impact.This would therefore require to consider lateral redistribution of water along the slopes, while adding a lateral flow component to the present model.

CFig. 2 .
Figure 2 Fig. 2. Data extracted from the BDSol-LR database: (a) Average soil depth (in cm).(b) Variability of the texture of the soils in the Languedoc -Roussillon.The texture classes are taken from the FAO classification.C stands for clay, SC sandy clay, SCL sandy clay loam, CL clay loam, L loam, Si silt, Sic silty clay, SiCL silty clay loam, SiL silt loam, LS loamy sand, SL sandy loam and S stands for sand.

Figure 4 Fig. 4 .
Figure 4 Fig. 4. Statistics of the porosity for different textures (after Brakensiek and Rawls, 1981).The bolded values stand for the mean value of the class whereas the light ones stand for the range within the class.

Fig. 9 .
Fig. 9. Maximum specific discharge as a function of catchment size.Peak specific discharges (empty and grey squares) estimated from post-flood investigations of the September 8-9th event(Gaume and Bouvier, 2004).Simulated peak discharges using LIQUID hydrological model for the 17 validation catchments (grey triangles) and the 4 studied catchments (black triangles).

Table 1 .
(Gaume and Bouvier, 2004) catchments used for the model evaluation(Gaume and Bouvier, 2004).The most probable values are in bold whereas the values in brackets give the range of uncertainty.

Table 4 .
Numerical configurations of the studied catchments.

Table 6 .
Rawls and Brakensiek (1985)nd initially available water storage capacities (θ s depth and (θ s −θ i ) depth) (m) estimated by theCosby et al. (1984)(refer as C84) andRawls and Brakensiek (1985)(refer as RB85) pedotransfer functions.VC is the Coefficient of Variation calculated as the Standard Deviation divided by the Average.