Finite element modeling of variably saturated flows in hillslopes with shallow water table

1427

Finite element modeling of variably saturated flows in hillslopes with shallow water table H. Beaugendre a, A. Ern a, J.P. Carlier b, I. Ginzburg b, and C.

Kao b

aCermics-ENPC, 6 et 8 avenue Blaise Pascal, cit6 Descartes, 77455 Marne la Vall~e Cedex 2, France bCEMAGREF, Parc de Tourvoie, B P 44, 92163 Antony Cedex, France During heavy rainfall episodes, subsurface flow can saturate the soil in various regions near the surface and, therefore, contribute to the production of overland flow. This paper investigates numerically the dynamics of water tables in partially saturated porous media and their role in the genesis of surface runoff. The water movement in variably saturated soils is modeled with Richards' equation. The water table position being an unknown of the problem, its intersection with the ground surface yields an unsteady obstacle-type problem. The governing equations are discretized by finite elements in space and an implicit Euler scheme in time. At each time step, the approximate solution is obtained using Newton's method embedded into a fixed-point iteration to determine those points lying on the soil surface where artesian conditions are met. Numerical results are presented for a simple hillslope test case. Particular attention is given to the impact of both the soil water retention curve and the unsaturated hydraulic conductivity function on model results, especially near saturation. 1. I N T R O D U C T I O N

The mechanisms leading to surface runoff in hillslopes exposed to heavy rainfall episodes have received significant attention over the last decades. Historically, the first approach considers streamflow to be generated by runoff over areally extensive portions of a watershed, where rainfall intensity overcomes soil infiltration capacity (Horton [9], 1933). This concept has been improved by the introduction of the partial area contribution approach initiated by Cappus [5] (1960) and Betson [3] (1964) whose work has been further developed by Dunne and Black [8] (1970). In this approach, the actual surface contributing to streamflow is restricted to saturated surfaces (where the water table reaches the ground). This contributing area might represent only a small portion of the watershed, and its extension may vary in time and space. These studies lead to a segmentation of the watershed into infiltration, runoff and exfiltration dominated zones. An important issue is to determine the relative importance of surface and subsurface water in storm hydrographs. Several field and laboratory studies have been conducted (Abdul and Gillham [1], Barros et al. [2]) and revealed that field results are strongly site-dependent, whereas laboratory

1428 tests are complex and time-consuming. The need for numerical experiments has rapidly increased (Ogden and Watts[12], Weiler and McDonnell [15]), allowing faster results and more flexibility. A widely used model to simulate partially saturated flow in porous media is the solution of Richards' equation. The complexity of this problem lies in the strong non-linearity of the equation resulting from the pressure-head-saturation relation and the non-linearity caused by the changing character of the boundary conditions owing to the movement of the water table. These difficulties have been challenged by a great number of researchers who have explored many ways: choice of the main unknown (Celia et al. [6]), of a nonlinear iterative solver, and of the spatial and time discretization scheme. Recently, this type of model has been used to investigate the mechanisms leading to surface runoff in hillslopes (Ogden and Watts [12], Cloke et al. [7]). These studies have emphasized the need for an accurate description of soil hydraulic parameters. Accurate measurement of the soil properties is generally tedious and costly. Consequently, many empirical equations have been used to describe the relations between water content, pressure and hydraulic conductivity (Brooks and Corey [4], Mualem [11], van Genuchten [13]). The choice of one of these models strongly conditions the hillslope hydrology [7], especially near saturation, and then influences the position of the water table. Moreover, the broadly used van Genuchten model may cause numerical instabilities near saturation for fine-textured soils. To overcome this difficulty Vogel et al. [14] proposed to modify the model by introducing a small minimum capillarity height, h~, in the soil water retention curve. One of the goals of this paper is a first assessment of the effect of this modification on the hydrological behavior of hillslopes.

2. G O V E R N I N G E Q U A T I O N S Modelling the soil as an isotropic variably saturated porous medium, we assume that the movement of liquid water can be described with Richards' equation

(1)

OtO(h) - COx. (k(h) . COxp) ,

