# A seepage face model for the interaction of shallow water tables with the ground surface: Application of the obstacle-type method

Héloïse Beaugendre, Alexandre Ern, Thomas Esclaffer, Eric Gaume, Irina Ginzburg, Cyril Kao

<figure><img src="image_1.png" /></figure>

## To cite this version:

Héloïse Beaugendre, Alexandre Ern, Thomas Esclaffer, Eric Gaume, Irina Ginzburg, et al. A seepage face model for the interaction of shallow water tables with the ground surface: Application of the obstacle-type method. Journal of Hydrology, 2006, 329 (1-2), pp.258-273. 10.1016/j.jhydrol.2006.02.019. inria-00337535

HAL Id: inria-00337535

https://inria.hal.science/inria-00337535v1

Submitted on 14 Nov 2008

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
A seepage face model for the interaction of
shallow water tables with the ground surface:
Application of the obstacle-type method

H. Beaugendre<sup>a</sup>, A. Ern<sup>b</sup>, T. Esclaffer<sup>c</sup>, E. Gaume<sup>c</sup>,
I. Ginzburg<sup>d</sup>, C. Kao<sup>d</sup>

<sup>a</sup> INRIA Rocquencourt, BP 105, 78153 Le Chesnay Cedex, France

Currently at: MAB, Université Bordeaux 1, 351 cours de la Libération, 33405

Talence Cedex, France

<sup>b</sup> CERMICS, Ecole Nationale des Ponts et Chaussées, 6 et 8 avenue Blaise Pascal,
77455 Marne la Vallée Cedex 2, France

<sup>c</sup> CEREVE, Ecole Nationale des Ponts et Chaussées, 6 et 8 avenue Blaise Pascal,
77455 Marne la Vallée Cedex 2, France

<sup>d</sup> CEMAGREF, Parc de Tourvoie, BP 44, 92163 Antony Cedex, France

**Abstract**

This paper presents a model to simulate overland flow genesis induced by shallow water table movements in hillslopes. Variably saturated subsurface flows are governed by the Richards equation discretized by continuous finite elements on unstructured meshes. An obstacle-type formulation is used to determine where saturation conditions, and thus seepage face conditions, are met at the ground surface. The impact of hillslope geometry, boundary conditions, and soil hydraulic parameters on model predictions is investigated on two-dimensional test cases at the metric and
hectometric scales. The obstacle-type formulation is also compared with a more detailed model coupling subsurface and overland flow, the latter being described by the shallow water equations in the diffusive wave regime.

***Key words:*** Finite elements, Richards equation, Seepage face, Obstacle-type model, Water table ground interaction, Numerical solutions.

# 1 Introduction

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 overland flow to be generated by runoff over those portions of the watershed where the rainfall intensity exceeds the soil infiltration capacity; see Horton (1933). This concept has been refined by the introduction of the partial area contribution concept initiated by Cappus (1960) and Betson (1964), whose work has been further developed by Dunne and Black (1970). In this approach, the actual surface contributing to overland flow is restricted to saturated surfaces, *i.e.*, where the water table reaches the ground level. This leads to a segmentation of the watershed into infiltration, exfiltration and runoff dominated zones. Since the contributing area represents only some portion of the watershed, a first important issue in hillslope hydrology is to determine the spatial extension of the contributing area and its temporal evolution resulting from the dynamics of subsurface flows. A further important issue is to determine the relative importance of surface and subsurface water in storm hydrographs and to compare these different contributions with results obtained by flow hydrograph separations (Loye-Pilot and Jusserand (1990); Iorgulescu (1997); Marc et al. (2001); McGlynn and
McDonnell (2003)).

To tackle these issues, several field and laboratory studies have been conducted (Abdul and Gillham (1984); Barros *et al.* (1999)). However, field experiments are often strongly site-dependent, and laboratory tests are often complex and time-consuming. This has prompted the development of numerical tools offering more flexibility and faster results in the simulation of hillslopes subjected to various operating conditions (Ogden and Watts (2000); Cloke *et al.* (2003); Weiler and McDonnell (2004)). These so-called virtual hillslope experiments have emerged as a valuable complement to field and laboratory measurements for improving process conceptualization in hillslope hydrology.

Variably saturated flows in porous media can be modeled with the Richards equation (Richards (1931)). This equation is based on two assumptions that, in some circumstances, can preclude its use. On the one hand, that the gas-phase remains connected so that a single air pressure can be defined. On the other hand, that there exist an algebraic relation linking the volumetric water content to the hydraulic head and that the water flow velocity can be expressed using a generalized Darcy law involving a relative hydraulic conductivity. Despite these limitations, the Richards equation is now a well-established approach to simulate water table dynamics; see, for instance, Celia *et al.* (1990) for a thorough discussion of numerical aspects including choice of the main unknown, of the nonlinear iterative solver and of the space and time discretization schemes. Using the Richards equation presents two advantages, namely the use of a single (nonlinear) partial differential equation that can be discretized by fairly standard or more advanced finite element or finite volume techniques (Woodward and Dawson (2000); Knabner and Schneid (2002); Basarian (2003); Bause and Knaber (2004); Fagherazzi *et al.* (2004); Manzini and
Ferraris (2004)) and the fact that the saturated and unsaturated portions of
the soil can be treated simultaneously as a single computational domain. Mod-
els based on the Richards equation can be used to simulate the response of
the water table to infiltration caused by rainfall.

A difficulty appears when the water table rises so much that some portion
of the ground is saturated. In this case, the intersection of the water table
with the ground surface produces seepage faces where infiltration, exfiltration,
overland flow and re-infiltration processes may occur, Esteves et al. (2000).
Simpler approaches neglecting flood routing and re-infiltration processes can
still improve the current understanding of hillslope hydrology; see, for instance,
Ogden and Watts (2000); Cloke et al. (2003, 2006) where the mechanisms lead-
ing to surface runoff in hillslopes are investigated. In particular, these studies
have emphasized the need for an accurate description of the seepage face area
(Rulon et al. (1985)). A mathematical way to determine the extent of the
seepage face is to consider that on the ground surface the boundary condition
for the Richards equation changes nature, e.g., a Dirichlet condition imposing
the head (saturated area) or a Neumann condition specifying the water flux
(unsaturated area). Since the water table position is a priori unknown, its in-
tersection with the ground surface leads to an unsteady obstacle-type problem,
in which the partition of the ground surface into saturated and unsaturated
zones and the water table dynamics are solved simultaneously.

The goal of this paper is threefold. First, to formalize the mathematical setting
for the obstacle-type model that can be used to describe the water table dy-
namics in a variably saturated soil where saturation conditions are met in an
*a priori* unknown part of the ground surface. Second, to elaborate on previous
work (Ogden and Watts (2000); Cloke *et al.* (2003)) to further investigate the
impact of soil hydraulic parameters, hillslope geometry, and boundary condi-
tions on model predictions, and thus to assess the suitability of the simulation
tool to carry on virtual hillslope experiments accurately and efficiently. Third,
to compare the obstacle-based approach to a more detailed model where the
overland flow, which is assumed to be governed by the shallow water equa-
tions in the diffusive wave regime, is coupled to the Richards equation for the
subsurface flow through the boundary conditions. Similar situations where
these couplings occur between the shallow water equations and the Navier-
Stokes equations or the Darcy equations have been recently investigated in
Discacciati et al. (2002); Miglio et al. (2003).

This paper is organized as follows. Section 2 presents the governing equations in the framework of the obstacle-type model and describes the numerical methods. Section 3 discusses the results obtained with the obstacle-type formulation on two-dimensional test cases at the metric and hectometric scales. The impact of hillslope geometry, boundary conditions, and soil hydraulic parameters on model predictions is studied. Section 4 investigates the more detailed model where the shallow water type equation describing the propagation of the overland flow is coupled with the water table dynamics. Section 5 draws the conclusions.

2 The Obstacle-Type Model (OTM)

