Effect of the spatial distribution of physical aquifer properties on modelled water table depth and stream discharge in a headwater catchment

Water table depth and its dynamics on hillslopes are often poorly predicted despite they control both water transit time within the catchment and solute fluxes at the catchment outlet. This paper analyses how relaxing the assumption of lateral homogeneity of physical properties can improve simulations of water table depth and dynamics. Four different spatial models relating hydraulic conductivity to topography have been tested: a simple linear relationship, a linear relationship with two different topographic indexes, two Ks domains with a transitional area. The Hill-Vi model has been modified to test these hypotheses. The studied catchment (Kervidy-Naizin, Western France) is underlain by schist crystalline bedrock. A shallow and perennial groundwater highly reactive to rainfall events mainly develops in the weathered saprolite layer. The results indicate that (1) discharge and the water table in the riparian zone are similarly predicted by the four models, (2) distinguishing twoKsdomains constitutes the best model and slightly improves prediction of the water table upslope, and (3) including spatial variations in the other parameters such as porosity or rate of hydraulic conductivity decrease with depth does not improve the results. These results underline the necessity of better investigations of upslope areas in hillslope hydrology. Correspondence to: C. Gascuel-Odoux (chantal.gascuel@rennes.inra.fr)


Introduction
In catchments underlain by crystalline bedrock with a deep, thick, weathered aquifer, runoff during storm events as well as during baseflow periods is often controlled by shallow groundwater.Water table depth and its dynamics control both water transit time within the catchment and solute fluxes at the catchment outlet.Water transit time depends on the slope position, from a few days in the riparian zone to a few years in uplands.Transit time in the uplands can be so long because of a much thicker unsaturated zone and also much lower hydraulic gradients on plateaus compared to midslope regions and riparian zones (Freer et al., 1997;Molenat and Gascuel-Odoux, 2002).While mean transit time can be estimated by water dating techniques or lumped models using input and output signals, modelling the internal flow within the catchment is necessary because of considerable uncertainty with these methods.Coupling water and solute transport requires a good knowledge of flow velocity and geochemical processes, which may vary laterally within the shallow groundwater.In such cases, solute fluxes such as nitrate are also controlled by hydraulic gradients within the groundwater that determine the contribution of different spatial domains and their connectivity (Ocampo et al., 2006a(Ocampo et al., , 2006b)).As an example, in upslope areas shallow groundwater can store nitrate, and in downslope positions can act as a nitrate buffer, with low N concentrations due to denitrification (Molenat et al., 2002;Steinheimer et al., 1998).Therefore, reasonable predictions of the water table depth and its dynamics in space and time are urgently needed for shallow groundwater catchments.
C. Gascuel-Odoux et al.: Effect of lateral variation of physical porperties on water table depth Experimental studies have shown that water table depth is more or less correlated with topography, in some cases strongly represented by topographic indexes (Moore and Thompson, 1996;Thompson and Moore, 1996) and in other cases less predicted by topographic indexes (Myrabo, 1997;Seibert et al., 1997).Some studies underline that this topographic dependence is not true for all storm events (Jordan et al., 1997) or for the whole hillslope (Molenat et al., 2005;Martin et al., 2006).The relationship is particularly weak in upslope areas where the water table depth is independent of topography (Molenat et al., 2005).Few studies have investigated the dynamics of water table depths in upslope areas because of the difficulties and costs associated with measuring relative deep water tables in areas not important for groundwater resources.The clear lack of correlation between topography and water table depth in some cases indicates that soil and aquifer properties cannot be assumed to be homogeneous for the entire hillslope and that the water table dynamics may also depend on variations of the soil and sub-soil physical properties.
Modelling studies have investigated the effects of vertical variations in soil properties on water table dynamics and discharge.Different conceptualisations were proposed for modelling vertical variation of soil properties.Different functions of hydraulic conductivity that decrease with depth have been tested by Ambroise et al. (1996).
More recently, assumptions about the drainable porosity (the pore space between field capacity and saturation) decreasing with depth and preferential flow have also been proposed and tested (Weiler andMcDonnell, 2004, 2006).The occurrence of shallow groundwater flow in a catchment is governed by a decrease of soil hydraulic conductivity and drainable porosity with depth.Under these conditions, percolation is stopped or limited and a saturated zone develops that generates lateral groundwater flow in response to a downslope hydraulic gradient.In these models and most other hillslope models, soil and aquifer properties were assumed to be laterally homogeneous, whereas many field observations show lateral variations in soil depth, porosity and conductivity.These variations are local (Tromp van Meerveld et al., 2007) or linked to topography with respect to soil (Curmi et al., 1998) or weathered materials (Dewandel et al., 2003(Dewandel et al., , 2006)).Few modelling studies have investigated the effect of these lateral variations.Saulnier et al. (1997) showed that variation in soil depth with topography does not really affect discharge.In contrast, local variation of soil depth can have a major effect on flow connectivity and thus discharge, as demonstrated by field observations and modelling approaches (Tromp van Meerveld and Weiler, 2008).Lateral variation in transmissivity was modelled using a Topmodel approach and calibrated with field observations (Lamb et al., 1997;Seibert et al., 1997).They produced slightly better agreement between simulated and observed water table depths.These studies focused on the bottom domain.
In saturated conditions, different pedogenic processes (e.g., lixiviation and redox) can induce topographic variation in the physical properties as observed in these studies and others (Curmi et al., 1998).Other kinds of processes have not been investigated.Particularly, detailed measurements of hydrodynamic properties as well as comprehensive and functional models of the hydrodynamic structure of weathered layers in crystalline rocks have just been started (Chilton and Foster, 1995;Taylor and Howard, 2000;Marechal et al., 2004;Dewandel et al., 2006).They show that hydrodynamic properties are not homogeneous in space but vary with the topography due to weathering processes.
There is still much debate about which modelling approach is appropriate to predict water table dynamics simultaneously with runoff in catchments dominated by saturated subsurface flow.One of the main questions is how far we can simplify the geometry of the system, the soil and sub-soil, as well as the hydraulics of groundwater systems, while still achieving an acceptable simulation of water table variations as well as transit time and solute fluxes.Fully distributed field studies and data are lacking for conceptualising these spatial variations and for testing them in subsurface flow models.Observations of water table dynamics should not be limited to the riparian zones or toeslopes as has primarily been the case so far.In upslope areas, we still lack a satisfactory understanding of the relationships between spatial and temporal variations in water table geometry and the geometric properties of catchments (topography, spatial distribution of soil depth and physical properties).
The aim of this work is to analyse how relaxing the assumption of homogeneity can improve simulations of water table dynamics, particularly for upslope position, by testing the effect of spatial dependence of hydraulic conductivity on topography.The study site is a catchment underlain by schist crystalline bedrock, where the shallow and perennial groundwater develops mainly in the weathered saprolite layer.On this site, different modelling approaches based on the homogeneity hypothesis were unable to efficiently simulate the water table depth over the whole hillslope (Molenat et al., 2005), thus it was considered necessary to search for an improved prediction of water table depth by considering different spatial models of physical properties in agreement with the observed structure of soil and bedrock properties.
The mean annual precipitation over the last 30 years is 909 mm, whereas the mean annual evapotranspiration and runoff from 1994 to 2004 (period of availability of local observations) are 709 and 400 mm, respectively.Three main formations are identified from the soil surface to the bedrock composed of Brioverian schist: the soil, the unconsolidated weathered bedrock, and the fissured and fractured weathered bedrock.The soils are silty loams, depths ranging from 0.5 to 1.5 m (Curmi et al., 1998).The soil system comprises an upland well-drained domain, with an average saturated hydraulic conductivity of 10 −5 m/s, and a poorly drained bottomland with an average hydraulic conductivity of 10 −6 m/s.The thickness of the unconsolidated material varies greatly in space from a few metres to 30 m.The boundary between weathered and fissured bedrock is poorly known and difficult to investigate.A model of the extension of the weathering processes was proposed for this region (Wyns et al., 1999;Dewandel et al., 2006).In this model, the unconsolidated weathering areas are more extended upslope than downslope.The hydraulic conductivity of the unconsolidated material is similar to that of the soil (Molenat and Gascuel-Odoux, 2002).A shallow, permanent and unconfined aquifer develops over the whole catchment area in the soil of bottomlands along the stream channel and in the unconsolidated weathered bedrock in the entire catchment.
The catchment is equipped with stream gauge stations at the outlet, a meteorological station and transects of wells that intercept the permanent and shallow groundwater that develops in the soil and unconsolidated weathered bedrock.Stream discharge was recorded every 1 to 6 min at the gauging station.The meteorological station recorded hourly rainfall, air temperature, soil temperature at 50 cm depth and other variables (wind speed, global radiation and relative humidity) required to calculate potential evapotranspiration using the Penman method.Wells dispatched on three transects and ranged from 1.5 to 20 m in depth are described in detail in Molenat et al. (2005Molenat et al. ( , 2008)).Water table depths were monitored with bubble sensors or shaft encoders with an integral data logger.Errors in water-table measurements with sensors ranged from 1 to 5 mm.The present work focuses on one transect named G consisting of six wells, from PG1 in the bottom to PG6 at the top (Fig. 1).This transect has been chosen because it has already been used to simulate water table depth without spatial structure of physical properties in a previous modelling study comparing water table simulations with three models -Topmodel, a kinematic model and a diffusive wave model -showed the difficulty in correctly predicting water table depth in the upslope area with all three models.Despite a slight improvement with the diffusive model, differences between observations and simulations remained large (Molenat et al., 2005).This transect presents a linear slope, without any slope effect of landscape elements such as hedge row or ditches which could interact with topography.The deeper wells (PG4, PG5 and PG6) remain located in unconsolidated weathered bedrock.They are not in the fissured layers: this has been clearly identified when drilling the wells.Slug tests have been performed on PG 2, 4 and 5 (Molenat, 1999).The hydraulic conductivity of the weathered layer was found to be variable over one order of magnitude from 4.10 −5 to 2.10 −6 m/s.Estimated saturated hydraulic conductivity at these three locations does not exhibit a spatial distribution but any rigorous analysis of the spatial distribution of Ks would require estimations at much more locations in space.
The observed stream discharge followed a classical pattern for humid and temperate climates underlain by shallow impervious bedrock (Fig. 2): the discharge is at its maximum in winter (December-January) and decreases until early autumn (September-October) when it dries out in most years.The water table depths follow different temporal patterns (Molenat et al., 2008) (Fig. 2).During a first period extending from late autumn to spring, the water table in the riparian zone ranges from the soil surface to around 0.5 m below the soil surface and water table rise occurs rapidly with each rainfall event.Positive values are recorded only in riparian zone, mainly at PG1, PG2 and PG3 locations, for short times.The positive values did not exceed 20 cm and correspond to flooding of the riparian zone.As a consequence of the flooding, the hydraulic heads in groundwater were larger than the soil elevation leading to Fig. 2. Daily precipitation, discharge and water table depth in the five wells for the water year 1999-2000; the first day is 1 September 1999 (Adapted from Molenat et al., 2005).
record positive values of water table depths.During winter periods, without any rainfall, the water table remains in the upper 50 cm and no recession can be observed.In the hillslope shoulder and plateau area, the groundwater table varies widely.The water table varies from 1 m to 3 m below the soil surface in the shoulder (PG4 and PG5) and from 4 to 6 m in the plateau area (PG6).Water table response to rainfall is much slower and smoother compared to the riparian zone.The water table falls as soon as rainfall ceases.During a second period which corresponds roughly to the summer period (July to September), the water table is far below the surface along the entire hillslope.Depletion in the riparian zone drives the water table down to 1 m below the soil, with very weak temporal dynamics in contrast to the winter period.In the shoulder and plateau area, recession of the water table is rapid and the final water table is very deep.An abrupt shift is observed from the first to the second period.Consequently, the hydraulic gradient can vary by a ratio of 1:5 during the year in upslope areas, whereas it is much lower in the midslope area and the riparian zone, close to the topographic slope.Thus, the hydraulic heads in the upslope areas were not proportional to the topographic slopes and no clear relationship could be observed.
The variations of the water table depth along the wet season are studied by analysing correlations between daily discharge and water table depth at the six wells (Table 1) and temporal distributions of the water depth (Fig. 3a), calculated on one wet season (from October 1999 to March 2000).The wells can be grouped in two sets: (1) PG2 and PG3, highly correlated (r 2 =0.97), water table depth always located above 1 meter, i.e. in soil layer, explaining our focus on only one of the two wells (PG3); PG1 can be assimilated to this group while the range of variations of water table depth is similar, but it has been moved aside because the water table depth shows little response by variations due to rainfall events, leading to a same and moderate correlation between discharge and water table compared to the other wells; (2) PG4, PG5, highly correlated (r 2 =0.95), water table depth ranging from the soil surface to 5 m deep, and PG6 more moderately correlated (r 2 =0.89 and 0.91 with PG4 and PG5, respectively) but the most correlated to the stream discharge (r 2 =0.80), water table depth ranging from 1 m below the soil surface to more than 7 m deep, explaining our choice of two wells (PG5 and PG6) in this second group.The water table in these six wells presents a high range of variations.Variations in PG6 directly control the seasonal variations of the stream discharge (Fig. 2).
The relation between the water table rise and the cumulative rainfall event amount, calculated for every rainfall event of five wet seasons (from January to March from 1998 to 2002) is linear (Fig. 3b).The water table recharge is quick and high.The slope of the linear relation allows us to deduce the drainable porosity.It is 6.4, 5.8 and 3.1% for PG4, PG5 and PG6, respectively, i.e. slightly decreasing from PG6 to PG4.