where t is the time (T), x is the space coordinate (L), 0 is the volumetric water content (dimensionless), k the unsaturated hydraulic conductivity (LT-~), h the soil water pressure head (L), and p - h + z the total hydraulic head (L) where z is the vertical coordinate (upwards). In particular, the use of Richards' equation assumes that there is no air pocket trapped by the flow. The Darcy flow velocity (LT -1) is defined by v(h) -

-k(h)

. Ox~ - - k ( h )

. Oxh - k ( h ) . ez ,

(2)

with ez the upward unit vector. Eq.(1) requires the knowledge of the soil water retention curve, O(h), and the unsaturated hydraulic conductivity function, k ( h ) . Let n > 1 and c~ be the van Genuchten soil parameters, and set m - 1 - 1 I n . Using the additional parameter h~, the modified van Genuchten model (modified VGM model) discussed in [14] is" O~(h) -

~/3 (1 + ( - o z h ) n ) - m ,

h < hs,

t 1,

h _> h~,

(3)

1429 where O'S(h) - O(h)-O~o,_o~ is the effective saturation, 0~ and 0~ being the residual and saturated water content, respectively. The parameter fl is defined as fl - (1 + (-c~h~)~) "~ so that O(h~) - 1. Letting k, be the saturated hydraulic conductivity, we set k(h) - k~k~(h) where k~ is the relative hydraulic conductivity function. Substituting Eq.(3) into the Mualem model [11] for k~ gives:

1-- (1--(/~-10)l/m) m] 2 [l_(l_/3_1/m)m]2

