Flood modelling with a distributed event-based parsimonious rainfall-runoff model: case of the karstic Lez river catchment

. Rainfall-runoff models are crucial tools for the statistical prediction of ﬂash ﬂoods and real-time forecasting. This paper focuses on a karstic basin in the South of France and proposes a distributed parsimonious event-based rainfall-runoff model, coherent with the poor knowledge of both evaporative and underground ﬂuxes. The model combines a SCS runoff model and a Lag and Route routing model for each cell of a regular grid mesh. The efﬁciency of the model is discussed not only to satisfactorily simulate ﬂoods but also to get powerful relationships between the initial condition of the model and various predictors of the initial wetness state of the basin, such as the base ﬂow, the Hu2 index from the Meteo-France SIM model and the piezometric levels of the aquifer. The advantage of using meteorological radar rainfall in ﬂood modelling is also assessed. Model calibration proved to be satisfactory by using an hourly time step with Nash criterion values, ranging between 0.66 and 0.94


Introduction
Rainfall-runoff models are crucial tools for the statistical prediction of flash floods and real-time forecasting. However, their efficiency is still limited by uncertainties related to both the spatial variability of rainstorms (Sangati et al., 2009) and the assessment of the initial wetness state in the hydrosystems (Brocca et al., 2008). In karstic basins, additional uncertainties concern the geometry and the hydrodynamics of the aquifer, and more generally the impact of the aquifers on the flood processes (Jourde et al., 2007;Bailly-Compte et al., 2008a).
Because of these uncertainties and poor knowledge of hydrological fluxes, which is often limited to rainfall and discharge at the outlet, rainfall-runoff models generally have to be calibrated beforehand, particularly in karstic catchments. In this case, models that require few hydrological data and no detailed description of the structure of the karstic catchment are therefore attractive. Black-box or reservoir models are therefore preferred to physical-based models and have thus been widely used for karstic catchments. Some authors use inverse modeling to determine rapid and slow flow impulse response of the karstic system and either identify some properties of karstic systems (Pinault et al., 2001) or propose groundwater level warning thresholds for flash floods (Maréchal et al., 2008). Parsimonious reservoir rainfall-runoff models with interconnected reservoirs are also used as management tools to predict spring discharges (Fleury et al., 2007(Fleury et al., , 2009 or regional water levels (Barrett and Charbeneau, 1997). Reservoir rainfall-runoff models are also used to represent the discharge of a river that crosses a karstic area. Rimmer and Salingar (2006) used HYMKE, a "gray box" model, to simulate the discharge of the Jordan river and its three tributaries. Le Moine et al. (2008) used a continuous Genie Rural (GR) model in the French catchment of Larochefoucault to simulate Tardoire and Bonnieure discharge and losses. Finally, Bailly-Comte Published by Copernicus Publications on behalf of the European Geosciences Union. et al. (2008b) uses an event-based reservoir semi-distributed rainfall-runoff model that reproduces karst-river interactions at a local scale to simulate Coulazou (Mosson tributary, near Montpellier, South of France) flash floods.
A critical point in flood modelling is the assessment of the wetness state of the catchment at the beginning of intense events. Continuous rainfall-runoff models are able to assess the soil moisture at the beginning of a rainfall event. Parsimonious models, like the French Génie-Rural (GR) family of models (Perrin et al., 2003), are very popular because they only require the assessment of the potential evapotranspiration, along with rainfall, to calculate soil moisture. However, continuous models can have certain disadvantages: (1) they require complete time-series of rainfall, discharge and meteorological data; (2) the wetness state of the catchment results from several combined processes such as evapotranspiration, deep aquifer percolation, or hillslope flow, which could not be described with only a few parameters. When considering only floods, one strategy consists in not taking into account the inter-flood period and using event-based models. However, in this case, the initial wetness state of the catchment has to be assessed from external information such as baseflow (Franchini et al., 1996;Fourmigué and Lavabre, 2005), surface model output (Goodrich et al., 1994;, remote sensing data (Quesney et al., 2000;Pellarin et al., 2006), or in situ measurement of soil water content (Brocca et al., 2008;Tramblay et al., 2010). Although continuous models could be found to be more efficient than event-based models (Berthet et al., 2009), event-based models should be very convenient for operational purposes, if the initial wetness state of the catchment would be known with good accuracy.
In addition, distributed models can greatly improve simulations, because they account for spatial variations in the climatologic or geographic inputs of the model (Arnaud et al., 2002). However, distribution can considerably increase the number of the parameters of the model making calibration difficult (Chaponnière, 2005). But using distributed input, e.g. radar rainfall measurements, is a good way to take advantage of the spatial information, without requiring additional parameters. In this sense, distributed models can still be parsimonious.
The concept of distributed parsimonious event-based rainfall-runoff models thus appears to be a good way of optimizing the model's efficiency in catchments where hydrological fluxes and features are poorly known. Such models seem to be rare for flash flood simulation in karstic catchments: it is used by Bailly-Comte et al. (2008b) to simulate the Coulazou karstic catchement flash floods. Thus, the main objectives of this paper are first to test an appropriate distributed event-based parsimonious model in a case study on a karstic basin, and second to look for relevant predictors to derive the initial condition of the event-based mode. An additional objective of the paper is to show how radar rainfall measurements can improve both flood simulation and param-eters calibration, focusing on the latter for the calibration of the initial condition and corresponding impact on the relationship with predictors. First we provide an overview of the model, and then describe the case study area and the available data. Model calibration is described with different rainfall inputs. The initial condition of the model is finally related to predictors like surface model output, baseflow and piezometric level.
2 Structure of the model The hydrological model operates over a regular grid mesh of cells, which are not connected. Rainfall from gauge stations is computed for each cell at any time, according to the method of the Thiessen polygons. The runoff from each cell is then calculated using a modified SCS runoff model, accounting for rainfall intermittency and subsurface runoff. Each cell produces an elementary hydrograph, routed at the outlet using a Lag and Route model. The complete hydrograph of the flood is finally obtained after addition of the elementary hydrographs.

Runoff model
The effective rainfall (or runoff) intensity i e (t) at a given time t is derived from the SCS instantaneous formulation (Gaume et al., 2004) where i b (t) denotes the precipitation intensity at the time t, P (t) the amount of rainfall since the beginning of the event, and S the water deficit at the beginning of the event. Please note that the runoff coefficient at a given time t only depends on both S (which does not change during the whole event) and the cumulated rainfall P (t) since the beginning of the event. P (t) can be considered as the level of the cumulated rainfall reservoir, whose capacity is infinite (Fig. 1). The S (L) value is assumed to be dependent of the landuse, the type of soil and the initial moisture conditions (Mishra and Sigh, 2003). For this last reason, the S value is expected to be event-dependent. It thus differs from the more general SMA-SCS procedure proposed by Michel et al. (2005), which added a non-zero initial level in the soil reservoir and considered S as a constant parameter for a given catchment.
In order to account that the potential runoff coefficient can decrease during periods when no rain occurs (rain intermittency), the cumulated rainfall reservoir was drained by a discharge which depends linearly on the level in the reservoir.
Nat. Hazards Earth Syst. Sci., 12, 1119-1133, 2012 www.nat-hazards-earth-syst-sci.net/12/1119/2012/ Thus the level P (t) in the cumulated rainfall reservoir is actually calculated by where the discharge coefficient Ds (T −1 ) is assumed to be constant for a given catchment. This allows that the runoff coefficient decreases during periods with no rain, because of near surface evaporation or deep drainage of the soils for example. The level of the reservoir at the beginning of each event is P (0) = 0. In order to account for slow discharge of soils or aquifers, an additional runoff was considered, as a part of the cumulated infiltration since the beginning of the event. The runoff coefficient divides the rainfall in runoff and infiltration. The infiltration then fills a "cumulated infiltration reservoir", whose finite capacity equals S (Steenhuis et al., 1995) (Fig. 1). The level stoc(t) of this reservoir is computed by where Ds, the coefficient of discharge of the reservoir, due to losses by evaporation, percolation to the deep aquifer and lateral flow. Both "cumulated infiltration reservoir" and "cumulated rainfall reservoir" have the same Ds coefficient of discharge because the two reservoirs must be empty at the same time. The SCS formulated with both discharge and delayed runoff indeed relies on duality between the cumulated rainfall reservoir and the infiltration reservoir. The runoff coefficient is directly related to the cumulated rainfall, but it is expected that it is also, from a physical point of view, related to the soil reservoir. This means that no runoff should correspond only to a cumulated rainfall equal to 0, but also to a cumulated infiltration equal to 0. In addition, both levels must be linked in the same way all along the event. In other terms, a given runoff coefficient value (i.e. a cumulated rainfall value) must correspond to a given saturation level (i.e. a cumulated infiltration) in the model. Thus, the model deals with the competitive relation between both saturated and unsaturated storage, as mentionned in Penna et al. (2011), since S expresses the unsaturated storage (as initial water deficit), but depends on the saturated storage, through the relationship with SIM or piezometric level. The same Ds coefficient ensures the two above conditions. The level of the reservoir at the beginning of the event is stoc(0) = 0. The additional delayed runoff i d (t) is finally expressed as a proportion, min(1,w/S) of the discharge of the "cumulated infiltration reservoir" which corresponds to the part of the cumulated infiltration which is, however, routed to the outlet of the catchment (soil water or groundwater) by lateral flow. The parameter w was added to quantify the part of the cumulated infiltration reservoir to simulate the delayed runoff of the observed hydrograph. The w value (L) is assumed to be constant for a given catchment and is used to fit the delayed runoff of the simulated hydrograph to the corresponding observed discharges. Because of the seasonal variations of the observed delayed runoff depending on the initial wetness state, min(1,w/S) (ranging between 0 and 1) was used to express the proportion of the reservoir discharge to be converted into delayed runoff. Given that S is constant for an event, but varies from one event to another, w/S also varies from one event to another. w/S increases when S decreases (i.e. soils are initially wet), and w/S decreases when S increases (i.e. soils are initially dry). In other terms, the contribution due to the drainage of the soils increases with the initial soil saturation. The total runoff i t (t) at the time t is thus given by: The runoff model finally deals with three parameters: S, the initial water deficit, Ds, the loss rate of both the "cumulated infiltration reservoir" and the "cumulated rainfall reservoir", and w, the part of the reservoir loss which emulates the delayed runoff (Fig. 1). The Eqs.
(3) and (4) are integrated at a given time step using an explicit centred scheme. This kind of model is able to suit a large range of flood processes. The water deficit can indeed have several meanings, referring for example to the surface moisture conditions for hortonian processes, the soil water content for dunnian processes or even the state of the deep aquifers in case of underground and karstic processes.

Routing model
The total runoff is then routed from each cell m to the outlet by a linear lag and route model (Maidment, 1992 (26-9); Bouvier and Delclaux, 1996;Lhomme et al., 2004). Each total runoff at time t 0 from each cell m generates an elementary hydrograph Q m (t) at the outlet of the catchment ( Fig. 2) where T m denotes the travel time from the cell m to the outlet (route), K m the diffusion of the initial input along the travel paths (lag), A the area of the cell. The Eq. (7) is integrated at a given time step using an explicit centred scheme.
The travel time T m is calculated from where l k are the lengths of the k-cells between the cell m and the outlet (Fig. 2), V 0 a velocity of travel, which is supposed to be uniform over the whole catchment (route parameter). The diffusion time K m is assumed to be proportional to the travel time T m , using a constant coefficient K 0 which is also supposed to be uniform over the whole catchment (lag parameter): The complete model finally deals with five parameters: S,w, Ds,V 0 ,K 0 , which remain the same for all the flood events, except for the initial condition S. Model parameters are also uniform in space. This model was implemented in the ATHYS modelling platform (www.athys-soft.org).

Lez catchment
The study site is the 114 km 2 Lez catchment at Lavalette, upstream from Montpellier (Fig. 3). Aquifer geology is dominated by Dogger, Malm and Berriasian limestones and dolomites with a thickness ranging from 650 to 1100 m. The hydrogeological catchment is supposed to extend about 380 km 2 (Avias, 1992) notably at the north of the topographic basin. Relief alternates between limestones plateaus named "causses" and plains. In the plains, limestones are covered by 200 to 800 m thick Valanginian marls, and soil whose thickness is generally less than 1 m. Altitudes range from 300 to 700 m in the calcareous causses and from 50 to 100 m in the marly plains. Vegetation is primarily scrub on the causses and crops (vineyard, olive trees) on the plains.
Floods occur mainly in autumn during intense rainstorms that can reach several hundred millimetres in 24 h (more than 300 mm at the Matelles in 1976). Roesch and Jourde (2006) showed that water storage in the hydrosystem (soil and karst) is much more important during the first autumn floods following the dry season than during floods at the end of autumn.

Rainfall and runoff data
Rainfall data were recorded at a one-hour time step and came from four recording rain gauges (Prades, St. Martin, Ensam, Mauguio) monitored by Meteo-France of which only the Prades station is located within the Lez basin, the other three being located at a distance of 5 to 10 km outside the catchment (Fig. 3). The mean annual rainfall at Prades is 909 mm (period 1994-2008). Rainfall data from the Nîmes radar (Calamar or Hydram treatment) were also available at a 5-min time step and a resolution of 1 km 2 . The main operations of these treatments consist on (i) a deletion of ground clutters and masks effects, (ii) an estimation of the vertical profile of reflectivity and (iii) a reflectivity conversion into a rainfall depth using the Marshall-Palmer Z-R relationship. Both Calamar and Hydram rainfall data are calibrated from the raingauges. Calamar calibration is uniform in sub-areas of the whole radar area, whereas it varies within the time: the correction factor is the mean error between the observations at each raingauge of a sub-area and the corresponding radar pixel, since the beginning of the rainy event; the date of beginning is updated after more than 3 time steps (of 5 min) without rain, or when the correction factor at the ith time step exceeds more than 15 % the (i-1) th one (Ayral et al., 2005). Hydram calibration is uniform over the whole radar area, and constant in time at monthly time step: the correction factor is computed as the difference between the spatially-averaged rainfall depths over the whole area since the beginning of the month; rainfall depth from raingauges derives from kriging interpolation (Cheze and Helocco, 1999).
The nearest available evapotranspiration data comes from Mauguio, at a daily time step; the mean annual potential evapotranspiration at Mauguio is 1322 mm (period 1996-2005), using the Penman-Monteith formula. Discharge time series at the basin outlet Lavalette were available at an hourly time step. The mean annual discharge at Lavalette is 2.1 m 3 s −1 (period 1987-2008).
The selected episodes correspond to twenty-one floods measured from 1994 to 2008 (Table 1). This sample covers a wide range of peak flow from 40 m 3 s −1 to 480 m 3 s −1 . Among these floods, ten of them present a peak flow of about 110 m 3 s −1 with a corresponding return period of two years, and four of them have a peak flow greater than 200 m 3 s −1 with a corresponding return period equal to or greater than five years (return periods have been derived from a Gumbel distribution fitted over the 1975-2011 period, see http://www.hydro.eaufrance.fr, Lez river). Five events occurred after a long dry period (September 2000, 2002, 2005and October 2008) and the sixteen others occurred after the first rains. Flood events have been initiated at 06:00 UTC to be coherent with the daily Hu2 index (see Sect. 3.4), which is computed at this time.
Lag-times (considered as the differences in time between centroids of both rainfall and runoff) are short: about 5 h with a standard deviation of 3 h for peak discharges greater than 50 m 3 s −1 . The runoff coefficients were computed as the ratio of the runoff volume Vr and the averaged total precipitation over the catchment; Vr was calculated as the difference between the total runoff and the base flow runoff; the baseflow discharge was set constant, equal to the minimum value before the rising of the flood). The runoff coefficients were about 0.4 after a long dry period (Table 1) when the aquifer showed low piezometric levels whereas they were about 0.8 after the first rains (Table 1) when the aquifer showed high piezometric levels. These runoff coefficients can be greater than 1 when the Lez spring overflows (i.e. piezometric levels greater than 65 m a.s.l.) and the aquifer piezometry is high at the beginning of the event (Table 1). This suggests that groundwater can actually contribute to flood volume in this latter case, mainly because of the extended hydrogeological area which has not been taken into account in the runoff coefficient computation (Table 1 main characteristics of the selected events).
To assess radar rainfall quality, linear regressions were established between the rain gauge measurements and the corresponding radar pixels. The corresponding radar pixel was derived from the mean value of the central pixel and its eight nearest neighbours, to normalize any local noise present in the central pixel. To check the spatial variability of radar rainfall, a first regression was established between the total storm depth of the 20 daily rain gauges and the 20 corresponding radar pixels in the 5000 km 2 area covered by the radar (Fig. 4). Then, a determination coefficient R 2 e was calculated.
To check the temporal variability of radar rainfall, the hourly rain gauge and corresponding radar pixel time series were correlated. This correlation was performed at each functioning continuous recording station (Prades, St. Martin, Ensam, Mauguio). A mean determination coefficient R 2 h was calculated for each event. Thus, R 2 e and R 2 h assessed the closeness of rainfall radar data to ground rainfall "reference" data from rain gauges both in space and time.
For each event, to assess the systematic residual error in radar rainfall estimations, the Mean Field Bias (MFB) coefficient (Vieux and Bedient, 2004), was calculated MFB = where n is the number of rain gauges taken into account in the calculation (here the 20 Meteo-France daily rain gauges); Gi is the rainfall recorded by the rain gauge at point i; Ri is the rainfall given by the radar at point i, both for the studied event (Table 2).
MFB values and determination coefficients R 2 e and R 2 h of the linear regression between rain gauges and radar depth Hydram and Calamar were calculated for each rainfall radar event.
Regardless of the method treatment on radar data, it was observed that the MFB is always greater than 1, i.e. that radar systematically underestimates rainfall. This error could come from the Z-R relation that is used and the treatments based only on radar data and applied to the reflectivity measurements (Borga, 2002). Except in 2002, we also noted that both determination coefficients R 2 e and R 2 h of the linear regressions are greater than 0.7 for storms at the beginning of autumn (September and October). These storms present high mean rainfall intensities (>7 mm h −1 ) that can indicate a rainfall event mainly convective. For the other storms (November, December and January), at least determination coefficient R 2 e is less than 0.5. These storms present lower mean rainfall intensities (<6 mm h −1 ) that can indicate a rainfall event with a stratiform component. Radar rainfall quality is better for convective rainfall events which occur at the beginning of autumn and is worse for more stratiform rainfall events which occur mainly at the end of autumn. Syst. Sci., 12, 1119-1133, 2012 www.nat-hazards-earth-syst-sci.net/12/1119/2012/ During this period, the poor quality of radar measurements is likely affected by both the weak vertical extension of the precipitation system and the low altitude of the 0 • C isotherm (Bourrell et al., 1994), because of the distance between the basin and radar (∼60 km). For each event, and for each pixel, radar rainfall intensities were corrected by the MFB coefficient.

Piezometric data
A set of 12 piezometers connected to the Lez spring (Karam, 1989) was retained as reference gauges of the aquifer level. These piezometers are evenly distributed throughout the Lez aquifer (Fig. 3). A group of seven piezometers covers the southern part of the aquifer, within the topographic basin. The five other piezometers are located in the northern part of the aquifer, outside the topographic basin. The correlation matrix calculated using piezometric data at the beginning of the events between 2000 and 2008 (Table 3) shows that all piezometers are strongly correlated, except the piezometer 10. This bad correlation can be explained by the important distance between this piezometer and the Lez spring. Some piezometers (1, 3, 6, 11) show small seasonal groundwater level variations (about 30 m) and are greatly influenced by the Lez spring overflow. The other piezometers show greater groundwater level variations (from 40 to 80 m) and are less influenced by the Lez spring overflow.

Hu2 soil moisture index
The Safran-Isba-Modcou (SIM) model was developed by Météo-France and enables the soil wetness index to be com-puted for the whole of France (Habets et al., 2008;Quintana Seguí et al., 2009). The ISBA component of the model (Noilhan and Mahfouf, 1996) combines elevation, land cover and soil characteristics with atmospheric input to estimate, among other variables, soil moisture conditions. The percentage of soil saturation is available daily at 06:00 UTC for cells of 8 × 8 km 2 at three different levels in the soils (Hu1, Hu2, Hu3 for surface horizon, root horizon, and deep horizon, respectively) Hu = θ θ s 100 (11) where θ denotes the soil moisture and θ S the saturated volumetric soil moisture. The value of the root layer Hu2 was selected here because it is a priori the most representative moisture index of the superficial deposits (Marchandise and Viel, 2008). The values used for Hu2 correspond to the average moisture of all the pixels that comprise the Lez topographic basin. The annual relative amplitude is about 40 % between the driest days (Hu2 = 50) and the wettest days (Hu2 = 90).

Model calibration
The model was calibrated from twenty-one available flood events. Model inputs are either radar rainfall values corrected by the MFB coefficient or rainfall measured by recording rain gauges interpolated by the Thiessen polygon method. The model was applied using 75 m cell size and 1 h time step. Parameters remain constant for all the events, except for S. The Ds parameter has been first calculated: as a mathematical property of the model, it can be considered as the coefficient of an exponential recession curve, which could be successfully fitted for the most part of the events. The mean value was found to be 0.28 d −1 , with a 0.12 d −1 standard deviation. Because V 0 and K 0 were found to be dependant, the K 0 parameter was first empirically set to 0.3. Then both V 0 and S parameters were calibrated for each event, by Nash criterion maximisation: The calibration was driven using w = 0 in first approximation, and considering only discharge values greater than 50 m 3 s −1 in order to focus the calibration over the flood peaks, where the influence of w is negligible. The V 0 values were found to be near the mean value V 0 = 1.3 m s −1 , (with a 0.2 m s −1 standard deviation), and this value was finally kept for all the events. This value is fairly in agreement with the high flows velocity in the main stream of the catchment, and gives a posteriori some justification for K 0 .
Then both w and S parameters were calibrated for each event, using the Nash criterion over the whole event and constant values for Ds, V 0 and K 0 : the mean value was found to be w = 101 mm, with a 31 mm standard deviation. Finally, S was calibrated again using mean values of Ds, V 0 and w over the whole event by Nash criterion maximisation on discharges.
Radar rainfall data inputs greatly improve the simulations of the September 2003 and September 2005 floods (Table 4) when the Prades recording rain gauge, located in the topographic basin, did not work. They also improve the simulation in October 2001 and October 2008. In addition, radar rainfall can lead to a significant change in the estimation of S, like in September 2000 because the Prades recording rain gauge did not represent the total storm depth fallen in the basin. This was not the case for the other events whose radar rainfall intensity values were of poor quality. We thus kept simulations obtained with radar rainfalls data for these five events at the beginning of autumn and the simulations obtained with rain gauges rainfalls for the other events. When both radar treatments were available (September 2003 andSeptember 2005), Hydram was selected for homogeneity sake. Table 5 sum up the performances of the selected simulation with the Nash criterion but also in terms of errors in peak discharge, with the E RPD criterion (Eq. 13), simulated time of peak with the Etp criterion (Eq. 14) and bias on discharge volume (Eq. 15): Generally speaking, for eighteen floods, the model simulates well (with a Nash ranging between 0.66 and 0.94) with a set of parameters that remains constant for all the events (Fig. 5), except the initial condition S.
On average, the model overestimates the peak discharge by about 15 % and underestimates the volume by about 12 %. The error on time to peak is generally about 1 h. Moreover, the model initial condition, S, shows high values for the September events after a dry period, and low values for other events occurring after the first autumn rainstorms. As expected, the model initial condition depends on the season and on previous rainfall.
The optimal initial condition of the model was then correlated with three predictors of the initial wetness state of the basin: Hu2 index from the Meteo-France SIM model, base flow, and piezometric levels. The values of these predictors were taken at 06:00 UTC at the beginning of the simulations.
The base flow at the outlet of the basin was also compared with the initial condition of the model. With a determination coefficient R 2 = 0.66, the correlation was significant for a threshold α lower than 0.05 (Fig. 7).
Finally, the initial condition of the model was compared with the aquifer water level measured at an hourly time step. Correlations were significant for a threshold α lower than 0.05 for nine of the 12 piezometers with a determination coefficient ranging between 0.50 and 0.81. The best correlations are obtained with the three piezometers (number 5, 9 and 12) located near the Matelles fault which would represent a main drain of the karstic aquifer (Karam, 1989) and the piezometers located near the Lez spring (number 1, 3, 6, 11) are greatly influenced by its overflow. The linear regression established between S and the piezometric level at "Claret" (piezometer 9) represents a mean behaviour with a R 2 = 0.71 (Fig. 8). Note that only 12 events are available to establish this linear regression.
The robustness of the model was controlled, by performing calibration and validation over 3 split-samples for each presented predictor. For Hu2 index and base flow predictors, each calibration split-sample was constituted from 14 events (events 1 to 14 for the 1st calibration split-sample, 1 to 7 and 15 to 21 for the 2nd calibration split-sample, 8 to 21 for the 3rd calibration split-sample) and each corresponding validation split-sample was made of the 7 events not included in the calibration split-sample. For the piezometric level at "Claret", given that the events 1 to 8 and 10 are not available, each split-sample was constituted from 8 events (events 9 to 17 (except event 10) for the 1st calibration split-sample, 9 to 13 and 18 to 21 for the 2nd calibration split-sample, 14 to 21 for the 3rd calibration split-sample) and each corresponding validation split-sample was made of the 4 events not included in the calibration split-sample (Table 6). The method consisted first in testing the robustness of the regression S-Hu2 derived from each calibration split-sample; and then in comparing the Nash criteria on both corresponding calibration and validation split-samples: it should be noted that the Nash criterion have been processed for each event using the S value estimated by the regression S-Hu2 derived from the calibration sample. This allows making significant the comparison between calibration and validation split-samples. more recent the events are the higher the determination coefficients. When the event 1 to 14 are used, determination coefficients are low for both Hu2 index and base flow (R 2 = 0.41 for Hu2 index and R 2 = 0.54 for base flow) whereas they are higher using more recent events. This can be explained by the improvement of rainfall spatialisation: between 1994 and 2002 (event 1 to 12), only the Prades hourly rain gauge is used to estimate the rainfall on the catchment. Since 2002, the use of the Saint-Martin hourly rain gauge has allowed a better estimation of the rainfall on the catchement and, consequently, of the initial condition S. For all the predictors, median Nash are ranging between 0.69 and 0.83 for calibration and between 0.49 and 0.85 for validation. These results give an acceptable confidence in the robustness of either the predictors regressions or the rainfall-runoff model. Syst. Sci., 12, 1119-1133, 2012 www.nat-hazards-earth-syst-sci.net/12/1119/2012/  6. Linear regression between model initial condition S and Hu2 index of SIM model. The dashed lines represent the 95 % confidence limits.

Model error and uncertainties
In this part, the objective is to use the model in an operative way, i.e. using calibrated parameters and S-predictors relationships. The Nash values associated to these simulations give a more realistic estimation of the expected model error, as it could be for future events. Thus, additional runs were performed using w = 101 mm, Ds = 0. (ii) S = a.Q b +b, (iii) S = a.P z + b for the 12 events which have available data for the three predictors. Corresponding Nash values have been reported in Fig. 9 for each event, the values being firstly classified in descendant order. Optimal Nash values have also been reported, using the same representation.
It can be seen that naturally, optimal Nash values are higher than those derived by using S-Hu2, S-Qb or S-Pz relationships. For the latter, the mean Nash values are respectively 0.75, 0.76, 0.72, and does not appear to be significantly different.
Piezometric levels at "Les Matelles" or "Bois Saint-Mathieu" appear to be the most efficient predictors to initialize the rainfall-runoff model but the number of events of these two regressions (14 events available for "Les Matelles" and 12 for "Bois Saint-Mathieu") is less than the 21 events available for the base flow or the Hu2 index regressions. The small number of events enabled significant correlations to be established with a signification threshold smaller than 0.05 but not a significant hierarchy between the different Table 6. Results of the split sample tests performed with the three predictors (Hu2 index, base flow Q b and Claret piezometric level). "a", "b" and "R 2 " are respectively the slope, the intercept and the determination coefficient of the regression established between the initial condition S and the given predictor using the events of the calibration sample.  predictors. Indeed, variance analysis showed that the correlation coefficients of the regressions R established with piezometry, Hu2 and base flow (Q b ) are not significantly different. The assumption "R(piezometry) = R(Hu2) = R(Q b )" is rejected in all cases at a 0.05 and 0.10 signification threshold.
Model errors result, in a large part, from the uncertainties associated with the regressions between S and external predictors, i.e. errors concerning the estimation of either predictors or of the optimal initial condition given by the model. Indeed: -The Hu2 index, as an output of the ISBA surface model, depends on meteorological or geographical forcings (assumptions concerning soil structure and texture or vegetation cover), which are subject to uncertainties; -During low flow periods, base flow is very sensitive to rating curve shifts which can occur after major floods.
Base flow may thus not be known precisely. Moreover, the base flow at Lavalette can be artificially modified because of the exploitation of the ground water that supports the Lez river. This modification is more important in summer when the base flow is weak; -The piezometric level is probably the local predictor known with the greatest accuracy. Its ability to account for the filling of the whole aquifer may be jeopardized by both the heterogeneity of the karstic medium and the measurement of local features. Nevertheless, the good correlations for most of the piezometers indicate good representativity of the water content for the whole aquifer; -Estimation of the initial condition S is also subject to different biases. Firstly, uncertainties concerning rainfall estimation for the basin have direct repercussions on estimation of the initial condition. Next, uncertainties concerning discharge benchmarking (notably those that are greater than 300 m 3 s −1 in the case of the Lez) directly affect estimation of the initial condition. Finally, model concepts and calibration of the parameters can also influence estimation of the initial condition.
Thus, there are many sources of uncertainties that must be carefully taken into account when interpreting results. However, correlations are mainly modified by random errors, and not by systematic errors. The latter are "filtered" during the regression calibration between S and the different predictors. To improve the relations, the main uncertainties we have to correct are thus random uncertainties: rainfall estimation as model inputs, possible rating curve shifts in low water, or meteorological forcing of surface models. Syst. Sci., 12, 1119-1133, 2012 www.nat-hazards-earth-syst-sci.net/12/1119/2012/

Contribution of the groundwater to the flood
As mentioned above, groundwater is assumed to contribute to runoff volume in December. Furthermore, hourly piezometric data enabled dynamics of the karst aquifer to be compared with the modelled store content stoc(t), which controls runoff. Both variables showed the same time-variation shape, but piezometric level always reacted some hours later (from 3 to 20 h) both in September and December (Fig. 10). This delay suggests that the groundwater contribution to surface runoff does not affect the first peak of the flood, but could have much more impact on following peaks in the case of a multi-storm event. This behaviour has already been observed in the region in the Nîmes karstic basin (Maréchal et al., 2008) and the Coulazou watershed (Bailly et al., 2008a).

Conclusions
The distributed event-based parsimonious rainfall-runoff model used in this study proved to satisfactorily reproduce the Lez karstic basin flash floods. It is based on a limited amount of data, mainly storm rainfalls, and contains a reduced number of parameters, which were constant for 21 events that occurred in the period from 1994 to 2008, except for the initial condition S. The distributed model makes it possible to take advantage of the radar rainfall data without increasing the number of parameters or the complexity of the model. These data were shown to significantly improve most of the simulation of floods that occur at the beginning of autumn, and to compensate for the low density of rain gauges in the basin, as well as the observation gaps. At the end of autumn, the weak vertical extension of the cloud and the low altitude of the vertical heterogeneity of the reflectivity field probably limit the efficiency of radar measurements.
The initial condition of the model is fairly connected to three predictors of the basin wetness state: the base flow, the Hu2 index from the Meteo-France SIM model and piezometric levels. Piezometric levels supply the best correlation (R 2 = 0.81) with the initial condition of the model, but the limited number of events available did not allow a significant hierarchy of the performances of these predictors to be established.
The small number of parameters required for this model, along with its ability to be initialized using readily available data means a wide range of operational hydrological applications can be envisaged such as long-term assessment or real-time forecasting of flash floods in karstic basins. The study now needs to be extended to other karstic catchments to confirm these conclusions.