Model
We used the physically-based hillslope model Hill-vi (Weiler and McDonnell, 2004;Weiler and McDonnell, 2006;Tromp-van Meerveld and Weiler, 2008;Anderson et al., 2009a,b) for runoff generation by simulating water fluxes in the saturated and unsaturated zone of the hillslope.The model is based on the concept that two compartments define the saturated and unsaturated zone for each hillslope grid cell, based on topography and soil depth.
The unsaturated zone is defined by the depth from the soil surface to the water table and its time-variable water content.The saturated zone is defined by the depth of the water table above an impermeable interface and total porosity n.Lateral flow in the saturated zone is calculated using the Dupuit-Forchheimer assumption.Routing is based on the grid cell-by-grid cell approach (Wigmosta and Lettenmaier, 1999).Hydraulic conductivity is defined by an exponential function to reproduce the changes of hydraulic conductivity with depth due to soil development.Since the active hydrologic zone at the studied hillslope is relatively deep (up to 10 m), we included a term for constant hydraulic conductivity at depth in our model.The transmissivity T is then given by where Ks is saturated hydraulic conductivity, K o is the saturated hydraulic conductivity at the soil surface, m is the rate of hydraulic conductivity decrease with depth and z is the depth into the soil profile (positive downward).
We include a constant hydraulic conductivity at deeper depths K c .In addition, the model includes a depth function for drainable porosity, taking into consideration that drainable porosity declines with soil depth (Weiler and McDonnell, 2004).Therefore, depth variations of K o and drainable porosity are both described by decreasing exponential functions.Each function comprises a parameter defining the change with depth as indicated in Table 2.The variation in depth is much slower for drainable porosity than for Ks.Actual evaporation from the unsaturated zone is calculated based on the relative water content in the unsaturated zone and potential evaporation.Drainage from the unsaturated zone to the saturated zone is controlled by a power law relationship between relative saturation in the unsaturated zone and saturated hydraulic conductivity at water table depth (Weiler and McDonnell, 2004).Further details on Hill-vi can be found in the work of Weiler andMcDonnell (2004, 2006) and Tromp van Meerveld and Weiler (2008).In total, 7 parameters were used to simulate the dynamics of the saturated and unsaturated zones.A factor of variation α(x) has been introduced into the Hill-vi model to simulate spatial variation in porosity, hydraulic conductivity or the m parameter.This factor is a multiplicative factor which can be added to these two variables: where x is the distance to the stream and Var o can be the hydraulic conductivity, the drainable porosity or the parameter that determines the exponential decline in hydraulic conductivity with depth (m).
We considered topographic variation in this factor α to be a first step in implementing heterogeneity.Topographic dependence can be due to saturated conditions in the toeslope.It can be positive, as previously described, with hydraulic conductivity or porosity decreasing towards the toeslope, such as in gleyic soils, or negative, instead increasing towards toeslope, such as in peat soils.This topographic dependence can also be due to spatial variation in weathering processes.
Four spatial models (i.e., four ways the factor α could vary with topography) have been tested: 1. the linear model, varying linearly with distance to the stream; the simplest model of the four; 2. the threshold model, which includes two areas with constant alpha joined by a short transitional domain, describing spatial variation linked to the plateau or bottom domain; 3. the mono index model, with variation according to a topographic index, computed with a unidirectional scheme, i.e., where the flow coming from one cell goes to a single neighbouring cell, mostly downslope; we assume that spatial variation of the physical properties follows water drainage flow pathways; 4. the multi index model, varying with a topographic index, computed similar to the mono index model but with a multi-directional schema, such that the flow coming from one cell goes to different neighbouring cells according to the slope gradient, assuming similar spatial variation due to water drainage.
For the two last models, α was calculated on the 3-D Digital Elevation Model (DEM).It was been normalised to compare to the four models.