(05 -

,

(4)

h < hs ,

h>h~.

1,

The parameter h~ < 0 is interpreted, in [14], as the minimum capillary height. Equations (3) and (4) reduce to the original VGM model when h~ - 0 (/3 - 1). Recall from [14] that, in the original VGM model, n < 2 guarantees continuity up to first order derivative in O(h) at saturation, while n >_ 2 ensures continuity up to second order derivative in O(h) at saturation and, hence first order continuity in the water capacity function C(h) - dh" dO With the modified VGM model for n < 2, the slopes k'~(h) and C'(h) remain finite but discontinuous when h reaches h~ < 0, as opposed to becoming infinite when h --+ 0 with the original VGM model. Note also that this modification leads to a loss of continuity of C~(h) even when n > 2. The soil parameters chosen for this paper and corresponding to SCL, YLC and SAND soils, are summarized in Table 1.The relative hydraulic conductivity functions and water capacity, C(h), as described with the original VGM model are plotted in Figure 1.

Table 1 Soil parameters [10], k~" saturated hydraulic conductivity; 0~" saturation water volume fraction; 0~" residual water volume fraction; c~ and n: van Genuchten soil parameters.

Soils

k~ (m/h)

c~ (Tft-1)

Silty Clay Loam (SCL)

0.0026

1.9

Yolo Light Clay (YLC)

0.018

Sandy Soil (SAND)

0.1

II (-)

Os (-)

Or (-)

1.31

0.41

0.10

3.6

1.9

0.55

0.23

3.7

5.0

0.50

0.05

3. T E S T P R O B L E M

The two-dimensional domain, f~ (its boundary Of~), selected to accomplish the study is sketched in Figure 2 and corresponds to the configuration considered by Abdul and Gillham [1]. The domain dimensions are 1.4 m length and 1 m to 0.8 m height. We assume symmetry condition on the left and right boundaries and impermeability at the bottom. Letting n be the normal vector outwards to the surface, we impose v . n = 0 at the upslope end (left boundary), at the downslope end (right boundary), and at the bottom surface of the domain. A constant rainfall intensity i (LT -~) is considered; the rainflow velocity is defined as v~ = - i e z where ez is the upward unit vector so that

1430

1.2

0.8

SCL YLC ............... -SAND .............

....................";,1

/I

7 ,7.

0.8

0.6 E

_~

SCL YLC ............... SAND .............

0.7

0.6

,,i-v

o 0.4

0.5 0.4 0.3 0.2

0.2 0 -35

0.1 I

i

i

i

i

i

-30

-25

-20

-15

-10

-5

0

t

0

5

6 -14 -12 -10

h (cm)

-8

-6

-4

-2

0

h (cm)

Figure 1. Hydraulic conductivity functions, k~(h) (left), and water capacity, for SCL, YLC and SAND soils using the original V G M model.

C(h)

(right),

Figure 2. Abdul and Gillham test problem.

i > 0. The upper surface of the model domain (top boundary: Ot2t) allows infiltration and exfiltration. The top boundary may be divided into two regions 9 0t2 + the non-saturated region" h _< 0 and v - ~ -

v~.

9 0t2~- the saturated region" h - 0 and v . r~ _> v~. The saturated region may also be divided into two subregions, see left panel of Figure 3" 9 the saturated zone that still allows some infiltration" h - 0 and v<. r~ _< v. n _< 0 9 the exfiltration region" h - 0 and v-r~ > 0. The feedback of surface runoff on the water table dynamics is neglected; in some sense, one assumes that surface runoff is instantaneously evacuated from the system.

1431

Figure 3. Left" scheme of the infiltration, partial infiltration and exfiltration lengths. Right" typical pressure head, h (m), and normal velocity, v . n (m/h), distribution along the top surface, obtained with the original VGM model (at saturation h - 0).

4. N U M E R I C A L

SCHEMES

The pressure head h is chosen as the main unknown of the problem. Richards' equation is discretized using a finite element method (p1 conforming) on unstructured, non-uniform triangulations. The water table position being an unknown of the problem, its intersection with the ground surface is treated as an unsteady obstacle-type problem. 4.1. S t e a d y p r o b l e m A steady-state solution exists whenever i < k~. For the sake of simplicity, we present the algorithm assuming that the water table intersects the top boundary at one point z0 at most. It implies that the region {z > z0} corresponds to the non-saturated zone. The weak formulation corresponding to the steady problem is the following: given z0 C R, let

l/zo -- { 0 E H l(f~); 0 - - 0 oil (;gf~tN {z < z0} },

(5)

and the form (non-linear in h, linear in O) azo (h,

f -

f

! k(h)(o h +

+!

J f~

(6)

Jo f~tn{z>zo}

Note that working with the space V~o implies h - 0 on the saturated area cgt2t N {z < z0}. Recall that H 1 is the space of square integrable functions whose distributional derivative is also square integrable. Then, we seek z0 E R and h c Vzo such that

(~)

azo (h, r - 0

(ii)

h<0

onOt2tn{z>zo},

(iii)

v(h).n

> v~.n

vr c ~ o , (7)

on cgt2t n {z < z0}.

Here zo represents the height of the water table. Note that the well posedness of (i) requires that Oftt N {z < z0} -r ~, i.e. that the water table has reached the top boundary.

1432 An approximate solution is sought using Newton's method embedded into a fixed-point iteration to determine those points lying on the soil surface where artesian conditions are met. The following iterative algorithm is proposed to solve the problem: 1. choose an initial z0 2. construct a triangulation T; 3. solve problem (i) using e.g. P~ conforming finite elements method; 4. check whether (ii)and (iii)are satisfied; 5. if (ii) is satisfied and (iii) is not, move z0 one mesh cell down; go to step 3; 6. if (iii) is satisfied and (ii) is not, move z0 one mesh cell up; go to step 3;

7. if both (ii) and (iii) are satisfied, then current pair (z0, one may refine the mesh and go back to step 3.

h(zo))

is the desired solution;

Note that from the maximum principle, both (ii) and (iii) can not be violated simultaneously. However, in numerical approximations, this may happen. In this case, we still consider that the water table has been correctly positioned. With this "loosened" convergence criterion, the final position of the water table depends from whether the converged position z0 has been approached from below or above. The two resulting values give lower and upper bounds for the water table position (typically differing from one or two mesh cells at the most).

4.2. U n s t e a d y problem The unsteady problem is solved using an implicit Euler scheme. For a time step k _> 0, given Zok and h k, we seek z0k+l and h k+l E Vzg+~ such that (i)

1 ff't (0( h k+l ) -- O (h) )k 5---t

(ii)

h k+l < 0

(iii)

v(hk+~).n >_v~.n

r Jr- az~)+l (h k+ 1 , r

-- 0

(8)

on 0 f ~ t N { z > z0k+l},

on

cggtt N {z

< z0k+l}.

The problem is solved using the same iterative algorithm as for steady problem. In step 1, the initial choice is z0k+l - z0k. Note also that in the unsteady case, problem (i) is well posed even if the water table has not reached the top boundary.

4.3. C o m p u t a t i o n of the infiltration and exfiltration fluxes To evaluate the infiltration and exfiltration fluxes associated with the discrete solution at a given time, we use the following approach based on test functions (this method is in fact the natural way of expressing mass conservation in the framework of the approximate scheme). The top boundary, O~t, is divided into two regions, an exfiltration region, cgf~x, where v . n > 0, and an infiltration region, cgf~i~, where v . n _< 0. Let 9~1 d and 9~d be the continuous piecewise linear functions defined on 0f~t such that pd _ 1 (respectively 9&d -- 1) at all the mesh nodes on Of~x (respectively on O f ~ ) . Then, there exist r d, Csd continuous piecewise linear functions defined on the whole computational domain, f~,

1433 such that" Od laa~- P~, r 10~,-+ infiltration flux Qd~ (t k) and the exfiltration flux -

V Ol

- 1 on For a time step k >_ i, the Qdex(tk) are defined as follows" 5t

(9)

5t

Clearly

~d(tk) + Qid(tk) --

5. N U M E R I C A L

f~ ok_ok--1 ~t ~

0 when k ~ +oc.

RESULTS

Numerical results are presented for the configuration shown in Figure 2 and for the three soils characterized in Table I. First we assess the level of grid refinement; results are presented for the YLC soil and the original VGM model. Then, we investigate the influence of rainfall intensity and the influence of the soil hydraulic parameters on the hillslope response still using the original VGM model for the three soils. Finally, we study the influence of the minimum capillary height on model results. 5.1. N u m e r i c a l a s s e s s m e n t We introduce L the length of the slope and Ls the portion of the hillslope that is saturated at steady-state. The time to reach equilibrium Te is defined as the time for the discrete water table to reach its equilibrium position; Te is obviously mesh-dependent. We consider the YLC soil subjected to a constant rainfall intensity i/k~-lO%. Here and below, the initial condition is a horizontal water table located at 0.7 m. The simulation is conducted through steady-state on two meshes to investigate the solution dependence to the mesh. Mesh 1 is obtained with an unstructured triangulation containing 4083 nodes and 7903 elements: the node spacing on the top surface is 1 cm and 5 cm at the bottom surface. Mesh 2 is obtained by imposing nodes spacing equal to 0.5 cm at the equilibrium water table position that results from the computation on Mesh 1. Mesh 2 contains 7429 nodes and 14474 elements. The agreement between the results produced by the two meshes is very satisfactory with relative error less than 1%; see Table 2.

Table 2 Mesh influence on the steady state

Te

(h)

Ls/L (%)

Mesh 1

4.18

60.28

Mesh 2

4.15

60.14

5.2. I n f l u e n c e of rainfall i n t e n s i t y on t h e h i l l s l o p e r e s p o n s e We consider the YLC soil subjected to four different rainfall intensities i/k~-2, 5, 10 and 15% respectively. The simulations are conducted through steady-state. Figure 4 a) plots the evolution in time of Ls/L for the four intensities. As expected, the relative hillslope saturated area increases with the rainfall intensity. This figure illustrates the

1434 nonlinear temporal response of this soil under different rainfall intensities, in agreement with Ogden and Watts [12]. This nonlinearity can be related to the soil hydraulic functions combined with the effect of the ratio i/k~. The higher this ratio, the faster the pressure head reaches small absolute values; note that for the YLC soil, k~ increases rapidly with pressure head near saturation. Letting Qrain - i , L be the rainfall rate, Figure 4 b) plots the temporal evolution of the relative exfiltration. Observe that this relative exfiltration decreases non uniformly in time as rainfall intensity increases; this can be related to a higher water table position at equilibrium.

100

1.2

i/ks=2% i/ks=5% ............ i/ks= 10% .............. i/ks= 15 % ............................

80

/

v

,~." .......".." ...../[ ."

J

J

0.8

/ .. .......................................

6O oO

Exfiltration i/ks=2% Exfiltration i/ks=5% ............ Exfiltration i/ks= 10% Exflitration i/ks=l 5% ............................

40

0

. . . ./. . iiii. ............................................

20

. . . . . . . iiiiiiiiiiiiill

0.6

0 0.4 0.2

0 0

' 0.2

0 0.4

0.6

0.8 t/Te

1

1.2

1.4

0

0.2

0.4

0.6

0.8

1

1.2

1.4

t/Te

b

Figure 4. a) Temporal evolution of the relative hillslope saturated area for three different rainfall intensities with YLC soil. b) Temporal evolution of the exfiltration for three different rainfall intensities with YLC soil.

5.3. I n f l u e n c e of t h e soil h y d r a u l i c p a r a m e t e r s on t h e d y n a m i c s of t h e w a t e r table 5.3.1. o r i g i n a l V G M m o d e l The evolution in time of the water table location has been studied for the three soils: SCL, YLC and SAND. Simulations were run with a constant rainfall intensity, i/k~ - 10 %, for a duration longer than the time to reach equilibrium for that rain rate. Table 3 gives the time to equilibrium for the three soils. The controlling parameter for Te appears to be k~ which varies in a very wide range. The rainfall intensity i increases as k~ does, and less time is then needed to fill the domain. Other parameters (n, c~, 0 ~ - 0~) will influence Te but to a smaller extent in our case. The SAND soil exhibits the fastest response and the SCL soil the slowest. Plots in Figure 5 a) of the fraction of hillslope saturated area, L s / L , versus t T e -1, reveal the temporal response of the hillslope for the three different sets of soil parameters. The shape of the three curves is a consequence of the respective shapes of the hydraulic functions for the three soils; see Figure 1. At steady-state, Figure 5 b) shows strongly different distributions of 0 on the unsaturated surface O f F , while L s / L is identical for the