The goal of this section is to describe the OTM and its numerical approxima-
tion. The OTM is based on the Richards equation to predict the water table
dynamics and on the assumption that the height of the overland flow as well
as re-infiltration processes have a negligible impact on the subsurface flow.
The numerical approximation is based on a finite element discretization and a Newton scheme embedded into a fixed-point iteration to determine those points at the ground where saturation conditions are met.

## 2.1 The Richards equation

The soil is modeled as an isotropic, variably saturated, porous medium in which air entrapment does not affect the flow. The movement of liquid water can be described by the Richards equation

$$ \frac{\partial}{\partial t}\theta(\psi) + \nabla \cdot v(\psi) = 0, \qquad (1) $$

where $\theta$ is the volumetric water content (dimensionless), $\psi$ the soil water pressure head (L) and $v(\psi)$ the Darcy flow velocity (LT$^{-1}$) defined by

$$ v(\psi) = -k(\psi)\nabla H, \qquad (2) $$

where $k$ is the hydraulic conductivity (LT$^{-1}$), $H = \psi + z$ the total hydraulic head (L) and $z$ the vertical coordinate (upwards). The Richards equation can be rewritten as

$$ \frac{\partial}{\partial t}\theta(\psi) = \nabla \cdot (k(\psi)\nabla\psi + k(\psi)e_z), \qquad (3) $$

where $e_z$ is the upward unit vector. The hydraulic conductivity is written in the form $k(\psi) = k_s k_r(\psi)$ where $k_s$ is the hydraulic conductivity at saturation and $k_r(\psi)$ is the relative hydraulic conductivity.

Eq. (3) requires the knowledge of the soil water retention curve, $\psi \mapsto \theta(\psi)$, and the relative hydraulic conductivity function, $\psi \mapsto k_r(\psi)$. Accurate mea-
surement of the soil properties is generally tedious and costly. Consequently, more or less empirical equations have been used to describe the relations between water content, pressure and hydraulic conductivity (Brooks and Corey (1964); Mualem (1976); van Genuchten (1980)). The choice of one of these models can have a sizeable impact on the flow dynamics especially near saturation and, thus, substantially modify the position of the water table (Cloke et al. (2003)). Virtual hillslope experiments must, therefore, be carefully assessed as far as their sensitivity to model parameters is concerned. Moreover, the broadly used van Genuchten Mualem model can cause numerical instabilities near saturation for fine-textured soils. To overcome this difficulty, Vogel et al. (2001) proposed to introduce a minimum capillarity height, $\psi_s$, in the soil water retention curve. In this work, the van Genuchten Mualem model (VGM model) and its modification proposed in Vogel et al. (2001) are considered.

## 2.2 Boundary conditions

The pressure head $\psi$ is chosen as the main unknown in the Richards equation. To specify suitable boundary conditions, we consider the two-dimensional domain $\Omega$ with boundary $\partial\Omega$ sketched on the left panel of Figure 1. Let $n$ be the outward normal vector. A constant rainfall intensity $i$ (LT$^{-1}$) is considered; the rainfall velocity is $v_r = -ie_z$.

The upper surface of the computational domain, $\partial\Omega_t$, allows for infiltration and exfiltration. The top boundary $\partial\Omega_t$ is split into $\partial\Omega_t^- \cup \partial\Omega_t^+$ where $\partial\Omega_t^-$ corresponds to the non-saturated surface and $\partial\Omega_t^+$ to the saturated surface. The key assumption made in the OTM is that the height of the overland flow as well as re-infiltration processes have a negligible impact on the water
table dynamics. As a result, the pressure head can be set to the atmospheric pressure on $\partial\Omega_t^+$, i.e $\psi = 0$. Following our notation, the flux through the top surface, $v(\psi) \cdot n$, is locally negative if infiltration occurs or positive if exfiltration occurs. The rainfall flux, $v_r \cdot n$, is then negative. On the unsaturated area, the infiltration flux is equal to the rainfall flux, and in the saturated area the infiltration flux is smaller than or equal to the rainfall flux.

Hence, the following holds:

$$ \psi \le 0 \text{ and } v(\psi) \cdot n = v_r \cdot n \quad \text{on } \partial\Omega_t^-, \qquad (4) $$

$$ \psi = 0 \text{ and } v(\psi) \cdot n \ge v_r \cdot n \quad \text{on } \partial\Omega_t^+. $$

The saturated surface $\partial\Omega_t^+$ can be further divided into:

* the partial infiltration zone in which $v_r \cdot n \le v(\psi) \cdot n \le 0$;

* the exfiltration zone in which $v(\psi) \cdot n > 0$.

This splitting is illustrated in the right panel of Figure 1. The saturated zone $\partial\Omega_t^+$ extends from point B down to the toe of the slope (point D). The exfiltration zone extends from point C down to D.

For the bottom and left boundaries, an impermeable layer is assumed, i.e., the condition $v(\psi) \cdot n = 0$ is imposed. On the right boundary, $\partial\Omega_r$, two boundary conditions can be considered: either a no-flow condition representing a symmetric configuration at the toe of the slope (henceforth referred to as BC1) or a specified total hydraulic head accounting for the presence of a stream (henceforth referred to as BC2).
## 2.3 Numerical schemes for the OTM

Since the water table position is a priori unknown, its intersection with the ground surface leads to an unsteady obstacle-type problem in which the partition of $\partial\Omega_t$ into $\partial\Omega_t^- \cup \partial\Omega_t^+$ and the water table dynamics must be solved simultaneously. Early work based on a similar approach includes that of Neumann (1973) and Rulon et al. (1985).

### 2.3.1 Stationary OTM

To formalize the OTM, consider first a stationary setting and, for the sake of simplicity, assume a BC1 boundary condition on $\partial\Omega_r$. Recall that $H^1(\Omega)$ is the space of square integrable functions in $\Omega$ whose distributional derivative is also square integrable in $\Omega$. The weak formulation corresponding to the steady problem is the following: given $\partial\Omega_t^+ \subset \partial\Omega_t$, let

$$V_{\partial\Omega_t^+} = \{ \phi \in H^1(\Omega); \phi = 0 \text{ on } \partial\Omega_t^+ \}, \qquad (5)$$

and let $a$ be the form (non-linear in $\psi$, linear in $\phi$)

$$a_{\partial\Omega_t^+}(\psi, \phi) = \int_{\Omega} k(\psi) \cdot \left(\frac{\partial}{\partial x}\psi + e_z\right) \cdot \frac{\partial}{\partial x}\phi + \int_{\partial\Omega_t^-} (v_r \cdot n) \phi. \qquad (6)$$
Then, we seek $\partial\Omega_t^+ \subset \partial\Omega_t$ and $\psi \in V_{\partial\Omega_t^+}$ such that