Numeric experiments
One transect, with six wells (transect G) and hourly data from the year 1999-2000 (Fig. 2), was selected for this study because previously used to compare different models of the water table assuming homogeneity of soil and bedrock properties for the entire hillslope (Molenat et al., 2005).The values of the parameters were chosen within a relatively small range according to the observed values (Table 2; Molenat et al., 2005).The transect is represented as a 10 m wide and 480 m long area with an homogeneous depth of 10 m.This is an approximate depth, deeper than the water table depth not to influence it, and averaging few observations.The period is 3650 h long (from 12:00 p.m. on 31 October 1999 to 12:00 p.m. on 31 March 2000), is representative of an average wet season and covers seasonal variations, between fall and spring, as well as the quick responses of the water table to rainfall events in winter.The time step of simulation was one hour.Thus, the transect and the period represent a well-investigated starting point to test the improvement of predictions of discharge and water table dynamics when introducing spatial variation in soil and bedrock properties.
A 2 m grid resolution was used for the Hill-vi modelling.We chose mean winter conditions (a daily rainfall of 3 mm per day) and applied them as constant rainfall over 2000 h to the model to develop the initial conditions for each model run.For the four spatial distributions of the physical properties, the factor of variation α was derived from a detailed topographic survey on a 2 m grid size.For the threshold model, the threshold was set to the slope break between the linear slope and the plateau (Fig. 4).For the two topographic index models, the DEM (with a 10 m grid size, plus manual measurements by digital Laser total station) was used with the general slope from the grid cell to the stream, according to the flow pathway, in place of a local slope, as proposed by Gascuel-Odoux et al. (1998) and Merot et al. (2006).This modification takes into account drainage of the shallow groundwater according to the slope position.The normalised values of the factor of variation are different, particularly in the upslope domain (Fig. 4).
The specific stream discharge measured at the outlet was considered as the specific discharge of the study hillslope, assuming that hillslope contribution was uniform along the stream network length.The assumption is commonly used in the literature since measuring hillslope runoff is very difficult or impossible.
The effect of the spatial distribution of physical properties was analysed with two methods: -Firstly, a Monte Carlo procedure (1000 parameter sets) was used to get the best set of the 7 parameters for each spatial model considering discharge on one hand and the water table on the other hand.We analysed the improvement of predictions when considering spatial variation of the physical properties based on the best set of parameters.A uniform probability distribution was assumed between the lower and upper limits of the parameters (Table 2).As a first step, only spatial variation in K o was studied, with its variation ranging over one order of magnitude and with K o increasing from downslope to upslope, according to the weathering processes in the upslope area observed by Wynns et al. (1999) and Dewandel et al. (2006).In addition, the combined lateral variation of K o with both m or porosity was investigated.No spatial variations in Kc was studied because it was assumed that its spatial variation is randomly distributed, according to the location of the fissures, and not topo dependent as for K o , m or porosity.
-Secondly, a sensitivity analysis was performed to explore the effects of the direction and magnitude of spatial variation in the physical properties.A fixed set of parameters was chosen which corresponds to the best predictions related to the water table depth considering no spatial variation (Table 2).Different ranges of variation in K o were tested, with a magnitude from 1 to 4, then the location of the threshold in the threshold model was moved up and down 10-40 m.
The performance of each model run was evaluated with the following objective functions: -For discharge, the Nash-Sutcliffe efficiency (Nash and Sutcliffe, 1970) and the efficiency of the logarithmic discharge values were calculated by: where σ g err and σ g obs are the variance of the simulation errors and observations, respectively, and λ 2 err and λ 2 obs are the variance of errors calculated from logarithmic discharge and the variance of logarithmic observed discharge, respectively.
-For the dynamics of the water table, the average difference in the average level (D) and the range of variations (R) between the observed and simulated water table variation for three wells (PG2, PG5 and PG6) were calculated with: where WTD is the water table depth, and sim and obs correspond to the simulated and the observed values.D and R were computed for all non-missing values during the study period.For R, logarithms are calculated to get two objective functions (R, D) with a target value equal to (0, 0).These two criteria are simple and keep a physical significance: they allow us to assess whether water table depth and its range of variation agree with the observed values, to evaluate the error in the studied period without focusing on detailed response dynamics.We have preferred to calculate D and R on the relative values rather than on the absolute values, despite a possible compensation of the errors.In our case, the compensation of errors cannot probably be involved because not involved in a previous work in which the errors were mainly due to a shift or a smoothing of the mean level of the water table depth (e.g., Molenat et al., 2005).Only behavioural simulations (i.e., simulations leading to a Nash-Sutcliffe efficiency greater than 0.4 and to a logarithmic discharge efficiency greater than 0.7) were retained for our analysis.