1435 three soils as it only depends on the rainfall intensity i/k~ for mass conservation reasons. More precisely, the value of the parameter n conditions the curvature of the plots. For instance, a sole modification of this parameter set to 5 for the YLC soil gives a similar response to that of the SAND soil.

Table 3 Time to equilibrium, Te, for the three soils using the original VGM model SCL

Te (h)

10.89

YLC 4.18

SAND 0.58

The evolution of the exfiltration fraction of the hillslope, Lezf/L, versus tTe -1, for the three soils is similar to Figure 5 a). We observe that only a part of the saturated area contributes to the exfiltration. At steady-state this part is independent of the soil and represents about 70 % of the saturated length. Figure 5 c) shows the fraction of the exfiltration and runoff versus tTe -1 for the three soils. Runoff is composed of exfiltration and of "direct" runoff which represents the water that never infiltrates. The dynamic of the exfiltration is similar to the evolution of the water table position for each soil. At equilibrium the exfiltration represents 45.9 % of runoff and direct runoff represents 54.1 % for each soil. Figure 5 d) represents the temporal evolution of the ratio exfiltration/runoff for the three soils. From a global point of view, the value of this ratio appears roughly identical for all soils: as runoff begins, it immediately reaches some value between 40 and 50 % and does not vary much till steady-state (t=Te). The differences between the curves essentially lie on an early occurrence of runoff (and hence exfiltration) for coarser textured soils. To a lesser extent, fine textured soils reach a maximum value of the ratio much faster after the beginning of the runoff and then slowly decreases to the steady-state value, while the SAND soil's exfiltration ratio slowly increases to a maximum value reached just before equilibrium. These behaviors are linked to the nature of the hydraulic functions. It is then quite remarkable that for all soils, exfiltration occurs as soon as runoff does, and that the ratio comes immediately to a nearly constant value, which means that at any time water exiting the system will have roughly the same composition. 5.3.2. modified VGM model: h~--2 cm The effect of the modified VGM model on the relative hydraulic conductivity function, k~, and water capacity for each soil is shown on Figure 6. The same value h~ = - 2 cm has been chosen for the three soils. The modification affects more significantly the soils with n close to 1, whereas for the SAND soil the curves are very similar. Simulations were run with a constant rainfall intensity, i/k~ - 10 %, for a duration longer than the time to reach equilibrium for that rain rate. The influence of the modified VGM model, with h~ = - 2 cm, on the evolution of the water table is shown on Figure 7. The modified VGM model does not affect the position of the water table at equilibrium, but impacts its evolution in time: saturation appears earlier and, hence, exfiltration and runoff occur faster; see Figure 8. This influence is more noticeable when n is small, since

