Spatial variability of the hydrodynamic properties of soil at hillslope scale in Cevennes area

Abstract:

Introduction

In the Cevennes mountainous area, flood processes are due to excess of soil water storage capacity depending on soil porosity, depth and water content at the beginning of the flood (Bouvier et al., 2007; Cosandey, 1994). For hydrological modeling, it is thus important to estimate the hydrodynamic properties of the whole vertical soil profile. In this study, the soils hydrodynamic properties were derived from inverse modeling of soil water content measurements. The efficiency of inverse modeling methods has indeed been proved to be satisfactory in the past years (Köhne et al., 2009; Vrugt et al., 2008; Wöhling and Vrugt, 2011). Assuming that soils have strong heterogeneity, water content was monitored along 4 hillslope transects (68 water content probes) on a small experimental catchment in order to grab the spatial variability of the hydrodynamic properties. The objectives of this study are i) to evaluate the robustness of the hydrodynamic parameters at the event scale depending on the different monitoring stations, ii) to study the spatial variability of the hydrodynamic parameters at the hillslope scale for a given transect, iii) to express this spatial variability as a statistical distribution, to be interpreted in terms of spatial organization.

1. Materials and methods

1.1 Study area and data collection

The Valescure area (44.1034°N, 3.8316°E) is a small headwater catchment of 3.83 km² located in the South of France at the southern boundary of the Cevennes mountain area. The altitude ranges from 244 m to 815 m ASL. The slopes are steep with an average value of 56%. The vadoze zone material is a Dystric Cambisol (2 mm, 25-45% coarse sand between 0.2 and 2 mm, 20-40% sand between 0.02 and 0.2 mm, 5-10% silt between 2 and 20 µm, and less than 2-3% clay < 2 µm.Four soil moisture monitoring transects were installed between 2009 and 2013, each transect dealing with 8 to 10 monitoring stations distributed every 10-20 m along a downstream-upstream axis. Transects vary in orientation, slopes, soil depth and land cover as shown in table 1. For each station, soil moisture was recorded from 2 ML2-X ThetaProbe probes (Delta-T Devices Ltd), one near the surface (15 to 20 cm depth depending on stations), one near the soil bottom position (25 to 40 cm depth depending on stations).

Rainfall was recorded in the center of the catchment with a tipping bucket rain gauge Precis-Mecanique 1000 cm².

Table 1: Main characteristics of the 4 hillslope transect

Transect	Orientation	Slope	Nb Stations	Monitored depths (min/max)	Land cover	Bedrock	Monitoring period
1	N.E	35%	10	20 / 40 cm	Chestnut-tree	Granite	2009 -2010
2	N.W	50%	8	15 / 40 cm	Chestnut/oak-tree	Granite	2010-2011
3	W	40%	8	18 / 35 cm	Chestnut/oak-tree	Granite	2011-2012
4	S.E	40%	8	15 / 35 cm	Chestnut/oak-tree	Orthogneiss	2012-2013

1.2 Soil water flow modelisation

The HYDRUS-1D (Šimůnek et al., 1998) model is used for soil water flow modelisation. It is based on the one dimensional Richards equation

[1]

where θ is the volumetric soil water content [L³.L^-3], t is time [T], z is the vertical coordinate (positive upward) [L], h represent the pressure head [L], K is the unsaturated hydraulic conductivity function [L.T^-1], and S is a sink term such as evapotranspiration [L³.L^-3.T^-1].

Soil hydraulic properties were parameterized by the Mualem–van Genuchten (MVG) (Mualem, 1976; van Genuchten, 1980):

if h≥hs if h<hs Â

[2]

[3]

with K_s  the saturated hydraulic conductivity [L/T], θ_r and θ_s [L³.L^-3] respectively the residual and saturated water content, α [L^-1] and n [-] both water retention curves shape parameters,  the effective saturation [-], h_s the air entry value and l [-] the porous connectivity factor which Mualem recommend to set as 0.5. In this study, we assumed, and.

1.3 Event scale data

An event-scale approach allowed considering both root water-uptake and evapotranspiration as negligible as regard to rainfall influence. For each transect, a set of rainfall (resumed in Table 2) events was selected according the following conditions:

- A minimum cumulated rain of 50 mm

- Rain intermittency during the event must not last more than 48 h

Table 2: Essential information about selected rainfall events for transects 1 to 4.

Â	Nb Episodes	Cumulated rain (min | max) [mm]	Rainfall duration (min | max) [d]	Max cumulated rain during 1h [mm]
Transect 1	13	50 | 928	2.5 | 21	42
Transect 2	10	70 | 928	2.6 | 20.6	46
Transect 3	5	54 | 225	2.6 | 14.2	19
Transect 4	10	50 | 315	3.4 | 27	33

1.4 Inverse modeling

Five parameters are thus necessary to define all the soil hydraulic properties: θ_r, the residual water content, θ_s, the saturated water content, K_s, the saturated hydraulic conductivity and α and n, two shape parameters.

Because of the interaction between model parameters, the inverse modeling methods are subject to non-uniqueness of solution, problem also known as equifinality (Beven, 1993; Minasny and Field, 2005). Therefore one possible solution is to reduce the number of parameters to calibrate as far as possible (Pollaco et al., 2008; Ritter et al., 2003).

Thus, it was decided to:

i) Consider the vadoze area as only one layer, involving the two probes at each station; this first layer ends 10 cm below the deeper probe. Underneath is considered a second layer corresponding to the weathered bedrock; this second layer is 150 cm deep.