Hill-Vi model application
The application of the Hill-vi model to the data set without considering any spatial variation of the physical properties results in a best fit with a similar efficiency for discharge as in previous studies.Molenat et al., (2005) calculated efficiencies of 0.87 (discharge calibration) and 0.82 (groundwater calibration) and of 0.76 (discharge validation) and 0.68 (groundwater validation) (validation period: December 2000 to April 2001), with a diffusive model compared to Eff-Q=0.59 and Eff log-Q=0.77(discharge) and to Eff-Q=0.43 and Eff log-Q=0.73 with Hill-vi (groundwater) (Table 3).The Hill-vi model seems to produce a similar behaviour, and the effects of spatial variation of the physical properties can be tested.
The best fit of the model to water table depth results in a slight decrease of the discharge efficiency, but the objective functions for the water table depth are improved (Mean D=1.83, Mean R=−0.05, fit to discharge; Mean D=1.67, Mean R=−0.04 fit to groundwater) (Table 3).The range of variation of the predicted water table depth is higher in the upslope area (D-PG6=1.52,R-PG6=-0.18)compared to the simulation with the best fit for discharge (D-PG6=1.71,R-PG6=−0.22)(Table 3).However, the simulated seasonal variations at PG5 and PG6 are much smaller than the observed variations for both objective functions (   domain, but are still not satisfactory.The Hill-vi model, as well as previous modelling approaches, fails to predict the observed dynamics of the water table, particularly in the upslope area (Fig. 5).This result also justifies the present approach of testing the effects of lateral variation of physical properties.

Effect of lateral variation of saturated conductivity on discharge and water table depth
The two criteria R and D, related to the water table depth, have been computed for the different spatial models.In Fig. 6, the behavioural model simulations with Eff>0.4 and Eff log-Q>0.7 are plotted for the three wells and compared with the different spatial models.The optimum value of zero for the two objective functions R and D cannot be reached when no spatial variations are taken into account, while the two criteria are closer to this target when considering any of the 4 other spatial models.However, a similar good agreement for all three wells can never be achieved simultaneously.The results are rather similar for the two topographic index models and the linear model.Thus, these three models cannot induce any significant differences in estimates of water table dynamics.Water table depth is best estimated with the threshold model, essentially because of a better prediction of PG5 and PG6.The number of behavioural simulations is lower for this spatial model compared to the three other models.The threshold model is the most selective model for estimating water table depth with a small number of parameter sets.
When analysing the result for the best set of parameters (Table 3), similar conclusions can be drawn: the best simulations of water table depth are obtained with the threshold model.The depth of the water table is always estimated well for PG2 (the mean error is lower than 0.1 m for all models), but only estimated with an error lower than 1 m for PG5 and PG6 using the threshold model (D=0.96Table 3. Mean error (D) and log of the relative range of variation (R) of water table depth during the testing period for three wells -PG2 (bottom), PG5 (midslope) and PG6 (plateau) -and their average values (Mean D and Mean R), as well as discharge efficiency and Log discharge efficiency, with different spatial models of variation in physical properties.7 synthesises these results and shows that any model is able to simulate discharge relatively well with similar efficiency.The spatial models have no significant effect on predicting discharge.The effect of the spatial models on water table depth is small for the downslope area (PG2).Conversely, the spatial models have a larger effect for the upslope area, particularly for estimating mean water table depth, when fitting to discharge, and for estimating the range of variation of water table depth when fitting to the water table.However, if the criteria are averaged, the effects are much smaller (Table 3).

K and drainable
This analysis indicates that variation in K o , according to the threshold model, results in the best prediction of water table depth in the upslope and downslope areas.When the model was fitted only against discharge, as it is generally done when no data on water table depth are available and no model of spatial structure is known, the water table depth cannot be reasonably predicted in the upslope area for the observed hillslope.Conversely, knowledge of the spatial structure can improve the estimate of water table depth.
Table 4 shows the effect of the different spatial models on the parameter values for the model with the best fit.The variations among the different models cover all possible ranges when the simulations are fitted to discharge, except for K o and drainable porosity which are rather constant.The range of the variation is slightly wider when the simulations are fitted to water table depth.
The dynamics of the predicted water table depth for all wells are illustrated in Fig. 8.In the riparian zone, the predicted water table depth is smoother than the observations (PG1, 2 and 3) independent of the spatial model.Some of the observed positive water table depths also cannot be reproduced by the model, since Hill-vi in this version does not consider exfiltration and flooding.The dynamics of water table depth are better predicted in the midslope and the upslope area if a spatial model is introduced.In particular for PG6, though, the large dynamics of water table variation are not captured by any of the models.
C. Gascuel-Odoux et al.: Effect of lateral variation of physical porperties on water table depth 3 6 computed for the selected simulations (Efficiency > 0.4; Log Efficiency > 0.7), during the tested period, for three wells -PG2 (bottom), PG5 (midslope) and PG6 (plateau) -for different spatial models of variation in hydraulic conductivity.computed for the best simulation fitted to: a) discharge and b) water table depth, during the tested period, for three wells -PG2 (bottom), PG5 (midslope) and PG6 (plateau) -for different spatial models of variation in hydraulic conductivity.
Fitted on Water table Fig. 7. Mean error (D) and log of relative range of variations (R) of water table depth, computed for the best simulation fitted to: (a) discharge and (b) water table depth, during the tested period, for three wells -PG2 (bottom), PG5 (midslope) and PG6 (plateau)for different spatial models of variation in hydraulic conductivity.