1436

100

1

S C L soil Y L C soil ............... S A N D soil .............

80

0.998 0.996 0.994

60 v

0.992

.,J

j

0.99

40

0.988 0.986 0.984

20 .... o.

0 0

0.2

0.4

0.6

0.8

1

i

i

1.2

1.4

S C L soil Y L C soil ............... S A N D soil .............

0.982 0

0.2

0.4

0.6

t/Te

Ex{iltration SCL Exfiltration YLC Exfiltration SAND Runoff SCL Runoff YLC R u n off S A N D

1.2

0.8

. . . . . ............. ........... ............. .................................. .............. ..." ..: ...............-:"".../ ::"

0.6

..

.::

/

yd'

:.; ..,

0.8

1

1.2

1.4

x(m)

;.

SCL YLC SAND

80

............... .............

o~

v

o

60

e-

...-

9

100

...,

rr

0.4

,-

" - "

40

O

....

0.2

20

;J

0

0

0.2

0.4

0.6

0.8 t/Te

1

1.2

1.4

L

0

0.2

0.4

0.6

0.8

1

J

1.2

1.4

t/Te

d

Figure 5. a) Temporal evolution of the relative hillslope saturated area for the three soils, original VGM model, b) Distribution of 0 along the top surface for the three soils, c) Temporal evolution of the relative outflow and runoff for the three soils, original VGM model: SCL soil, YLC soil and SAND soil. d) Temporal evolution of the ratio, Exfiltration/Runoff, for the three soils SCL, YLC and SAND using original VGM model.