ii) Preset as far as possible the parameters: by setting them to observed or measured values (Saturated water content θ_s for layer 2); because the model is considered little sensitive to the parameter (residual water content θ_r for layer 1 and layer 2); from a previous calibration of the model, retaining the mean calibrated value of the parameter as preset value for the final calibration (θ_s for layer 1, nfor layer 2),

iii) Analyze the dependency between the couples of parameters, and preset one of the two parameters in case of dependency. Both α and n parameters were found to be highly dependent at the event scale, thus αwas preset for both layers 1 and 2.

Top boundary condition is only controlled by precipitation (L.T^-1). Boundary condition at the bottom of layer 2 was considered as free drainage. Soil initial water content has been set as the observed value at the beginning of each event.

Finally, only 3 parameters are free to be calibrated: n and K_s for layer 1, K_s for layer 2.
For each event and each station, those parameters were obtained by minimizing the RMSE between both observed and simulated water content values at both probes positions of the considered station:

where N denotes the total number of time steps of the event, θ^surf the water content of the upper probe, θ^deepthe water content of the deeper probe. In order to avoid problem of parameter estimation due to a rough surface of the RMSE, a global algorithm was used, consisting in generating a lot of sets of parameters within min and max boundaries (Tab.2). For each event 10,000 parameters sets have been tested to find the best parameters. In order to enlarge the possible range of parameters, the 10 best RMSE have been retained in order to give an idea of the numerical variability of each parameter, around the optimal value.

Â

Table 2: Mualem - van Genuchten parameters min-max intervals used for generating the random parameters sets

		Layer 1					Layer 2
	ϴr [cm3/cm3]	ϴs [cm3/cm3]	α [mm]	n [mm-1]	Ks [mm/h]	ϴr [cm3/cm3]	ϴs [cm3/cm3]	α [mm]	n [mm-1]	Ks [mm/h]
min	0.005	0.36	0.012	1.1	200	0.05	0.25	0.008	3	1
max	0.005	0.52	0.012	5	2000	0.05	0.25	0.008	3	500

In bold the parameters to be calibrated. The saturated water content of layer 1 has been previously calibrated, and set to the mean value for each station

2. Results

2.1. Estimation of MvG Parameters

First, we analyzed the robustness of the estimation of the parameters for a given station. Robustness is expressed as a low variability of the parameters from one event to another. An example is given for the transect 1, highlighting the event variability of n₁, K_s1 and K_s2 (Fig.1)

It can be seen that n ranges from 1.35 to 1.55 [mm^-1], and that common values K_s1=1000 [mm/h] and K_s2=10 [mm/h] are convenient for the whole set of events, when considering the min/max possible values corresponding to the 10 best RMSE. Such common value is not possible for n, but the variability of nseems to be acceptable.

RMSE are found to be 0.034 in average, ranging from 0.020 to 0.047 from the whole set of event.

An example of simulation is given in Fig.2.

The results are similar in the other stations, in the sense that the event variability of the parameters is quite reduced.

2.2. Spatial variability of SHP

The distribution of the parameters have been studied at the transect scale, considering for each station the averaged, min and max values calculated for the set of events of a given station (Fig.3). It can be seen that the averaged values (colored bars) are different for the stations, ranging from 1.34 to 1.74 [mm^-1] for the n parameter, and from 4 to 320 [mm/h] for the K_s2 parameter (parameters range for transect 1). In the case of K_s2, those significant variations could be considered as an effect of the grade of alteration of the bedrock.

As the numeration of the station follows the axis downstream/upstream of each transect, no spatial organization appears to be detected along this axis. The spatial variation thus seems to be randomly distributed at the hillslope scale.