$$ \left\{ \begin{array}{ll} (i) & a_{\partial\Omega_t^+}(\psi, \phi) = 0 \quad \forall \phi \in V_{\partial\Omega_t^+}, \\ (ii) & \psi \le 0 \quad \text{on } \partial\Omega_t^-, \\ (iii) & v(\psi) \cdot n \ge v_r \cdot n \quad \text{on } \partial\Omega_t^+. \end{array} \right. \qquad (7) $$

Note that working with the space $V_{\partial\Omega_t^+}$ implies $\psi = 0$ on the saturated area $\partial\Omega_t^+$. The well posedness of the stationary flow problem (i) requires that $\partial\Omega_t^+ \neq \emptyset$, i.e., that the water table reaches the top boundary. In addition, one can prove that a solution to (i)–(iii) exists whenever $i \le k_s$, i.e., whenever the rainfall intensity does not exceed the soil infiltration capacity.

An approximate solution $\{\partial\Omega_t^+, \psi\}$ of Eq. (7) is sought using a fixed-point iterative scheme as follows:

(1) choose an initial $\partial\Omega_t^+$;

(2) solve problem (i);

(3) check whether (ii) and (iii) are satisfied;

(4) if (ii) is satisfied and (iii) is not, move $\partial\Omega_t^+$ one mesh cell (or more) downhill; go to step 2;

(5) if (iii) is satisfied and (ii) is not, move $\partial\Omega_t^+$ one mesh cell (or more) uphill; go to step 2;

(6) if both (ii) and (iii) are satisfied, then the current pair $\{\partial\Omega_t^+, \psi\}$ is the desired solution.

In this algorithm, the set of points lying on the soil surface where saturation conditions are met is determined iteratively and the boundary between $\partial\Omega_t^+$
and $\partial\Omega_t^-$ always coincides with a grid node. Nothing prevents the use of the OTM and the previous algorithm in situations with multiple seepage faces, i.e when there are multiple intersections of the water table with the ground surface.

Several numerical methods can be employed to solve problem (i) in step 2. In this work, the Richards equation is discretized using linear and continuous finite elements on an unstructured mesh covering the whole computational domain. The discrete equations are nonlinear owing to the soil hydraulic functions. These equations are solved approximately using Newton's method. The convergence criterion is based on a weighted $L^2$-norm of the solution update. In the results presented in the following sections, a tolerance of $10^{-7}$ is used for Newton's method and a tolerance of $10^{-8}$ is used for the linear solver within each step of Newton's method. Since the Newton solver is called several times within the fixed-point iterative scheme, it is efficient to use the previous pair $\{\partial\Omega_t^+, \psi\}$ as a starting estimate once the saturated region $\partial\Omega_t^+$ has been updated. Once the fixed-point iterative scheme has converged, it is possible to refine the mesh and go back to step 2 to improve the quality of the solution.

Owing to the maximum principle, both conditions (ii) and (iii) in Eq. (7) cannot be violated simultaneously. However, in numerical approximations, this can 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 on whether the converged position $\partial\Omega_t^+$ has been approached from below or above. The two resulting values yield lower and upper bounds for the water table position that typically differ from one or two mesh cells at the most.
### 2.3.2 Unstationary OTM

The unstationary OTM is solved using an implicit Euler scheme with an adaptive time step. For a time step $k \ge 0$, given $(\partial\Omega_t^+)^k$ and $\psi^k$, we seek $(\partial\Omega_t^+)^{k+1} \subset \partial\Omega_t$ and $\psi^{k+1} \in V_{(\partial\Omega_t^+)^{k+1}}$ such that

$$ \left\{ \begin{array}{ll} (i) & a'_{(\partial\Omega_t^+)^{k+1}}(\psi^{k+1}, \phi) = 0 \quad \forall \phi \in V_{(\partial\Omega_t^+)^{k+1}}, \\ (ii) & \psi^{k+1} \le 0 \quad \text{on } (\partial\Omega_t^-)^{k+1}, \\ (iii) & v(\psi^{k+1}) \cdot n \ge v_r \cdot n \quad \text{on } (\partial\Omega_t^+)^{k+1}, \end{array} \right. \qquad (8) $$

where

$$ a'_{(\partial\Omega_t^+)^{k+1}}(\psi^{k+1}, \phi) = \frac{1}{\delta t} \int_{\Omega} (\theta(\psi^{k+1}) - \theta(\psi^k)) \phi + a_{(\partial\Omega_t^+)^{k+1}}(\psi^{k+1}, \phi). \quad (9) $$

This problem is solved using the same iterative algorithm as that employed for the stationary OTM. In step 1, the initial choice is $(\partial\Omega_t^+)^{k+1} = (\partial\Omega_t^+)^k$. Note that in the unstationary case, problem (i) is well-posed even if the water table has not reached the top boundary.

The finite element code with Newton's method to approximate the Richards equation has been tested against analytical solutions for steady and unsteady infiltration columns. The test cases are described in Ginzburg et al. (2004).

## 3 Virtual Hillslope Simulations

In this section, the OTM is used to perform virtual hillslope simulations at the metric and hectometric scales. The metric-scale problem is used to investigate
the influence of rainfall intensity and soil hydraulic parameters within the VGM model on the dynamic response of the water table. The metric scale for these simulations allows for a very fine representation of the soil. The hectometric-scale problem focuses on the effects of geometrical properties and initial conditions since at large scales the choice of the initial condition may engender difficulties, such as downward infiltration front, depending on the soil texture.

## 3.1 Preliminary definitions

### 3.1.1 The saturated ground fraction

Let $L$ be the length of the slope and let $L_s$ be the portion of the hillslope that is saturated. Henceforth, the ratio $L_s/L$ is termed the saturated ground fraction.

### 3.1.2 Infiltration, exfiltration, and runoff fluxes

Let $Q_{\text{rain}} = iL(e_z \cdot n)$ be the rainfall rate. We define $Q_{\text{in}}$ to be the infiltration flux and $Q_{\text{not in}}$ to be the “direct runoff” flux, i.e., the water that never infiltrates. Hence,

$$Q_{\text{rain}} = Q_{\text{in}} + Q_{\text{not in}}. \qquad (10)$$

The exfiltration flux $Q_{\text{exf}}$ corresponds to the top surface exfiltration $Q_{\text{exf}|\partial\Omega_t}$ plus the exfiltration into the stream $Q_{\text{exf}|\partial\Omega_r}$ (subsurface flow through the right surface) if any. Hence,

$$Q_{\text{exf}} = Q_{\text{exf}|\partial\Omega_t} + Q_{\text{exf}|\partial\Omega_r}. \qquad (11)$$
The quantity $Q_{runoff}$ is defined to be the exfiltration flux plus the direct runon flux. Hence,

$$Q_{runoff} = Q_{exf} + Q_{not in} = Q_{exf} + Q_{rain} - Q_{in}. \qquad (12)$$

Note that $Q_{runoff}$ is not the instantaneous water flux at the toe of the slope and into the stream but the instantaneous water flux into the overland flow and into the stream.

To evaluate the infiltration and exfiltration fluxes associated with the discrete solution to the Richards equation 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 when continuous finite elements are used. For simplicity, the method is described when no exfiltration occurs through the right boundary. The top boundary $\partial\Omega_t$ is divided into two regions, an exfiltration region, $\partial\Omega_{ex}$, where $v(\psi) \cdot n > 0$, and an infiltration region, $\partial\Omega_{in}$, where $v(\psi) \cdot n \le 0$. Let $\varphi_1$ and $\varphi_2$ be the continuous piecewise linear functions defined on $\partial\Omega_t$ such that $\varphi_1 = 1$ (respectively $\varphi_2 = 1$) at all the mesh nodes on $\partial\Omega_{ex}$ (respectively on $\partial\Omega_{in}$). Then, there exist two continuous piecewise linear functions defined on the whole computational domain, say $\phi_1$ and $\phi_2$, such that: $\phi_1|_{\partial\Omega_t} = \varphi_1$, $\phi_2|_{\partial\Omega_t} = \varphi_2$, and $\phi_1 + \phi_2 = 1$ on $\Omega$. For a time step $k \ge 1$, the infiltration flux $Q_{in}(t^k)$ and the exfiltration flux $Q_{exf}(t^k)$ are defined as follows:

$$Q_{exf}(t^k) := \int_{\partial\Omega_t} (v(\psi^k) \cdot n) \varphi_1 = \int_{\Omega} \nabla\phi_1 \cdot \nabla\psi^k + \int_{\Omega} \phi_1 \frac{\theta^k - \theta^{k-1}}{\delta t}, \qquad (13)$$

$$Q_{in}(t^k) := \int_{\partial\Omega_t} (v(\psi^k) \cdot n) \varphi_2 = \int_{\Omega} \nabla\phi_2 \cdot \nabla\psi^k + \int_{\Omega} \phi_2 \frac{\theta^k - \theta^{k-1}}{\delta t},$$
so that the following mass balance equation holds:

$$Q_{\text{exf}}(t^k) + Q_{\text{in}}(t^k) = \int_{\Omega} \frac{\theta^k - \theta^{k-1}}{\delta t}. \qquad (14)$$

At steady-state, $Q_{\text{in}} = Q_{\text{exf}}$ and, hence, $Q_{\text{runoff}} = Q_{\text{rain}}$. Finally, the time to reach equilibrium, $T_e$, is defined as the lowest time for which $|Q_{\text{in}} - Q_{\text{exf}}| \le 5 \times 10^{-3} Q_{\text{in}}$.

### 3.1.3 The VGM and the modified VGM model

Let $n > 1$ and $\alpha$ be the van Genuchten soil parameters, and set $m = 1 - \frac{1}{n}$. Using the additional parameter $\psi_s < 0$, the modified VGM model discussed in Vogel et al. (2001) is

$$\tilde{\theta}(\psi) = \begin{cases} \beta (1 + (-\alpha\psi)^n)^{-m}, & \psi < \psi_s, \\ 1, & \psi \ge \psi_s, \end{cases} \qquad (15)$$

where $\tilde{\theta}(\psi)$ is the effective saturation defined by Eq. (15) and $\beta = (1 + (-\alpha\psi_s)^n)^m$ so that $\tilde{\theta}(\psi_s) = 1$. We set $k(\psi) = k_s k_r(\psi)$ where the relative hydraulic conductivity is specified as a function of $\tilde{\theta}(\psi)$ as follows:

$$k_r(\tilde{\theta}) = \begin{cases} \tilde{\theta}^{1/2} \frac{\left[1 - \left(1 - (\beta^{-1}\tilde{\theta})^{\frac{1}{m}}\right)^m\right]^2}{\left[1 - \left(1 - \beta^{-\frac{1}{m}}\right)^m\right]^2}, & \psi < \psi_s, \\ 1, & \psi \ge \psi_s. \end{cases} \qquad (16)$$

The parameter $\psi_s$ is interpreted in Vogel et al. (2001) as a minimum capillary height. Eqs. (15) and (16) reduce to the original VGM model when $\psi_s = 0$.
Five soils, whose hydraulic parameters are collected in Table 1, are selected for our virtual hillslope simulations: Silty Clay-Loam (SCL), Yolo Clay-Loam (YLC), Sand 2 from Kao *et al.* (2001), Sand 1 from Vogel *et al.* (2001), and Sand OW from Ogden and Watts (2000).

## 3.2 Metric-scale problem

The two-dimensional domain selected to perform the study is sketched in Figure 2; it corresponds to the configuration considered by Abdul and Gillham (1984). The domain dimensions are 1.4 m in length and 1 m to 0.8 m in height. For the left boundary we assume an impermeable layer (BC1). The initial condition is a horizontal water table located at 0.7 m and the hydrostatic pressure profile. Simulations are run using a constant rainfall intensity for a duration longer than the time necessary to reach equilibrium. After assessing grid refinement issues in the next section, we investigate the influence of rainfall intensity on model predictions for YLC. Then we study the influence of the soil hydraulic parameters for SCL, YLC, and Sand 2. Particular attention is given to the shape of both the soil water retention and the relative hydraulic conductivity functions, especially near saturation.

### 3.2.1 Mesh density assessment

We consider YLC subjected to a constant rainfall intensity such that $i/k_s = 10\%$. Results obtained on two meshes are compared. The reference mesh, Mesh 1, is an unstructured triangulation consisting of 4083 nodes and 7903 elements: the node spacing is uniform and 1 cm on the top surface (142 nodes on the top surface) and 5 cm on the bottom surface. To test local mesh re-
finement, a new mesh, Mesh 1*, is produced from Mesh 1 by imposing a node
spacing of 0.5 cm at the equilibrium water table position that results from
the computation on Mesh 1. Mesh 1* consists of 7429 nodes and 14474 ele-
ments (240 nodes on the top surface). The agreement between the steady-state
results produced by the two meshes is very satisfactory with relative errors
lower than 1%; see Table 2. Since local mesh refinement improves the results
only marginally in the present simulations, quasi-uniform meshes are used in
the rest of this work. Furthermore, to test mesh coarsening, a series of meshes
(Mesh 2, Mesh 3, and Mesh 4) is produced by diminishing the number of nodes
at the top surface. Results, which are collected in Table 2, indicate that the
ratio Ls/L is only mildly affected by mesh coarsening while the equilibrium
time is somewhat more sensitive to it.

3.2.2 Influence of rainfall intensity

We consider YLC subjected to four different rainfall intensities such that
$i/k_s = 2, 5, 10,$ and $15\%$. Figure 3a plots the time evolution of the saturated
ground fraction for the four rainfall intensities. As expected, the saturated
ground fraction increases with the rainfall intensity. The nonlinear temporal
response of the soil under different rainfall intensities is related to the soil
hydraulic functions and, in particular, to the fact that for YLC, the relative
hydraulic conductivity increases rapidly with pressure head near saturation.
Figure 3b plots the time evolution of the relative exfiltration flux, *i.e.*, the
exfiltration flux normalized by the rainfall rate. We observe that it is not *a
priori* possible to rank at all times the exfiltration fluxes monotonically in
terms of the rainfall intensity. The same phenomenon is observed with SCL
and Sand 2.
### 3.2.3 Influence of the soil hydraulic parameters

The time evolution of the water table location is studied for SCL, YLC, and Sand 2. The corresponding relative hydraulic conductivity and soil water capacity, $C(\psi) = \frac{d\theta}{d\psi}$, are plotted in Figure 4 as a function of pressure head for the original VGM model. Simulations are run with a constant rainfall intensity such that $i/k_s = 10\%$ for each soil. Table 3 collects the predicted equilibrium times for the three soils. Sand 2 exhibits the fastest response and SCL the slowest. The parameter that exerts the major influence on the equilibrium time appears to be the hydraulic conductivity at saturation. Since the ratio $i/k_s$ is kept fixed, the rainfall intensity increases as the hydraulic conductivity at saturation does, and less time is then needed to fill the domain. Other parameters (such as $(\theta_s - \theta_r)$, $n$, and $\alpha$) also influence the equilibrium time, but to a smaller extent. Further studies should include a sensitivity analysis to assess quantitatively the relative importance of each parameter.

Figure 5a presents the time evolution of the saturated ground fraction. The shape of the three curves is a consequence of the respective shapes of the hydraulic functions for the three soils. At steady-state, Figure 5b shows strongly different distributions of the effective saturation along the unsaturated ground surface, although the saturated ground fraction is essentially identical for the three soils as it is controlled by the rainfall intensity for mass conservation reasons. The value of the parameter $n$ appears to condition the curvature of the plots. For instance, a sole modification of this parameter (by setting it to 5 for YLC) yields a similar response to that observed for Sand 2. For the three soils, the time evolution of the saturated ground fraction is plotted in Figure 5a. As sketched in Figure 1b (not in scale), only a part of the saturated ground (from C to D) contributes to the exfiltration, the saturated zone from B to C
still allowing for some infiltration. In this case, at steady-state, the exfiltration part (CD) is independent of the soil texture and represents about 70% of the saturated length Ls. The figure showing the time evolution of the fraction of the exfiltration part (CD/L) for the three soils is similar to Figure 5a and is therefore not presented. Figure 5c presents the time evolution of the relative exfiltration and runoff fluxes for the three soils. For each soil, the dynamics of the exfiltration are similar to that of the water table position. Furthermore, at equilibrium, the exfiltration represents 45.9% of runoff and direct runon represents 54.1%. Finally, Figure 5d presents the time evolution of the ratio of exfiltration to runoff fluxes for the three soils. The value of this ratio is roughly identical for the three soils: as soon as runoff begins, the ratio reaches a value between 40% and 50% and does not vary much until steady-state is reached. The differences between the curves result from 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; this value then slowly decreases to the steady-state value. This differs from the results for Sand 2 where the exfiltration ratio slowly increases to a maximum value that is reached just before equilibrium. This type of behavior is linked to the shape of the soil hydraulic functions.

To conclude this section, we consider the modified VGM model. In particular, we address the influence of the minimum capillary height on model results for SCL, YLC, and Sand 2. The effect of the modified VGM model on the relative hydraulic conductivity and water capacity functions is shown in Figure 6 for each soil. The same value, ψs = -2 cm, has been chosen for the three soils. Simulations are run with a constant rainfall intensity such that i/ks = 10%. The time evolution of the water table is shown in 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. As discussed by Vogel *et al.* (2001), the modification of the VGM model affects more significantly those soils for which the parameter *n* takes values close to 1. For Sand 2, the hydraulic functions obtained with the modified VGM model are very similar to those obtained with the original VGM model. This fact is reflected in the results presented in Figures 7 and 8.

### 3.3 Hectometric-scale problem

The geometry studied by Ogden and Watts (2000) is sketched in Figure 9. A no-flow boundary condition is imposed at the bottom and left surfaces, simulating an impermeable layer. The initial condition corresponds to a horizontal water table located at the toe of the slope and a hydrostatic pressure profile. The soil consists of Sand OW to facilitate comparisons with the results reported in Ogden and Watts (2000). Simulations are run with a constant rainfall intensity of *i* = 30 mm/h (*i*/*k*<sub>s</sub> = 0.6%).

The effects of hillslope geometry on model predictions is investigated first. The land surface slope *So* is always set to 10%, while the depth to impermeable layer, *D*, and the slope length, *L*, are varied. Three geometries are studied: Geometry 1 (*L* = 50 m, *D* = 1 m), Geometry 2 (*L* = 50 m, *D* = 2 m), and Geometry 3 (*L* = 30 m, *D* = 1 m). A constant head boundary condition representing a stream is imposed at the right surface (BC2). Figure 10 presents the dynamic response of the hillslope for the three geometries. The mesh consists of 11996 elements for Geometry 1, of 7196 elements for Geometry 2,
and of 14877 elements for Geometry 3. At equilibrium, Ogden and Watts have derived a simple expression to evaluate the saturated ground fraction based on the assumption that the hydraulic gradient is equal to the land surface slope So wherever the water table intersects the ground surface; the resulting formula is (Ogden and Watts (2000))

$$ \frac{L_s}{L} = 1 - \frac{Dk_sSo}{iL}. \qquad (17) $$

Our results, see Table 4, agree well with the above formula, indicating that at equilibrium, the velocity field is mainly tangential to the top surface. Since BC2 is imposed, the exfiltration flux essentially occurs through the right surface and is proportional to D; for instance, $Q_{exf}$ is found to be twice as large when D is doubled. The time to reach equilibrium naturally increases with increasing soil depth and slope length.

Figure 11 compares the two boundary conditions (BC1 and BC2) for Geometry 1. At this space scale, the impact of the boundary condition is negligible. The differences between the two solutions are localized at the toe of the slope (in the area near the stream) and locally change the nature of the exfiltration process: through the top surface for BC1 and directly into the stream for BC2.

As can be seen from Figure 12a, the relative hydraulic conductivity of Sand OW decreases very slowly compared to that of YLC or Sand 1. With this type of soil and the initial condition considered above, a numerical solution can be obtained on a reasonable mesh size (mesh spacing of around 0.5 m). The choice of another soil, such as YLC for example, requires a much finer mesh if the same initial condition is used. Indeed, the low hydraulic conductivity of the soil on the left part of the geometry induces a downward infiltration
front. To compute a steady-state solution for YLC without excessive computational effort, the initial condition illustrated in Figure 12b is considered. Figures 12c and 12d present, respectively, the time evolution of the saturated ground fraction and the relative infiltration and exfiltration fluxes for YLC with the rainfall intensity $i/k_s = 10\%$ and Geometry 1. The relative exfiltration flux is in agreement with the formula derived in Ogden and Watts (2000), namely $\frac{Q_{\text{exf}}}{Q_{\text{rain}}} = \frac{k_s SoD}{i(e_z \cdot n)L} \approx 0.02$.

## 4 Coupling with Overland Flow

The OTM examined so far is suitable to describe the coupling between the subsurface water table and surface runoff provided the height of the overland flow as well as re-infiltration processes can be neglected. A more sophisticated approach described in Esclaffer (2003) and Beaugendre et al. (2004) consists of modeling the overland flow using a simplified form of the shallow water equations, the so-called diffusive wave approximation which is widely used to describe flood routing (Moussa and Bocquillon (1996); Beauquillon (1978); Daluz-Vieira (1983)). This section first describes briefly the diffusive wave approximation for the shallow water equations and how this model can be coupled to the Richards equation. Then, numerical results obtained with the OTM and the coupled model are compared at the metric and hectometric scales.
4.1 The Coupled Subsurface/Overland Flow Model (CSOFM)

Transient flow of shallow water can be described by the following equations
(Esteves et al. (2000); Zhang and Cundy (1989))

$$
\left\{
\begin{array}{l}
\displaystyle \frac{\partial}{\partial t} y + \frac{\partial}{\partial x} (yV) = w, \\[1em]
\displaystyle \frac{\partial}{\partial t} V + V \frac{\partial}{\partial x} V + g \frac{\partial}{\partial x} y + g(S_f - S) = 0,
\end{array}
\right.
\qquad (18)
$$

were $y$ is the water depth (L), $V$ the $x$-component of the flow velocity (LT$^{-1}$),
$w$ the mass source term, $g$ the gravity (LT$^{-2}$), $S$ the river bed slope, and $S_f$
the energy line slope ($S$ and $S_f$ are dimensionless).

We assume that the first and second terms in the momentum equation can
be neglected in comparison with the third and fourth terms. This yields the
diffusive wave approximation in the form

$$
\left\{
\begin{array}{l}
(i) \quad \displaystyle \frac{\partial}{\partial t} y + \frac{\partial}{\partial x} (yV) = w, \\[1em]
(ii) \quad \displaystyle \frac{\partial}{\partial x} y + S_f - S = 0.
\end{array}
\right.
\qquad (19)
$$

The Manning-Strickler uniform flow formula (Viollet *et al.* (1998)) is chosen to
describe the energy line slope $S_f$ (Panday and Huyakorn (2004); VanderKwaak
(1999); Liu *et al.* (2003))

$$
V = K_S R^{\frac{2}{3}} \operatorname{sgn}(S_f) |S_f|^{\frac{1}{2}}, \qquad (20)
$$

where $K_S$ is the Strickler coefficient of roughness (L$^{-\frac{1}{3}}$T$^{-1}$) and $R$ the hydraulic radius (L) defined as the ratio between the cross-sectional flow area $A$ (L$^2$) and the wet perimeter $\chi$ (L). Assuming that the overland flow occurs as
a thin layer with a wide rectangular section, the relation $R \approx y$ holds. The Manning-Strickler uniform flow formula is also considered in the InHM model, VanderKwaak and Loague (2001).

Letting $q(y) = yV$, Eqs. (19)(ii) and (20) readily yield

$$q(y) = K_S y^{\frac{5}{3}} \operatorname{sgn}(S_f) |S - \frac{\partial}{\partial x} y|^{\frac{1}{2}}. \qquad (21)$$

This is conveniently rewritten as follows:

$$q(y) = K_S y^{\frac{5}{3}} \frac{S - \frac{\partial}{\partial x} y}{|S - \frac{\partial}{\partial x} y|^{\frac{1}{2}}}. \qquad (22)$$

The continuity equation then becomes

$$\frac{\partial}{\partial t} y + \frac{\partial}{\partial x} \left( K_S y^{\frac{5}{3}} \frac{S - \frac{\partial}{\partial x} y}{|S - \frac{\partial}{\partial x} y|^{\frac{1}{2}}} \right) = w. \qquad (23)$$

One advantage of the diffusive wave approximation is to reduce the shallow water equations to a single partial differential equation. However, by doing so, the differential order of the equation is increased, and it is now necessary to supply two boundary conditions for Eq. (23). The condition $y = 0$ is imposed on $\partial\Omega_t^-$ and expresses the fact that the overland flow only exists on $\partial\Omega_t^+$ where it flows downhill. At the hillslope outlet, another simplification, the so-called kinematic wave approximation, i.e., $\frac{\partial}{\partial x}y = 0$, is made. In other models like MODHMS (Panday and Huyakorn (2004)) and HydroGeoSphere, based on earlier work of Therrien and Sudicky (1996), it is possible to choose between the Chezy and the Manning-Strickler formula.

The CSOFM is formulated as follows. On $\partial\Omega_t^+$, the boundary condition for the Richards equation becomes $y = \psi$. This condition expresses the fact that the vertical pressure profile in the overland flow is hydrostatic. Moreover, the
source term $w$ in Eq. (23) is evaluated from the rainfall and exfiltration fluxes
in the form $w = v_r \cdot n - v(\psi) \cdot n$. More details on the computational method
used to approximate the CSOFM are given in the appendix.

Since the CSOFM assumes the relation $y = \psi$, this model does not account for
changes in water storage and flow area due to micro-topographic depressions
and/or obstructions as proposed in more general models such as MODHMS,
HydroGeoSphere, and InHM. Since the coupling between the overland and
subsurface flows relies on the continuity of the hydraulic head, only one pres-
sure value is defined at each node of the ground surface. On the contrary,
MODHMS, HydroGeoSphere, and InHM allow for a discontinuity of the pres-
sure field at the ground surface nodes. Then, an evaluation of the exchanged
flux from the overland flow to the subsurface flow is needed (Gerke and van
Genuchten (1993)).

## 4.2 Numerical results

This section compares the numerical results obtained with the OTM and the
CSOFM at the metric and hectometric scales.

### 4.2.1 Metric-scale problem

The dynamics of the water table are studied for Sand 1 using both BC1 and
BC2. Here and in subsequent results, the Strickler coefficient of roughness $K_S$
is set to $10 \text{ m}^{-\frac{1}{3}}\text{s}^{-1}$. The geometry is that described in Section 3.2. The initial
condition is a horizontal water table located at 0.8 m (the toe of the slope) and
a hydrostatic pressure profile. Figure 13 and 14 present the time evolution of
the saturated ground fraction and of the relative infiltration and exfiltration
fluxes, respectively. The OTM and the CSOFM yield similar results for both
boundary conditions. Both models predict that the final water table position
is lower with BC2 than with BC1. This is related to the fact that BC1 and
BC2 yield different flow velocity fields near the right boundary: with BC2,
the velocity vectors on the saturated zone are nearly tangent to the surface
so that most of the exfiltration occurs through the right surface. Figure 14c
presents the time evolution of the relative runoff flux for the two models and
the two boundary conditions. At this space scale, the boundary condition
clearly affects the dynamics of runoff. Again, both models yield similar results.

4.2.2 *Hectometric-scale problem*

We consider the geometry of Ogden and Watts described in Section 3.3. The initial condition is a horizontal water table located at the toe of the slope and a hydrostatic pressure profile. Simulations are run with a constant rainfall intensity such that *i* = 30 mm/h, and BC2 is imposed at the right surface. Figure 15 compares the results obtained with the two models on Geometry 1. Under the present conditions, both models yield similar results with the exception of a slight difference in the steady-state values predicted for the exfiltration flux. This is related to the fact that because the pressure takes slightly larger values along the top boundary when predicted with the CSOFM, its normal derivative is larger, and, hence, the exfiltration flux is also larger.

**5 Conclusions**

The Richards equation combined with different models for hydraulic functions
(BCM model (Brooks and Corey (1964)); original (van Genuchten (1980)) and
modified (Vogel *et al.* (2001)) VGM model) has been discretized by finite elements in space. The nonlinearities induced by the soil hydraulic functions can be efficiently handled by Newton's method. From a numerical viewpoint, the soil labelled SCL is the most difficult to deal with because of its strong non-linear behavior especially near saturation. A second nonlinearity stems from the fact that the movement of the water table induces a changing character of the surface boundary condition for the Richards equation. This yields a unsteady obstacle-type problem for which a robust solution procedure based on a fixed-point iteration has been proposed.

The obstacle-type model has been used to investigate the hydraulic behavior of virtual hillslopes under constant rainfall conditions. Results on the small-scale geometry have shown that the rainfall intensity has a significant influence on the spatial extension of the saturated area of the hillslope and on the exfiltration flux. For a fixed value of the ratio $i/k_s$, the structure of the soil, characterized by the van Genuchten hydraulic parameters, strongly affects the dynamic response of the system, whereas the steady-state values remain fairly identical. For the three soils addressed in the metric-scale problem, exfiltration occurs as soon as runoff does; this means that at any time, water exiting the system will have roughly the same origin, either from direct runoff or through subsurface flow. In addition, a driving parameter for the time evolution of the saturated ground fraction and of the relative exfiltration flux appears to be the exponent $n$ in the VGM model. Further investigations should be performed to determine whether its influence expresses itself through the relative hydraulic conductivity or the soil water capacity or both functions. Furthermore, results similar to those of Vogel *et al.* (2001) for infiltration in very dry soils are observed with the modified VGM model, especially for fine-textured soils where
the introduction of a minimum capillary height has a significant influence on runoff and exfiltration genesis. Based on the above findings, a sensitivity analysis can be conducted as the next step to achieve a better understanding of the parameters conditioning surface runoff in hillslopes.

On the hectometric-scale problem, the numerical results are in very good agreement with the simple closed formulas derived by Ogden and Watts. Moreover, at this space scale, boundary conditions on the right surface have a marginal impact on the solution. This medium-scale problem also underlines the importance of initial conditions in virtual hillslope simulations.

Finally, we have compared the obstacle-type model with a coupled subsurface/overland flow model in which the dynamics of the overland flow are governed by the shallow water equations in the diffusive wave approximation. Both models yield similar results for the present test cases, both at the metric and hectometric scales, because the mechanisms leading to saturation are mainly controlled by subsurface flows and, more specifically, by the hydraulic properties of the soil. In other words, neglecting the feedback of surface runoff on the water table dynamics leads altogether to reliable results for the present test cases. Subsequent work should also consider models including re-infiltration processes under different rainfall conditions to further investigate the impact of overland flow on the genesis of surface runoff. Some preliminary test cases involving multiple seepage faces have been run. More realistic situations involving heterogeneous soil conductivity are being investigated to study the performance of the OTM code when dealing with perched water tables.
Acknowledgements

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 the Estime project at INRIA for stimulating discussions. We
also thank the referees for their careful reading of the manuscript and their
suggestions.
Appendix : CSOFM

This appendix briefly describes the computational method used to approximate the CSOFM. After time-discretization with the implicit Euler scheme, the subsurface flow problem consists of seeking $\psi^{k+1} \in H^1(\Omega)$ such that for all $\phi \in H^1(\Omega)$,

$$ \frac{1}{\delta t} \int_{\Omega} (\theta(\psi^{k+1}) - \theta(\psi^k)) \phi + a_{(\partial \Omega_t^+)^{k+1}}(\psi^{k+1}, \phi) = R(\psi^k, \psi^{k+1}, \phi). \quad (24) $$

Eq. (23) and the boundary condition $y = \psi$ on $\partial\Omega_t^+$ are used to approximate the right-hand side as follows:

$$ R(\psi^k, \psi^{k+1}, \phi) = - \int_{\partial\Omega_t^+} \left( \frac{\psi^{k+1} - \psi^k}{\delta t} \phi + \frac{\partial}{\partial x} q(\psi^{k+1}) \phi \right). \quad (25) $$

Mass lumping is used for the unsteady term, i.e., for the vertex $N_i$ located on $\partial\Omega_t^+$,

$$ \int_{\partial\Omega_t^+} \left( \frac{\psi^{k+1} - \psi^k}{\delta t} \right) \phi \approx \frac{\psi^{k+1}(N_i) - \psi^k(N_i)}{\delta t} L_i, \quad (26) $$

where $L_i$ is the distance separating the midpoints of the boundary cells on the left and on the right of the vertex $N_i$. An upwind scheme is used to approximate the spatial derivative in the right-hand side of Eq. (25). Assuming that the nodes on $\partial\Omega_t^+$ are numbered downwards and defining the numerical flux

$$ Q_{i-1,i}^{k+1} = K_s(\psi^{k+1}(N_{i-1}))^{\frac{5}{3}} \frac{S - \frac{\psi^{k+1}(N_i) - \psi^{k+1}(N_{i-1})}{L_{i-1/2}}}{\left| S - \frac{\psi^{k+1}(N_i) - \psi^{k+1}(N_{i-1})}{L_{i-1/2}} \right|^{\frac{1}{2}}}, \quad (27) $$
where $L_{i-\frac{1}{2}}$ denotes the distance between $N_{i-1}$ and $N_i$, we set

$$ \int_{\partial \Omega_t^+} \frac{\partial}{\partial x} q(y) \phi \approx Q_{i,i+1}^{k+1} - Q_{i-1,i}^{k+1}. \qquad (28) $$

Eq. (28) is used for all nodes in the interior of $\partial \Omega_t^+$. At the outlet node, say $N_{i_0}$, we set

$$ Q_{i_0,i_0+1}^{k+1} = K_S(\psi^{k+1}(N_{i_0}))^{\frac{5}{3}} S^{\frac{1}{2}}. $$

With this choice, the boundary condition $\frac{\partial}{\partial x}y = 0$ is weakly enforced at the hillslope outlet.

Considering BC2 amounts to enforce a constant total hydraulic head along the right boundary. At each time $t$, this value can be evaluated using the node located at the toe of the slope.

## References

Abdul A. S. and Gillham R. W., Laboratory studies of the effects of the capillary fringe on streamflow generation, *Water Resour. Res.*, 20(6):691–708, 1984.

Arfken G., Hypergeometric Functions, Mathematical Methods for Physicists, 3rd ed. Orlando, Academic Press, 13.5, 748–752, 1985.

Barros A. P., Knapton D., Wang M. C. and Kuo C. Y., Runoff in shallow soils under laboratory conditions, *J. Hydraul. Eng.*, 4(1):28–37, 1999.

Barry D. A., Parlange J.-P., Sander G. C. and Sivaplan M., A class of exact solutions for Richards’ equation, *J. of Hydrology*, 142:29–46, 1993.

Bastian P., Higher order discontinuous Galerkin methods for flow and trans-
port in porous media, Challenges in scientific computing – CISC 2002 (E. Baensch, ed.), Springer. Oxf. Texts Appl. Eng. Math. 35:1–22, 2003.

Bause M. and Knabner P., Computation of variably saturated subsurface flow by adaptive mixed hybrid finite element methods, *Advances in Water Resources*, 27:565–581, 2004.

Beaugendre H., Ern A., Esclaffer T. and Gaume E., Numerical investigation of surface runoff in hillslopes with variably saturated flows, ECCOMAS Proceedings, 2004.

Beauquillon C., Propagation des écoulements transitoires intermittents dans les systèmes hydrauliques à surfaces libres, *Rapport technique, Université Montpellier II*, Laboratoire d'Hydrologie Mathématique, Montpellier, 1978.

Betson R. P., What is watershed runoff? *J. Geophys. Res.*, 69(8):1541–1551, 1964.

Beven K., Hillslope hydrographs by the finite element method, *Earth Surface Processes*, 2:13-28, 1977.

Brooks R. H. and Corey A. T., Hydraulic properties of porous media, *Hydrol. Pap.*, Colo. State Univ., Fort Collins, 3:27, 1964.

Cappus P., Etude des lois de l'écoulement. Application au calcul et à la prévision de débits. Bassin expérimental de l'Alrance. La Houille Blanche A:493-520, 1960.

Celia M. A., Boulotas E. T. and Zarba R. L., A general mass-conservative numerical solution for the unsaturated flow equation, *Water Resour. Res.*, 26(7):1483–1496, 1990.

Cloke H. L., Renaud J.-P., Claxton A. J., McDonnell J. J. Anderson M. G., Blake J. R. and Bates P. D., The effect of model configuration on modelled hillslope-riparian interactions, *J. of Hydrology*, 279:167–181, 2003.

Cloke H. L., Anderson M. G., McDonnell J. J., Renaud J.-P., Using numerical
modelling to evaluate the capillary fringe groundwater ridging hypothesis of
streamflow generation, *J. of Hydrology*, 316:141–162, 2006.

Daluz-Vieira J.-H., Conditions governing the use of approximations for the Saint-Venant equations for shallow surface water flow, *J. of Hydrology*, 60:43–58, 1983.

Discacciati M. , Miglio E. and Quarteroni A., Mathematical and numerical models for coupling surface and groundwater flows, *Appl. Num. Math.*, 43:57–74, 2002.

Dunne T. and Black R. D., Partial area contributions to storm runoff in a small New England watershed, *Water Resour. Res.*, 6:1296–1311, 1970.

Esclaffer T., Etude théorique de la formation des débits de crues à l'échelle des versants, *Rapport de stage de DEA*, Université de Paris 12, ENGREF, ENPC-CEREVE, 2003.

Esteves M., Faucher X., Galle S. and Vauclin M., Overland flow and infiltration modelling for small plots during unsteady rain: numerical results versus observed values, *J. of Hydrology*, 228:265–282, 2000.

Fagherazzi S. , Furbish D. J., Rasetarinera P. and Hussaini M. Y., Appli-
cation of the discontinuous spectral Galerkin method to groundwater flow,
*Advances in Water Resources*, 27:129–140, 2004.

Gerke H. H. and van Genuchten M. T., A dual-porosity model for simulating the preferential flow movement of water and solutes in structured porous media, *Water Resources Research*, 29(2):305–319, 1993.

Ginzburg I., Carlier J.-P. and Kao C., Lattice Boltzmann approach to
Richards' equation, *CMWR Proceedings*, 2004.

Horton R. E., The role of infiltration in the hydrologic cycle, *Eos Trans. AGU*, 14:446–460, 1933.

Iorgulescu I., Analyse du comportement hydrologique par une approche
intégrée à l'échelle du versant, application au bassin de la Haute Mentue,
*PhD Thesis*, Ecole Polytechnique Fédérale de Lausanne, 1997.

Kao C., Bouarfa S. and Zimmer D., Steady state analysis of unsaturated flow
above a shallow water-table aquifer drained by ditches, *J. of Hydrology*,
250:122-133, 2001.

Knabner P. and Schneid E., Adaptive hybrid mixed finite element discretiza-
tion of instationary variably saturated flow in porous media, *High Per-
formance Scientific and Engineering Computing* (M. Breuer et al., eds.),
Springer Verlag Berlin, 37–44, 2002.

Liu Y. B., Gebremeskel S., De Smedt F., Hoffmann L. and Pfister L., A
diffusive transport approach for flow routing in GIS-based flood modeling,
*J. of Hydrology*, 283:91–106, 2003.

Loye-Pilot M-D. and Jusserand C., Décomposition chimique et isotopique d'un
hydrogramme de crue d'un torrent méditerranéen montagnard. Réflexions
méthodologiques, *Revue des Sciences de l'Eau*, 3:211-231, 1990.

Manzini G. and Ferraris S., Mass-conservative finite volume methods on 2-D unstructured grids for the Richards' equation, *Advances in Water Resources*, 27: 1199-1215, 2004.

Marc V. and Didon-Lescot J-F. and Michael C., Investigation of the hydrolog-
ical processes using chemical and isotopic tracers in a small Mediterranean
forested catchment during autumn recharge, *J. of Hydrology*, 247:215-229,
2001.

McGlynn B. L. and McDonnell J. J., Quantifying the relative contributions
of riparian and hillslope zones to catchment runoff, *Water Resour. Res.*,
39(11):1310, 2003.

Miglio E., Quarteroni A. and Saleri F., Coupling of free surface and ground-
water flows, *Comput. Fluids*, 32:73–83, 2003.
Moussa R. and Bocquillon C., Criteria for the choice of flood routing methods in natural channels, *J. of Hydrology*, 186:1–30, 1996.

Mualem Y., A new model for predicting the hydraulic conductivity of unsaturated porous media, *Water Resour. Res.*, 12:513–522, 1976.

Neuman S. P., Saturated-unsaturated seepage by finite elements, *J. of the hydraulics div.*, 99, (HY2), 2233–2250, 1973.

Ogden F. L. and Watts B. A., Saturated area formation on nonconvergent hillslope topography with shallow soils: A numerical investigation, *Water Resour. Res.*, 36(7):1795–1804, 2000.

Panday S. and Huyakorn P. S., A fully coupled physically-based spatially-distributed model for evaluating surface/subsurface flow, *Advances in Water Resources*, 27(4):361–382, 2004.

Richards L. A., Capillary conduction of liquids through porous mediums, *Physics*, 1:318–333, 1931.

Rulon J. J., Rodway R. and Freeze R. A., The development of multiple seepage faces on layered slopes, *Water Resour. Res.*, 21(11):1625–1636, 1985.

Therrien R. and Sudicky E. A., Three-dimensional analysis of variably-saturated flow and solute transport in discretely-fractured porous media, *J. of Contaminant Hydrology*, 23:1–44, 1996.

VanderKwaak J. E., Numerical simulation of flow and chemical transport in surface-subsurface hydrologic systems, *Ph.D. thesis*, Department of earth sciences, University of Waterloo, Ontario Canada, 1999.

VanderKwaak J. E. and Loague K., Hydrologic-response simulations for the R-5 catchment with comprehensive physics-based model, *Water Resour. Res.*, 37(4):999–1013, 2001.

van Genuchten M. Th., A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, *Soil Sci. Soc. Am.*, 44:892–898, 1980.
Viollet P.-L., Chabard J.-P., Esposito P. and Laurence D., *Mécanique des fluides appliquée*, Presses de l'Ecole Nationale des Ponts et Chaussées, Paris, 1998.

Vogel T., van Genuchten M. Th. and Cislerova M., Effect of the shape of the soil hydraulic functions near saturation on variably-saturated flow predictions, *Advances in Water Resources*, 24(6):133-144, 2001.

Weiler M. and McDonnell J., Virtual experiments: a new approach for improving process conceptualization in hillslope hydrology, *J. of Hydrology*, 285:3-18, 2004.

Woodward C. S. and Dawson C. N., Analysis of expanded mixed finite element methods for a nonlinear parabolic equation modeling flow into variably saturated porous media, *SIAM J. Numer. Anal.*, 23:517-530, 2000.

Zhang W. and Cundy T., Modeling of two-dimensional overland flow, *Water Resour. Res.*, 25(9):2019-2035, 1989.

Zhu J. and Mohanty B. P., Analytical solutions for steady state vertical infiltration, *Water Resour. Res.*, 38:1145-1149, 2002.
<figure><img src="image_3.png" /><figcaption>Fig. 1. (a) Schematic representation of the domain $\Omega$. (b) Distribution along the top surface of typical pressure head, $\psi$ (m) (solid line), and normal velocity, $v \cdot n$ (m/h) (dashed line).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 2. Abdul and Gillham geometry.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 3. (a) Time evolution of the saturated ground fraction for four different rainfall intensities with YLC. (b) Time evolution of the relative exfiltration flux for four different rainfall intensities with YLC.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 4. (a) Relative hydraulic conductivity $k_r(\psi)$ and (b) water capacity $C(\psi)$ for SCL, YLC, and Sand 2 using the original VGM model.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 5. Original VGM model. (a) Time evolution of the saturated ground fraction for SCL, YLC, and Sand 2. (b) Distribution of the effective saturation along the top surface for the three soils. (c) Time evolution of the relative exfiltration and runoff fluxes for the three soils. (d) Time evolution of the ratio between exfiltration and runoff fluxes for the three soils.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 6. (a) Relative hydraulic conductivity $k_r(\psi)$ and (b) water capacity $C(\psi)$ for the modified ($\psi_s = -2$ cm) VGM models; SCL, YLC, and Sand 2.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 7. Effect of the modified VGM model ($\psi_s = -2$ cm) on the time evolution of the saturated ground fraction for (a) SCL, (b) YLC, and (c) Sand 2.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 8. Effect of the modified VGM model ($\psi_s = -2$ cm) on the time evolution of exfiltration and runoff fluxes for (a) SCL, (b) YLC, and (c) Sand 2.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 9. Ogden and Watts geometry (not in scale).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 10. (a) Time evolution of the saturated ground fraction and (b) of the relative exfiltration flux for the three geometries.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 11. Effect of boundary condition on the time evolution of (a) the saturated ground fraction, (b) the relative exfiltration flux, and (c) the relative runoff flux.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 12. (a) Relative hydraulic conductivity for Sand 1, YLC, and Sand OW. (b) YLC initial water table level (not in scale). Time evolution of (c) the saturated ground fraction for YLC and (d) the relative infiltration and exfiltration fluxes.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 13. Comparison at the metric scale of OTM and CSOFM predictions for the time evolution of the saturated ground fraction using BC1 or BC2 for Sand 1.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 14. Sand 1: comparison of OTM and CSOFM predictions for the relative in-filtration flux and exfiltration fluxes using (a) BC1 or (b) BC2. (c) Comparison of OTM and CSOFM predictions for the relative runoff flux using BC1 or BC2.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 15. Comparison at the hectometric scale of the two models for the BC2 boundary condition. Time evolution of (a) the saturated ground fraction and (b) the relative exfiltration flux.</figcaption></figure>
<caption>Table 1</caption>

VGM soil hydraulic parameters.

<table><thead><tr><th>Texture</th><th>θ<sub>r</sub> (-)</th><th>θ<sub>s</sub> (-)</th><th>α (1/m)</th><th>n (-)</th><th>k<sub>s</sub> (m/h)</th></tr></thead><tbody><tr><td>Sand OW</td><td>0.069</td><td>0.435</td><td>0.326</td><td>3.9</td><td>5.0</td></tr><tr><td>Sand 1</td><td>0.045</td><td>0.430</td><td>14.5</td><td>2.68</td><td>0.297</td></tr><tr><td>Sand 2</td><td>0.05</td><td>0.5</td><td>3.7</td><td>5</td><td>0.1</td></tr><tr><td>YLC</td><td>0.23</td><td>0.55</td><td>3.6</td><td>1.9</td><td>0.018</td></tr><tr><td>SCL</td><td>0.1</td><td>0.41</td><td>1.9</td><td>1.31</td><td>0.0026</td></tr></tbody></table>
<caption>Table 2</caption>

Mesh influence on model predictions at steady state.

<table><thead><tr><th></th><th><i>T</i><sub>e</sub> (h)</th><th><i>Ls</i>/<i>L</i> (%)</th><th>Nodes on the top surface</th></tr></thead><tbody><tr><td>Mesh 1*</td><td>4.15</td><td>60.14</td><td>240</td></tr><tr><td>Mesh 1</td><td>4.18</td><td>60.28</td><td>142</td></tr><tr><td>Mesh 2</td><td>4.11</td><td>59.15</td><td>71</td></tr><tr><td>Mesh 3</td><td>4.33</td><td>60.00</td><td>35</td></tr><tr><td>Mesh 4</td><td>4.81</td><td>61.11</td><td>18</td></tr></tbody></table>
Table 3

Equilibrium time for three soils using the original VGM model.

<table><thead><tr><th></th><th>SCL</th><th>YLC</th><th>Sand 2</th></tr></thead><tbody><tr><td>T<sub>e</sub> (h)</td><td>10.89</td><td>4.18</td><td>0.58</td></tr></tbody></table>
Table 4

Equilibrium time and saturated ground fraction predicted numerically or using the empirical formula of Ogden and Watts.

<table><thead><tr><th></th><th>Geometry 1</th><th>Geometry 2</th><th>Geometry 3</th></tr></thead><tbody><tr><td>T<sub>e</sub> (h)</td><td>5.98</td><td>11.97</td><td>1.43</td></tr><tr><td>L<sub>s</sub>/L (%)</td><td>67.3</td><td>34.7</td><td>45.6</td></tr><tr><td>1 - <sup>Dk<sub>s</sub>So</sup>/<sub>iL</sub> (%)</td><td>66.7</td><td>33.3</td><td>44.4</td></tr></tbody></table>