Sensitivity analysis of the magnitude and direction of the spatial variations of the physical properties
Since all results for testing different spatial models with variation in saturated hydraulic conductivity were not satisfactory, a sensitivity analysis was performed to study the influence of possible additional effects.We tested the effects of the magnitude and direction of the variations in the spatial model, considering the different models as shown in Fig. 9 and using the simulation with no spatial variation as a reference.The value of saturated conductivity was fixed 4, 3 and 2 times higher going upslope, such that it increased from bottom to top, and 2 times higher going downslope, such that it increased from top to bottom, compared to a reference value.Fig. 9 shows that the magnitude and direction of the spatial model do not affect the prediction of discharge.For the threshold model, better predictions are obtained when K o increases from downslope to upslope than the reverse.For this model, larger magnitudes of spatial variation result in better prediction of water table dynamics in the upslope area (PG6).For the other three models, the improved prediction in one well is compensated for by worse predictions of the midslope area (PG5).
In addition, we also tested whether the location of the change between the two domains of the threshold model has any effect on performance, using the best-fit model for the  threshold model as a reference.Figure 10 shows that the best estimate of water table depth is obtained with a threshold located just before the break of the slope, 60 m downslope of the initial location.However, the performance increase is small and only visible for PG6.

Effect of lateral variation of K o , m and porosity on discharge and water table depth
The improvement of the water table depth estimation, considering spatial variation only of K o , is real, but the estimates are still far from the observations.Two hypotheses can be generated: (1) spatial variation of K o does not have a large effect because it concerns layers that are too superficial, compared to the depth of variation in groundwater in the plateau domain, due to a too small fitted value of m; (2) the reactivity is induced by both spatial variations of K o and drainable porosity.These two hypotheses have been tested C. Gascuel-Odoux et al.: Effect of lateral variation of physical porperties on water table depth 3 9 computed for a reference case considering a different range of spatial variation of saturated hydraulic conductivity from the top to the bottom and reverse, during the tested period, for three wells -PG2 (bottom), PG5 (midslope) and PG6 (plateau) -and different spatial models of variation in hydraulic conductivity.
-2 (Fig. 11, Tables 3 and 4), by decreasing m and drainable porosity from upslope to downslope, as previously for K o , to simulate deeper and more intense weathering processes upslope.
No improvement is observed.The discharge is always predicted as well as previously, when fitting to the discharge or to the water table.Water table depth is predicted better than without the spatial model, but not clearly better than when considering spatial variation of K o only.Therefore, the actual results due to only variation of K o can be considered to be the best solution in cases where no measurements are available.