The empirical statistical distributions of the parameters can give an idea of the comparison of the parameters from one transect to another (Fig.4). It can be seen that transect 1, 2, 3 have close distributions for the K_s2 parameter, while transects 2 and 4 are similar for the n parameter. The higher K_s2of T4 could suggest that orthogneiss substratum is less permeable than granite substratum. But both the interpretation and the comparison of the distributions are difficult because of the low number of points of each transect. Those assumptions will have to be further studied with complementary experiment to be confirmed. Beyond the possible physical interpretation, the distribution at the catchment scale shall be used to account for the spatial variability of the soils in modeling applications.

Conclusion

Inverse modeling of soil water content was performed on a dense experimental network consisting in 4 hillslope transects, that is 34 stations and 68 probes. The vadoze area was considered as a 2-layers system, whose layer 1 is soil and layer 2 is weathered area/bedrock. Hydrodynamical parameters were extracted at the event scale, to avoid a complex consideration of evapotranspiration and redistribution of the water content between the most significant events. Three parameters (n for layer 1, K_sfor layer 1 and 2) were finally calibrated at the event scale, the other being previously set as constant in order to avoid equifinality problems.

The estimation of these parameters was proved to be robust enough at the monitoring station scale, for one event to another, as highlighted in an example. Common values could be found for K_s1 and K_s2 for all the events, and the event variability of the nparameter is acceptably reduced. Significant differences of Ks2 value depending on the stations could be representative of the grade of alteration of the bedrock but this hypothesis would required further study to be validated.

At the hillslope scale, no spatial organization could be found along the downstream-upstream axis, and the parameter values appear to be randomly distributed within the hillslope. The comparison of the distribution makes that some differences could be seen for the n₁ and K_s2parameters according a given transect. However, it does not seem possible to conclude to significant differences in the distributions at the transect scale.

As a first attempt of spatial interpretation, this work should be improved to get more accurate interpretation of the spatial variability of the parameter, whether numerical (due to the inverse modeling method) or physical (due to the local heterogeneity of the soils). It is thus planned considering more elaborated assumptions about the inverse modeling protocol (multiple soil layers, multi-objective function, consideration of continuous data) to gain in reliability.

Bibliography

Beven, K.J., 1993. Prophecy, reality and uncertainty in disturbed hydrological modeling. Advances in Water Resources, 16: 41-51.

Bouvier, C. et al., 2007. Recent advances in rainfall-runoff modeling: extrapolation to extreme floods in southern France. In: F.-A. Proceedings of the 1st International Workshop on Hydrological Extremes (Editor), Observing and modeling excpetional floods and rainfalls, Cosenza (It), pp. 229-238.

Cosandey, C., 1994. Formation des crues"cévenoles" dans des bassins élémentaires du Mont-Lozère. Revue des Sciences de l'Eau, 7: 377-393.

Köhne, J.M., Köhnne, S. and Šimůnek, J., 2009. A review of model applications for structured soils: a) Water flow and tracer transport. Journal of Contaminant Hydrology 104: 4-35.

Minasny, B. and Field, D.J., 2005. Estimating soil hydraulic properties and their uncertainty: the use of stochastic simulation in the inverse modeling of the evaporation method. Geoderma, 126: 277-290.

Mualem, Y., 1976. A new model predicting the hydraulic conductivity of unsaturated porous media. Water Resources Research, 12: 513-522.

Pollaco, J.A.P., Soria-Ugalde, J.M., Angulo-Jaramillo, R., Braud, I. and Sauguier, B., 2008. A linking test to reduce the number of hydraulic parameters necessary to simulate groundwater recharge in unsaturated soils. Advances in Water Resources, 31: 355-369.

Ritter, A., Hupet, F., Munoz-Carpena, R., Lambot, S. and Vanclooster, M., 2003. Using inverse methods for estimating soil hydraulic properties from field data as an alternative to direct methods. Agricultural Water Management, 59: 77-96.

Šimůnek, J., Šejna, J.M. and van Genuchten, M.T., 1998. The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media. IGWMC - TPS - 70, International Ground Water Modeling Center, Colorado School of Mines, Golden, Colorado, pp. 186.

van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal, 44: 892-898.

Vrugt, J.A., Stauffer, P.H., Wöhling, T., Robinson, B.A. and Vesselinov, V.V., 2008. Inverse modeling of subsurface flow and transport properties: A review with new developments. Vadose Zone Journal, 7(2): 843-864.

Wöhling, T. and Vrugt, J.A., 2011. Multiresponse multilayer vadose zone model calibration using Markov chain Monte Carlo simulation and field water retention data. Water Resources Research, 47(W04510): 1-19.