the influence of h~ on the hydraulic conductivity functions is more important. 6. C O N C L U S I O N Richards' equation combined with the original [13] and modified [14] van Genuchten's model for hydraulic functions has been discretized by finite elements in space. The changing character of the surface boundary conditions owing to the movement of the water table yields an unsteady obstacle-type problem for which an algorithm is proposed. This model has been used to investigate the hydraulic behavior of hillslopes under constant rainfall conditions. Numerical results have first shown that rainfall intensity has a large influence on the length of the saturated hillslope and on the exfiltration ratio which decreases as i increases.

1437

1.2

0.8

modified VGM - SCL modified VGM - YLC ............. modified VGM - SAND ........... ,...............,

modified VGM - SCL modified VGM - YLC modified VGM - SAND

0.7 0.6

0.8

,,-,..

\

0.5

E 0.6

v

o 0.4

............... .............

0.4 0.3 0.2

0.2

0.1 ,

0

";;"2 " ; " ~ " " ~ ": .....

-35

-30

,

,-""~

......... , ......... , ....

"]

0

i

J

i

i

i

i

-25

-20

-15

-10

-5

0

_.

5

6 -14

-12

-10

h (cm)

-8

-6

-4

-2

0

2

h (cm)

Figure 6. Hydraulic conductivity functions, k~(h) (left), and water capacity, C(h) (right), for the modified (h~ - - 2 cm) VGM models on the SCL, YLC and SAND soils.

100

100

original VGM SCL modified VGM SCL .............

80

original VGM YLC modified VGM YLC .............

80

60

60

v

v

J

40

40

20

20 /

.,,,,/ 2

4

6

8

10

12

1

14

2

t(h)

3

4

6

7

t(h)

100

original VGM SAND modified VGM SAND ........

80 60 v _.1 or)

._j

40 20

0

0.2

0.4

0.6

0.8

1

t(h)

Figure 7. Effect of the modified VGM model (h~ - - 2 cm) on the temporal evolution of the relative hillslope saturated area for the three soils: SCL soil (top left); YLC soil (top right) and SAND soil (bottom).

For a fixed value of the ratio i/k~, the structure of the soil, characterized by the van Genuchten hydraulic parameters, strongly affects the dynamic response of the system,

1438

1.4

SCI': Exfiltr&tion original SCL: Runoff original SCL: Exfiltration modified SCL: Runoff modified

1.2

0.8

1.4

VGM ' VGM .......... VGM ............ VGM ............

...";;; .'ii[,"

0.8

.: //' ;." .." // /,

0.6

YLC: Exfiltration driginal VGM YLC: Runoff original .......... YLC: Exfiltration modified VGM ............ YLC: Runoff modified VGM ............

1.2

....;::."