Effect of spatial variability of physical properties on discharge and water table predictions
Previous studies of the Kervidy-Naizin catchment have shown the difficulties in predicting water table depth and its dynamics in the upslope area (Molenat et al., 2005).Soil physical properties have previously only been modelled as they change with depth but not laterally along the hillslope.This hypothesis of lateral homogeneity has been relaxed here in an attempt to better predict water table dynamics by including different spatial models describing variations of the physical properties.Intensive spatial studies concern the soil (the 1 m top layer) but do not account generally for the weathered layers (deeper that 1 m) due to the difficulties to measure hydraulic properties in such material.Since we knew that a spatial structure did exist but could not choose a-priori the right one, we chose to test different spatial structures considering the proposed four models.To establish a spatial structure, it would be necessary to get numerous slug tests, which is difficult from a practical point of view.All tested spatial models were related to topography to implement relatively simple spatial variations of the physical properties.The addition of such spatial models has improved predictions of water table dynamics in the studied hillslope.The threshold model, which defines two lateral domains of K o (one for the plateau and one for midslope and toeslope), has been shown to particularly improve prediction of water table dynamics.This spatial model agrees with Dewandel model (Dewandel et al., 2003(Dewandel et al., , 2006) ) proposing a structural model 4 1 spatial variation of Ks, m and drainable porosity and using two spatial models of saturated hydraulic conductivity, during the testing period, for three wells -PG2 (bottom), PG5 (midslope) and PG6 (plateau).Threshold variation K and porosity Fig. 11.Mean error (D) versus log of relative range of variations (R) of water table depth, computed for acceptable simulations (Efficiency>0.4;Log Efficiency>0.7),considering spatial variation of Ks, m and drainable porosity and using two spatial models of saturated hydraulic conductivity, during the testing period, for three wells -PG2 (bottom), PG5 (midslope) and PG6 (plateau).
for highly weathered hillslope.The weathering is more developed upslope and on the plateau, and therefore the weathered layer is thicker and its hydraulic conductivity is higher.Otherwise different spatial models had been previously compared, taking into account vertical variations of physical properties by including different vertical layers on this hillslope (Molenat and Gascuel-Odoux, 2002;Martin et al., 2006).These different hypotheses have been found to be not relevant for improving prediction of the water table.Thus, investigating and testing lateral variations of the physical properties appears to be a necessary step to better describe the spatial structure of this hillslope and to test the effect of different spatial domains of K o in relation to observed weathering processes along the hillslope.
However, the overall modelling results still remain unsatisfactory, and different reasons can be found.One explanation could be the effect of small-scale local variations of the physical properties.Local characteristics could explain the high reactivity of the water table following intensive rainfall events and also influence lateral transfer in the shallow groundwater.A second reason could be the too low range of spatial variation in K o .It would be interesting to enlarge the range of variation of the physical properties.However, since the combination of possibilities is huge, this would need to be done in a very structured way.The last potential reason is that the connectivity of the studied hillslope to the stream is not representative of the other hillslopes of the catchment and is therefore not related well to the discharge dynamics of the whole catchment.The studied hillslope presents a constant slope except a small plateau (after PG6), and a small convex domain (PG1 to PG3) which corresponds to a riparian zone (PG1 to PG3).It is included in a headwater catchment, with hillslopes presenting similar geometry and same land-use.Consequently it can reasonably be assumed that the hillslope contribution is uniform along the stream-length, including the riparian zone, and differences in both magnitude and timing are small.If not, we would have obtained bad simulations of discharge while achieving good simulations of water table depth when fitting the simulations to them.
Otherwise, internal processes, such as preferential flow or spatial variation in evapo-transpiration (ET) and wrong boundary conditions could also be evoked as reasons to explain such difficulties in modelling the water table depth.The cumulated values of rainfall and potential evapotranspiration (PET) over the study period are about 340 and 80 mm, respectively.Therefore the recharge and lateral flow processes are clearly dominant.If spatial variation of ET would have to be considered, it would increase from upslope (deeper water table) to downslope (saturated conditions), and be rather similar from PG1 to PG3 (saturated conditions), and from PG4 to PG6 (soil at C. Gascuel-Odoux et al.: Effect of lateral variation of physical porperties on water table depth field water capacity).But, firstly, the difficulties are to simulate the water table during the recharge periods, where ET is not an active process, as well during the recession periods, where ET is an active process.Secondly, the simulated values are rather correct in PG4, much higher in PG5 and PG6 while ET can be considered similar for these three sites.Therefore spatial variations of ET could exist but could not explain the difficulties in modelling water table depth along the hillslope.Concerning the preferential flow, this process is completely outside of our focus.However, if preferential flow occurred, it would be considered as similar for these three wells PG4, PG5 and PG6.At least, spatial variations of others parameters would have to be tested, particularly the boundary conditions.Aquifer boundary could not be similar to the catchment boundary defined from the topography.Furthermore, connection with a deep and regional aquifer we did not take into account in our modelling, could also affect the water table dynamics.The spatial variations of K c as well lateral preferential flow pathways in the weathered layers could significantly influence the spatial variations of the water table dynamics, as shown by Tromp-van Meerveld and Weiler (2008) in a heterogeneous soil/bedrock interface and by Anderson et al. (2009b) in the soil macropores, using similarly Hill-vi model.
Whatever the reasons and additional assumptions to be tested, we can obtain two major conclusions from this study: (1) There is a real lack of physical measurements of the weathered layers (e.g., depth, transmissivity and porosity), particularly on the plateau domain.These measurements are important because of their potential effect on the dynamics of the water table, as shown in this study.(2) Hill-vi is an interesting tool to investigate the effect of the spatial structure of the hillslope and its effects on hydrological processes.The Hill-vi model has been applied here in different conditions than previously (a 10-m deep transect, a high seasonal reactivity of water table depth, etc.).Finally, the prediction of discharge was acceptable, and the model has included relevant spatial structure for the hillslope, which have to be confirmed.

Interest for hillslope hydrology
As very often in hydrology, our experiments and their results are related to a specific place.However, the studied hillslope represents common characteristics of many deeply weathered watersheds, and more general assessments can also be drawn from this study.
As previously observed, the prediction of discharge was not greatly affected by spatial variations in physical properties, reversely to the prediction of the water table depth or other internal processes (see also Weiler and McDonnell, 2006).A complex spatial model to describe the physical properties is often not necessary for improving predictions of discharge, but it may be necessary for prediction of internal processes, water pathways and water transit time.This study showed that lateral variation in hydraulic conductivity seems to be as important as variation in soil depth for correctly representing the fluxes and dynamics in a hillslope.Therefore, in the studied environment, as well as in similar environments where the water table dynamics are different in the upslope and downslope areas, a simple, process-driven spatial model representing changes in saturated hydraulic conductivity could be a good solution to better predict discharge and water table dynamics when detailed data of water table depth are missing.However, it is necessary to infer a relevant spatial model and to parameterise it adequately.An inverse model approach (Carrera et al., 2005) cannot be the solution as long as we do not know the correct spatial model.The possible combinations of parameters are too numerous, and the unknown values of the parameters to describe the spatial model amplify the problem.In the past, inverse models have been applied to estimate groundwater flow by fitting them to many observed heads; however, spatial variability has been fixed to the different facies, and only their absolute parameterisation has been changed (Fienen et al., 2009).This spatial model can only be reliable if the correct spatial representation from soil or geological surveys is available.But, particular parameters like hydraulic conductivity are already highly variable within one soil or geological unit, so we can hardly use this information to represent absolute differences among the areas.Only a combination of water table measurements (or other spatial observations, such as soil moisture, if other processes are studied) and an appropriate spatial model, together with discharge observations, can provide the necessary information to parameterise a distributed model (Yeh et al., 2008).As can be seen for the studied hillslope, the predictions based on the different spatial models are rather different from each other and different from predictions without a spatial model.
In summary, we would like to stress the real need for more water table observations in the entire watershed.Generally, wells are drilled in the lower hillslope zone, while we would need deeper wells in the upslope area.Observations in the downslope area are often poor event responses as indicated in table 1, and therefore not sufficient to provide independent observations for benchmarking a model (Lyon et al., 2006).This study also shows that predictions of PG2 are close to the observations, because they are closely related to the event dynamics of discharge.Upslope information provides us additional information to calibrate a model while it is correlated to seasonal variations of discharge.This independent information is lacking in most experimental catchments (see Lyon et al., 2006 for an example).Until now, many hillslope studies have focused on an extent of 20-50 m.Extending the area of these studied slopes may help us to better characterise the behaviour of the entire hillslope.It is also important to ensure that our distributed Hydrol.Earth Syst.Sci., 14, 1179Sci., 14, -1194Sci., 14, , 2010 www.hydrol-earth-syst-sci.net/14/1179/2010/ models can predict the behaviour of the entire hillslope since the resulting transit times could be very different if the unsaturated zone in the upslope area is 1 m or 10 m deep.