0.6

...;)" 0.4

0.4

0.2

0.2

. . r 0

2

4

6

8

10

12

14

0

1

1.4

,

i

4

5

SAND! Exfiltration'original VGM SAND: Runoff original . . . . . SAND: Exfiltration modified VGM , SAND: Runoff modified VGM .........

1.2

a3

3 t(h)

t(h)

._.G

2

0.8 0.6 0.4 0.2

0

0.2

0.4

0.6

0.8

t(h)

Figure 8. Effect of the modified VGM model (h~ - - 2 cm) on the temporal evolution of exfiltration and runoff for the three soils" SCL soil (top left); YLC soil (top right) and SAND soil (bottom).

whereas the steady-state values remain identical. Results similar to those of Vogel et al. [14] for infiltration in very dry soils are observed while using the modified VGM model, especially for fine textured soils, where we have found that the introduction of a minimum capillary height h~ has a significant influence on the runoff and exfiltration genesis. The driving parameter for the shape of the curves plotting the evolution of saturation and exfiltration ratio is suspected to be the exponent n. Further investigations are being performed to determine whether its influence expresses itself through the function k~(h) or C(h) or both. From the basis of this primary findings, a sensitivity analysis of this type of problem should be conducted as a next step, in order to have a better understanding of the parameters conditioning surface runoff in hillslopes. ACKN OWL ED G EMENT S

This work has been supported by INRIA under the Cooperative Research Project DYNAS (http'//www-rocq.inria.fr/estime/DYNAS). We thank our colleagues from CEREVE-

1439 ENPC and INRIA for stimulating discussions. REFERENCES

1. A.S. Abdul and R.W. Gillham, Laboratory studies of the effects of the capillary fringe on streamflow generation, Water Resour. Res. No. 20 691-708 (1984). 2. A.P. Barros, D. Knapton, M.C. Wang and C.Y. Kuo, Runoff in shallow soils under laboratory conditions, J. Hydraul. Eng. No. 4(1) 28-37 (1999). 3. R.P. Betson, What is watershed runoff?, J. Geophys. Res. No. 69(8) 1541-1551 (1964). 4. R.H. Brooks and A.T. Corey, Hydraulic properties of porous media, Hydrol. Pap. 3, Colo. State Univ., Fort Collins, pp. 27 (1964). 5. P. Cappus, Etude des lois de l'6coulement. Application au calcul et ~ la prdvision de d6bits. Bassin exp6rimental de l'Alrance. La Houille Blanche A 493-520 (July-August 1960). 6. M.A. Celia, E.T. Boulotas and R.L. Z~rba, A general mass-conservative numerical solution for the unsaturated flow equation, Water Resour. Res. No. 26(7) 1483-1496 (1990). 7. H.L. Cloke, J.-P. Renaud, A.J. Claxton, J.J. McDonnell, M.G. Anderson, J.R. Blake and P.D. Bates, The effect of model configuration on modelled hillslope-riparian interactions, J. of Hydrology No. 279 167-181 (2003). 8. T. Dunne and R.D. Black, Partial area contributions to storm runoff in a small New England watershed, Water Resour. Res. No. 6 1296-1311 (1970). 9. R.E. Horton, The role of infiltration in the hydrologic cycle, Eos Trans. AGU No. 14 446-460 (1933). 10. C. Kao, S. Bouarfa and D. Zimmer, Steady state analysis of unsaturated flow above a shallow water-table aquifer drained by ditches, J. Hydrol. No. 250 122-133 (2001). 11. Y. Mualem, A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. No. 12 513-522 (1976). 12. F.L. Ogden and B.A. Watts, Saturated area formation on non convergent hillslope topography with shallow soils : a numerical investigation, Water Resour. Res. No. 36(7) 1795-1804 (2000). 13. M.Th. van Genuchten, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. No. 44 892-898 (1980). 14. T. Vogel, M. Th. van Genuchten and M. Cislerova, Effect of the shape of the soil hydraulic functions near saturation on variably-saturated flow predictions, Advances in Water Resources No.24(6) 133-144 (2001). 15. M. Weiler and J. McDonnell, Virtual experiments: a new approach for improving process conceptualization in hillslope hydrology, J. of Hydrology No. 285 3-18 (2004).