Conclusions
Water table depth and its dynamics are often poorly predicted in upslope areas.We have analysed how relaxing the assumption of lateral homogeneity of physical properties can improve simulations of water table dynamics.Four different spatial models relating hydraulic conductivity to topography were tested for a well-studied catchment (Kervidy-Naizin, Western France).We used the Hill-vi model to represent the shallow and perennial groundwater that develops in the weathered saprolite layer.The results indicate that discharge and water table depth in the riparian zone are similarly well predicted by the four models, as well as with a model not considering the spatial variability.However, a spatial model including higher conductivity in the upslope area improves the prediction of the water table in this area.There could be more hypotheses tested with this approach, but we should constrain them with field observations.This study underlines the real need to better investigate the upslope areas in watersheds and the hydraulic properties of the weathered layers, particularly when questions of residence time and coupling water with solute transport are involved.Upslope information provides additional, independent information to calibrate a distributed model.This additional information has to be better analysed, particularly with a cross-analysis of discharge and water table dynamics in different hillslope positions.However, the purpose should always be to choose the simplest spatial structure of the hillslope, taking into account soil, weathered layers and bedrock.

Figure 3 .Figure 3 .Fig. 3 .
Figure 3. Variations of the water table depth: A) Temporal d depth in the six wells, calculated on five water years (1998 to water table depth versus rainfall amount, calculated for all th seasons (January to March, from 1998 to 2002 ) (N=88).

Figure 4 .Fig. 4 .
Figure 4. Normalized factor of spatial variation of the hydraulic conductivity from the b to the top domain, for different spatial models.

Fig. 5 .
Fig. 5. Temporal variation of watertable depth at the six wells in the water year 1999-2000, from 31 October to 31 March: observed values, simulated values fitted to discharge and fitted to water table depth considering no spatial variation of physical properties.

Fig. 6 .
Fig.6.Mean error (D) versus log of relative range of variations (R) of water table depth, computed for the selected simulations (Efficiency>0.4;Log Efficiency>0.7),during the tested period, for three wells -PG2 (bottom), PG5 (midslope) and PG6 (plateau)for different spatial models of variation in hydraulic conductivity.

Figure 8 .
Figure 8. Observed and modelled daily variation of water table depth related to the best fits regarding discharge, for six wells and different spatial models of variation in saturated hydraulic conductivity, from 31 October 2002 to 31 March 2003

Fig. 9 . 0 Figure 10 .Fig. 10 .
Fig.9.Mean error (D) and log of relative range of variations (R) of water table depth, computed for a reference case considering a different range of spatial variation of saturated hydraulic conductivity from the top to the bottom and reverse, during the tested period, for three wells -PG2 (bottom), PG5 (midslope) and PG6 (plateau) -and different spatial models of variation in hydraulic conductivity.

Table 1 .
Coefficient of correlation between discharge and water table depth in the six wells (PG1 to PG6) during the study period(October 1999 to March 2000).

Table 2 .
Values of the parameters in the Monte Carlo analysis and in the sensitivity analysis.

Table 4 .
The number of acceptable simulations of observations (Eff-Q>0.4and Efflog-Q>0.7) and the values of the parameters of the best-fitted simulation regarding discharge and water table depth, with different spatial models of variation in physical properties which all consider K 0 , m and drainable porosity increase from downslope to upslope.Gascuel-Odoux et al.: Effect of lateral variation of physical porperties on water table depth Figure 7. Mean error (D) and log of relative range of variations (R) of water table depth, Hydrol.Earth Syst.Sci., 14, 1179Sci., 14,  -1194Sci., 14,  , 2010www.hydrol-earth-syst-sci.net/14/1179/2010/ C.
table depth m