subsurface flow modelling

Journal of Hydrology 366 (2009) 9–20

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

A generalized Richards equation for surface/subsurface flow modelling S. Weill, E. Mouche *, J. Patin Laboratoire des Sciences du Climat et de l’Environnement, UMR CEA-CNRS-UVSQ, Centre d’Etudes de Saclay, Orme des merisiers, 91191 Gif-sur-Yvette, France

a r t i c l e

i n f o

Article history: Received 21 February 2008 Received in revised form 1 December 2008 Accepted 9 December 2008

Keywords: Surface/subsurface interactions Integrated modelling Saturation excess runoff Infiltration excess runoff Numerical simulation

s u m m a r y An approach for modelling surface/subsurface flow and transport in a fully integrated way is presented. The diffusive wave approximation used to model runoff is formulated as a nonlinear Darcy law and coupled with the Richards equation for variably saturated flow. The water dynamics in the three physical domains, land surface, vadose zone and saturated zone, is described through a single Richards type equation with domain-dependent parameters. This approach naturally provides pressure and fluxes continuity along land surface. An advective–diffusive transport equation is incorporated in the model in order to follow and identify pre-event and event water. This multi-domain Richards equation is solved with a mixed hybrid finite element formulation. The time discretization is implicit and the nonlinear equations are solved within a sequential iterative Picard scheme. This model can describe both infiltration excess or saturation excess runoff. In order to assess this modelling approach, several classical validation, verification and application test cases are presented. For overland flow solely, the model is compared to an analytical solution and to commonly used hydrological models. The integrated model is then validated with a sandbox laboratory experiment. Finally, the dynamic response of a saturated contributing area at the toe of a 3D virtual hillslope is simulated to illustrate the applicability of the model to 3D situations. Ó 2008 Elsevier B.V. All rights reserved.

Introduction Surface processes, such as runoff or overland flow, are strongly coupled to subsurface ones and surface water in streams or rivers is in constant interaction with subsurface water from the vadose zone and shallow groundwater systems. These interactions and the mechanisms leading to surface runoff have received particular attention over the last decades (Kirkby, 1978; Beven, 2001) and are considered by several scientific communities as processes controlling water budget at the catchment scale. They are also of primary importance when studying runoff and erosion from the plot scale (e.g., Gimenez and Govers, 2001) to the field scale (e.g., Nearing et al., 2005). Thus, numerical models, which allow to evaluate velocity fields, flowpaths and relative contributions of different sources to stream, must accurately take into account surface/subsurface interactions to provide reliable predictions and to help for a better understanding of hillslope and catchment hydrological processes. Runoff generation processes are highly dependent on topography (Anderson and Burt, 1978; Ambroise, 1995), on spatial and temporal distribution of rainfall, on vegetation and on hydraulic properties of the soil (Dunne et al., 1991; Woolhiser et al., 1996). Historically, the first to provide a conceptual model for runoff generation was Horton (1993). Hortonian runoff, also called infiltra-

* Corresponding author. Present address: LSCE, Orme des Merisiers, CE de Saclay, 91191 Gif sur Yvette, France. Tel.: +33 1 69 08 22 54; fax: +33 1 69 08 77 16. E-mail addresses: [email protected] (S. Weill), [email protected] (E. Mouche), [email protected] (J. Patin). 0022-1694/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2008.12.007

tion excess runoff, may occur when rainfall intensity exceeds the soil’s infiltration capacity, which physically describes the soil’s ability to transmit rainfall water. Under these conditions, the higher part of the soil is not able to infiltrate the whole amount of rainfall. Water accumulates on the land surface and then eventually flows to the stream if it does not re-infiltrate on unsaturated zones encountered along its flow path. The time needed to reach runoff depends on the initial moisture of the land surface and on rainfall intensity (Rubin, 1969). Betson (1964) modified Horton’s concept and showed that in most of the cases only a small part of a catchment area generates Hortonian overland flow (partial area concept). The second major runoff generation concept, the saturation excess runoff, arose later on. Cappus (1960) and Dunne and Black (1970) showed that only the saturated areas of a watershed contribute to surface runoff. In this case, water table reaches the land surface level and exfiltration occurs. Two main domains are then to be considered, the infiltration and the exfiltration/runoff domains. Even if these two concepts were considered separately for many times, they may occur simultaneously if the field configuration allows it. The dynamics of water in a catchment and particularly at the surface/subsurface interface is still poorly understood. The determination of the saturated contributing surfaces and their evolution in time and space, the determination of the relative contributions of surface and subsurface water to streamflow (Sklash and Farvolden, 1979; Ribolzi et al., 2000; Kirchner, 2003; McGlynn and McDonnell, 2003) are key issues in hillslope hydrology. Some experiments have been completed to help for a better

10

S. Weill et al. / Journal of Hydrology 366 (2009) 9–20

understanding but they are expensive and complex (e.g., Abdul and Gillham, 1984; Barros et al., 1999). Recently, the concept of ‘‘virtual hillslope experiment” has appeared as an alternative approach to improve our understanding of processes in hillslope hydrology (Ogden and Watts, 2000; Weiler and McDonnell, 2004). Numerical tools to model variably saturated flow and overland flow independently are widely available. Efficient models and robust algorithms have been developed to study variably saturated– unsaturated flows (Celia et al., 1990; Paniconi and Putti, 1994; Huyakorn et al., 1996; Therrien and Sudicky, 1996) and overland flow (Gottardi and Venutelli, 1993; Giammarco et al., 1996; Singh and Bhallamudi, 1998). The numerical coupling of surface and subsurface flows in integrated hydrological models has been a widely recognized research area since Freeze and Harlan (1969) laid out the guidelines in 1969. Many attempts were made to develop models that could handle both surface and subsurface flow dynamics. Smith and Woolhiser (1971) used a 1D unsaturated flow equation coupled to a 1D, kinematic overland flow model. Freeze (1972), Freeze (1972) and Govindaraju and Kavvas (1991) coupled a 3D variably saturated groundwater equation with a 1D channel and overland flow model to examine the general response of watershed and especially the role of saturated surface areas. The Système Hydrologique Européen (SHE) original model (Abbott et al., 1986a; Abbott et al., 1986b) could describe 1D channel flow, two dimensional overland flow, 1D vertical unsaturated flow and 2D saturated groundwater flow. Esteves et al. (2000) coupled a 2D channel flow equation with a 1D Green and Ampt infiltration solution to model Hortonian runoff. Beaugendre et al. (2006) coupled 2D variably saturated flow with an obstacle-type algorithm mainly to study the influence of soil’s properties and geometry on runoff generation. In most of these works, surface/subsurface processes coupling is accomplished using matching boundary condition techniques and water routing at the land surface. When ponding occurs, boundary condition at the land surface is switched from a specified rainfall flux to a prescribed head equal to the water depth. Water that cannot infiltrate is then routed to the stream by the overland flow model. As subsurface response is both variable in time and space, this method may sometimes lead to algorithmic and computational difficulties. For the last few years, hydrological modeller have been developing integrated models that take into account these processes altogether and fully couple overland flow and subsurface flows (VanderKwaak and Loague, 2001; Panday and Huyakorn, 2004; Kollet and Maxwell, 2006). Traditionally, variably saturated flow is described using Richards equation (Richards, 1931). Overland flow is described either with a diffusive wave or a kinematic wave equation. In the works of VanderKwaak (1999) and Panday and Huyakorn (2004), coupling is done via a first order exchange flux between surface and subsurface domains (conductance concept). The coupling flux, characterized by the head gradient between the two domains and a proportionality constant, is then treated as a general source/sink term and the proportionality constant can be used as a fitting parameter. Kollet and Maxwell (2006) modified the matching boundary method presented before so that the surface water equation is used as a boundary condition to the Richards equation. This method provides pressure and flux continuities at the ground surface. This boundary condition is head-dependent and allows to remove the conductance term and the exchange flux between surface and subsurface domains. This study presents a new modelling approach which describes surface and subsurface flows in an integrated way. This paper focuses on the physics and the numerics of the new formulation. The application of the model at the catchment scale and the related issues such as, equifinality (e.g., Beven, 2001) or the impact of soil

and landcover heterogeneities (e.g., Woolhiser et al., 1996) are not considered here, as the model has been applied only at the plot scale. For overland flow modelling, we use the diffusive wave approximation and the resulting equation is considered mathematically as a Richards equation. Thus, we introduce in the computational domain a porous layer, called runoff layer, at the land surface in which surface processes are simulated. A single Richards equation with domain-dependent parameters is then built to describe all surface and subsurface processes. The whole flow dynamics is described in a single continuum, where all the water fluxes are proportional to the hydraulic gradient, from saturated zone to land surface, and pressure and velocities are continuous through the land surface. In a similar way a transport equation describing tracer migration in this continuum is proposed. Using this formulation, flow and transport in surface and subsurface are described in a unified manner. This approach was implemented in the finite element code Cast3M (website http://www.cast3M.fr). The single nonlinear flow and transport equations are solved respectively with mixed hybrid finite element and finite volume formulations (Mosé et al., 1994; Dabbenne, 1998; Michel et al., 2004). The time discretization scheme is fully implicit and nonlinearities are solved with a Picard algorithm (Paniconi and Putti, 1994; LePotier et al., 1998). We first present validation examples for overland flow solely. Then our integrated surface–subsurface flow and transport model is validated with the Abdul and Gillham sandbox laboratory experiment (Abdul and Gillham, 1984). Finally, the dynamics of a saturated area at the toe of a 3D virtual hillslope (e.g., Govindaraju and Kavvas, 1991) is simulated to illustrate the applicability of the model to 3D situations. Development of the governing equations The equations governing water dynamics in surface and subsurface domains are presented in this section. These equations have been discussed in great detail in the literature, we therefore only provide the general context and the mathematical development made in the model. The advective–dispersive equation for transport modelling is also described. More details are given in Weill (2007). Unsaturated and saturated flows Variably saturated flow in deformable porous media is modelled with the Richards equation (Richards, 1931):

rðhÞ

! oh !  r:ðKðhÞ r ðh þ zÞÞ ¼ q ot

ð1Þ

where t is time [T], h is the pressure head [L], z is the elevation [L], rðhÞ is the specific volumetric storativity [L1], KðhÞ is the hydraulic conductivity [LT1] and q represents sink and/or source of water [T1]. If we neglect hydro mechanical coupling in the unsaturated zone, which is justified for many soils, the specific volumetric storativity reduces to the soil capillary capacity CðhÞ = oh/oh in the unsaturated zone [L1], where h is the water content in the unsaturated zone. In the saturated zone, the specific volumetric storativity reduces to the specific volumetric storage S (Bear and Bachmat, 1991). The specific storativity and the hydraulic conductivity may be therefore defined in the whole pressure domain as follows:

rðhÞ ¼ KðhÞ ¼



CðhÞ for h < 0 for h P 0 S  K sat kr ðhÞ for h < 0 K sat

for h P 0

ð2Þ ð3Þ

where K sat is the saturated hydraulic conductivity of the medium [LT1] and kr ðhÞ its relative permeability. The main advantage using

S. Weill et al. / Journal of Hydrology 366 (2009) 9–20

this pressure head formulation is that both saturated and unsaturated flows can be simulated in a single computational domain. This formulation has been implemented in the finite element code Cast3M. It combines ease of implementation and good performances in term of stability and convergence (LePotier et al., 1998) and allows to deal naturally with the saturated/unsaturated transition. In this work the soil water retention and relative permeability curves are given by van Genuchten models (van Genuchten, 1980):

Sw  Swr Se ðhÞ ¼ ¼ 1  Swr

(

1=½1 þ ðahÞn m

for h < 0

1

for h P 0

Þm 2 kr ðSe Þ ¼ Se1=2 ½1  ð1  S1=m e

ð4Þ ð5Þ

where Sw is the water saturation, Swr the water residual saturation, Se the effective saturation. a [L1] and n are the van Genuchten fitting parameters, with m ¼ 1  1=n. Some other models, like Brooks and Corey (1964), can also be used. Both Neumann (imposed flux) and Dirichlet (imposed head) boundary conditions (BC) can be imposed at the domain boundaries. Surface flow Runoff on land surface is usually described by depth averaged flow equations commonly referred as Saint Venant equations (e.g. in Beven, 2001). These equations consist of mass and momentum balance equations. The mass balance equation used here is

ohs ! ! þ r:ðhs U s Þ ¼ qs ot

ð6Þ

! where t is time [T], hs is the water depth [L] and U s is the depth 1 averaged flow velocity [LT ]. Index s stands for surface. qs represents the total sink and/or source of water [LT1] including rainfall, water infiltrating in the subsurface compartment or subsurface water exfiltrating from it. This equation is valid for shallow water and gentle slopes only, i.e. when water depth measured perpendicular to the slope may be approximated by the vertical water depth (Streeter et al., 2003). Following the diffusion wave approximation inertia terms of the momentum balance equation may be neglected, and the friction slope Sf ;i in i direction is therefore reduced to Sf ;i ¼ ri ðzl þ hs Þ (e.g. in Wasantha Lal, 1998), where Zl is the land surface elevation [L]. Many equations relates the velocity U s;i to the friction slope Sf ;i (Kirkby, 1978). In this study, we choose the Manning–Strickler uniform flow formula, which has been widely used to describe surface water flow:

U s;i ¼

2=3 hs

nman

qffiffiffiffiffiffi Sf ;i

ð7Þ

where nman [L1/3T] is the Manning’s coefficient, i stands for the x and y horizontal directions, the z axis being vertical upward. Assuming that the water depth gradient is much smaller than the surface elevation gradient, Hromadka et al. (1985) and Wasantha Lal (1998) have shown that Eq. (7) may be written as a function of mean local slope S as follows: 2=3

U s;i ¼ 

hs

pffiffiffi ri ðzl þ hs Þ nman S

ð8Þ

Introducing this equation in Eq. (6), the mass balance equation for surface flow is finally written in a similar way as in Gottardi and Venutelli (1993), Giammarco et al. (1996) or Panday and Huyakorn (2004): 5=3

! ohs ! hs pffiffiffi rðzl þ hs Þ  r: ot nman S

! ¼ qs

ð9Þ

11

This 2D nonlinear diffusion equation looks like a Richards equation pffiffiffi 5=3 with a horizontal hydraulic conductivity T s;xx ¼ T s;yy ¼ hs =nman S. As this equation is averaged along the vertical direction, these conductivities have the dimension of a transmissivity [L2T1]. Authors like Gottardi and Venutelli (1993) or Panday and Huyakorn (2004) express these conductivities as T s;ii ¼ hs ki where ki is called conductance. Eq. (9) shows that, in the diffusive wave approximation of Saint Venant equations, runoff flowrate is proportional to the runoff water hydraulic gradient. It suggests simulating runoff as a flow in a porous medium that has particular properties. Therefore, in our modelling approach, a 3D porous layer, called runoff layer, is added at the top of the subsurface domain to model surface flows (see Fig. 1). In order to define vertical water transfer in the runoff layer, one has to introduce a vertical hydraulic conductivity component K zz . The choice of this component must be guided by the two following physical principles: (i) when overland flow occurs, water pressure has to be hydrostatic along the vertical direction so that the pressure at the bottom of the runoff layer, i.e. at the land surface, can be set equal to the water depth hs and (ii) when there is no runoff, which means that water pressure along the land surface and water depth in the runoff layer are negative, all rainfall water has to infiltrate and vertical water transfer in the runoff layer has to be instantaneous. To fulfil these principles the vertical transfer along the runoff layer thickness must be much more rapid than both the vertical transfer in the soil and the runoff process itself. Therefore K zz must be much larger than the soil saturated hydraulic pffiffiffi conductivity K sat 2=3 and the horizontal conductance ki ¼ hs =nman S. Introducing the runoff layer thickness e [L], Eq. (9) may be rewritten in a 3D form:

! 1 ohs !  r :ðK s ðhs Þr ðzl þ hs ÞÞ ¼ 0 e ot

ð10Þ

where the components of the hydraulic conductivity tensor K s ðhs Þ are:

pffiffiffi 5=3 K s;xx ðhs Þ ¼ K s;yy ðhs Þ ¼ hs =ðe nman SÞ K s;zz ðhs Þ ¼ K zz

ð11Þ ð12Þ

and the apparent runoff velocity components:

U s;xx ðhs Þ ¼ K s;xx ðhs Þrx ðzl þ hs Þ

ð13Þ

U s;yy ðhs Þ ¼ K s;yy ðhs Þry ðzl þ hs Þ

ð14Þ

U s;zz ðhs Þ ¼ K s;zz rz ðzl þ hs Þ

ð15Þ

This velocity is not the depth averaged runoff velocity of Eq. (8), but the apparent runoff velocity in the runoff layer. The runoff flow rate [L2T1] components are

Q s;ii ¼ e U s;ii

ð16Þ

where i ¼ x; y. Let us remark that the vertical integration of the 3D equation Eq. (10) leads to Eq. (9). The source/sink term qs in Eq. (9) which accounts for rainfall and infiltration/exfiltration sources is removed from Eq. (10) as rainfall becomes a boundary condition at the top of the runoff layer and infiltration/exfiltration is naturally calculated, using fluxes and pressure continuity, at the land surface. The 3D runoff layer is a mathematical abstract construction which does not bring any new physical effects. Considering surface/subsurface coupled dynamics, one should naturally introduce at the land surface a 2D surface domain representing soil surface and mapping the 3D subsurface domain. One needs then to define 2D surface elements that allows to describe both horizontal transfer, i.e. runoff propagation, and vertical transfer, i.e. rainfall, infiltration and exfiltration. Finite elements, or even finite volumes, are not naturally suited for such an operation. This is why the initial 2D runoff model has been extended to a 3D model and that a

12

S. Weill et al. / Journal of Hydrology 366 (2009) 9–20

e Top of th

yer

runoff la

DIUM

OUS ME

OR YER = P OFF LA

RUN

High vertical permeability

ttom ace = bo

r

noff laye

of the ru

rf

Land Su

Pressure and fluxes continuity along the land surface

Subsurface domain

Fig. 1. The runoff layer concept.

high vertical conductivity was defined in the runoff layer. The introduction of a high hydraulic conductivity vertical component K zz allows to equalize water heads at the bottom and top of the runoff layer, and then to use a 3D computational domain to describe a 2D process. The same approach has been developed for the study of flow and transport in 2D fractures connected with 3D matrix blocks in 3D fractured media (Bath and Jackson, 2002). As for the subsurface flow Eq. (1), the runoff Eq. (10) is extended to negative values of pressure head, that is negative water depth. In this domain of negative pressure runoff must not occur. Therefore, by analogy with Richards equation for subsurface, we define for the whole pressure domain a water content hs :

 hs ðhs Þ ¼

for hs < 0

0

hs =e for hs P 0



0

for hs < 0 pffiffiffi 5=3 hs =ðe nman SÞ for hs P 0

K s;zz ðhs Þ ¼ K zz

ð18Þ ð19Þ

! oH !  r :ðK s ðhs Þ rHÞ ¼ 0 ot

ð20Þ

where

rs ðhs Þ ¼

ohs ¼ ohs



0 for hs < 0 1=e for hs P 0

ð21Þ

Following VanderKwaak (1999) or Panday and Huyakorn (2004), the water depth BC at the runoff layer outlet can be either a zero water depth gradient or a critical water depth. The first BC is equivalent to impose an outlet apparent runoff velocity equal to 5=3

U s;n ¼

hs

pffiffiffi rn ðzl Þ e nman S

~ and speWe define water content ~ h, hydraulic conductivity K ~ cific volumetric storativity r in the whole surface/subsurface domain as functions of water pressure head h in the subsurface domain and water depth hs in the surface domain:

ð17Þ

The water content in the runoff layer does not have any physical sense. The only meaningful variable which allows to quantify the volume of runoff water is the water depth hs . As the hydraulic head H in the runoff layer and the water depth hs are related by the relationship hs ¼ H  zl , the final runoff equation describing flow in the surface domain is written:

rs ðhs Þ

Generalized Darcy type equation for surface/subsurface flows

~hðHÞ ¼

and hydraulic conductivity components K s;xx ðhs Þ; K s;yy ðhs Þ and K s;zz ðhs Þ

K s;xx ðhs Þ ¼ K s;yy ðhs Þ ¼

for small scale processes such as water storage in microtopography. Like Panday and Huyakorn (2004) or VanderKwaak (1999), the runoff layer concept may be used to model streams and rivers.

ð22Þ

where index n stands for normal to the outlet and rn is the normal derivative. This zero water depth gradient BC is imposed at the outlet of the runoff layer. In this model, other water depth–velocity relationships, like Chezy or Darcy-Weisbach ones, may be implemented. Hydraulic conductivity of the runoff layer may be also modified to account

~ KðHÞ ¼



hðhÞ

in the subsurface domain

hs ðhs Þ in the runoff layer  KðhÞ in the subsurface domain

K s ðhs Þ in the runoff layer  r ðhÞ in the subsurface domain r~ ðHÞ ¼ rs ðhs Þ in the runoff layer

ð23Þ ð24Þ ð25Þ

Thus, a single Richards type equation with domain-dependent parameters describes surface and subsurface flows:

r~ ðHÞ

! oH ! ~ ~  r:ðKðHÞ r HÞ ¼ q ot

ð26Þ

~ is a global source/sink term equal to the source/sink term q where q of the subsurface domain (Eq. (1)). The combined use of this equation and of the mixed hybrid finite element formulation, described later in the paper, provides a continuous surface/subsurface coupling. As we said in the introduction, some authors have already modelled surface/subsurface interactions in an integrated way using either a first order coupling law or a matching boundary technique (VanderKwaak, 1999; Panday and Huyakorn, 2004; Kollet and Maxwell, 2006). Using a matching boundary technique, the pressure head along the land surface has to be checked during the simulation to switch the boundary condition from an imposed flux to an imposed head boundary when ponding occurs. The first order coupling law implemented in Panday’s or VanderKwaak’s works stipulates that the exchange flux between surface and subsurface is proportional to the head difference between the two domains. This method relies on the introduction of an interface connecting both domains and on a proportionality parameter. The coupling strategy presented here is in a way similar to the first order law strategy as it relies on controlling surface/subsurface dynamics through land surface flux without switching boundary condition at the land surface. Nevertheless, in our modelling approach, the total flow dynamics is fully described by a single Richards equation. Pressure and

S. Weill et al. / Journal of Hydrology 366 (2009) 9–20

fluxes continuity at the land surface is provided naturally. Pressure along the land surface is an ordinary unknown that, nonetheless, controls surface/subsurface interactions and rainfall/ runoff partition. No interface between surface and subsurface domains is introduced. Surface and subsurface waters interact directly at the physical land surface.

Space and time discretizations

Tracer transport equation

! ! ~ U ¼ KðHÞ rH oH ! ! r~ ðHÞ þ r  U ¼ q~ ot

An advective–dispersive transport equation has been implemented in the model. The aim is to follow rainfall or subsurface water and determine their relative contributions to a hydrograph. To reach this goal, we suppose that either rainfall or subsurface water is marked, i.e. contains a tracer concentration. The resolution of the transport equation gives access to the time evolutions of tracer concentrations and fluxes. In particular, the proportion of marked water that exit the system is directly connected with hydrograph separation issues. The transport equation in the subsurface domain can be written (e.g., De Marsily, 1986):

! oðhcÞ ! ! þ r :ð U c  D r cÞ ¼ qc ot

ð27Þ

! where h is the water content, c the concentration, U the Darcy 1 velocity [LT ], D the diffusion–dispersion tensor [L2T1] and qc the source/sink concentration term [T1]. The diffusion–dispersion tensor is split into a diffusive and a dispersive component:

! D ¼ d þ aj U j

ð28Þ 2 1

where d is the diffusive component [L T ], a the dispersivity tensor ! [L] and j U j the Darcy velocity norm. The transport equation to be solved in the surface domain is averaged along the vertical direction and can be written (VanderKwaak, 1999):

! oðhs cÞ ! ! þ r:ðQ s c  Ds r cÞ ¼ qcs ot

ð29Þ

! where hs is the water depth, c the depth averaged concentration, Q s 2 1 the runoff flow rate [L T ] which components are given in Eq. (16), Ds the surface diffusion–dispersion tensor [L3T1] and qcs the source/sink concentration term [T1] which accounts for tracer transfer between runoff and soil and eventually from rain to runoff. As for the subsurface domain, the diffusion–dispersion tensor can be written:

! Ds ¼ hs ds þ as jQ s j

ð30Þ 2 1

where ds is the diffusion tensor [L T ], as the surface dispersivity ! tensor [L] and jQ s j the runoff flow rate norm. ~ defined in Eq. (24) and (25), transport in the ~ Using h and K surface/subsurface continuum is described in a unified manner by the following equation:

~! oð~hcÞ ! ! þ r :ð U c  D r cÞ ¼ q~c ot

ð31Þ

! ! ~ ~ where U ¼ KðHÞ rH, and D and q~c equals respectively D and qc in the subsurface domain and Ds and qcs in the surface domain. Numerical solution In this section, the numerical schemes are first presented. Then, we focus on the type of boundary condition imposed at the runoff layer outlet. Finally details on convergence and time step strategy are given. A more complete description can be found in Weill (2007).

13

The single nonlinear Darcy Eq. (26) for flow is discretized with the mixed hybrid finite element (MHFE) formulation (Chavent and Jaffré, 1986; Mosé et al., 1994). This formulation solves simultaneously in each cell the following generalized Darcy equation and mass balance equation:

ð32Þ ð33Þ

Unknowns are defined at nodes different of those of the classical finite element (FE) formulation. In a cell the nodes are the centres of cell faces and the cell centre (see Fig. 2), whereas in FE they are the cell vertices. MHFE ensures mass conservation in each cell of the mesh and therefore is particularly interesting for the modelling of heterogeneous media. Moreover, fluxes and pressure are continuous at the interface between cells. This property is of particular importance for our modelling approach as it provides fluxes and pressure continuity at the land surface. The mixed finite elements used in Cast3M are the Raviart Thomas elements of lowest order ! (Brezzi and Fortin, 1991): in a given cell the velocity U is defined by the fluxes across the cell faces and the head pressure h is defined by its mean value in the cell, defined at the cell centre, and its mean values along each cell face, defined at the cell face centres. It has been pointed out recently (Le Potier, 2005; Younes and Fontaine, 2008) that accuracy problems may arise when using MHFE on skewed cells, or regular cells with highly anisotropic conductivity tensor. As studied slopes are gentle and the runoff layer has a constant thickness the runoff layer cells are not skewed and have parallel faces (see Fig. 2). It can be shown that, for this type of cell, the vertical transfer is practically decoupled from horizontal one as long as the anisotropy ratio K s;xx =K s;zz is much smaller than one. This is confirmed by the validation and verification test cases studied in the next section. The time discretization is implicit and the nonlinear terms are solved with an iterative Picard algorithm:

r~ ðhnþ1;i Þ !nþ1;iþ1 U

Hnþ1;iþ1  Hnþ1;0 ! !nþ1;iþ1 ~nþ1;iþ1 ¼q þr U Dt ~ nþ1;i Þ! ¼ Kðh rHnþ1;iþ1

ð34Þ ð35Þ

~nþ1;iþ1 the where n is the time index, i the iteration index and q source/sink term. Though Picard iterative algorithm leads to a slower convergence than the Newton algorithm (Paniconi and Putti, 1994), our experience is that Picard method has been efficient to solve many nonlinear problems: unsaturated flow, multiphase flow, thermo-hydraulics or hydro-mechanics (LePotier et al., 1998; Genty et al., 2000; Barnel et al., 2002). As transport in runoff is mainly advective the transport equation is solved with a Finite Volume (FV) formulation. It has been shown that for this type of problem the FV implemented in Cast3M is more accurate than the MHFE (Michel et al., 2004). The FV formulation solves simultaneously in each cell the following mass flux and mass balance equations:

! ~! J ¼ D rc ~ ! ! ~ oðhcÞ þ r  J ¼ Qc ot ! ! ~c ¼ q ~c  r  ð U cÞ Q

ð36Þ ð37Þ ð38Þ

Eqs. (36) and (37) are discretized with a cell centred FV scheme proposed initially by Aavtsmark in 2D geometry (Aavstmark et al., 1998) and extended by Le Potier (2005) to 3D geometries, unstructured meshes and highly anisotropic diffusion operators. In the scheme improved by Le Potier the gradient on each cell

14

S. Weill et al. / Journal of Hydrology 366 (2009) 9–20

r

ff laye

e runo

f th Top o

Water depth profile at iteration i

Outlet

runoff of the Bottom nd surface la layer =

o o

Calculation nodes in MHFE Fig. 2. Schematic illustration of the iterative boundary condition concept (time step (n + 1)).

vertex is computed using the cell centred variable c and its values on the cell edges; then these last ones are eliminated by imposing flux continuity condition and the resulting global matrix is symetric definite positive. The method appears to be quite robust and precise (Michel et al., 2004). After integration in a cell the advec!! ~ c , Eq. (38), is expressed as tive term r ð U cÞ of the source term Q a sum of the advective fluxes across each face. Depending on the cell Peclet number Pe this sum is discretized with an upwind scheme ðPe > 2Þ, a centred scheme ðPe < 2Þ, or a mixed one called ‘‘upwicent” (Michel et al., 2004). Time discretization is always implicit. The flow and transport equations are solved sequentially. First, we solve the flow equation to obtain the water content and the velocity fields. Then, we solve the transport equation. For time step ðn þ 1Þ, the semi-discrete equation solved for transport can be written:

~hnþ1 cnþ1  ~hn cn ! ~! r cnþ1 Þ ¼ q~c nþ1 þ r  ð U nþ1 cnþ1 Þ  ðD Dt

ð39Þ

All the unknowns in this equation, except the advective flux, are defined at the mesh centres. The model and these algorithms have been implemented in the finite element code Cast3M (website www-cast3m.cea.fr). It is a general multi purpose computational tool developed, since two decades at the Commissariat à l’Energie Atomique (CEA), for solid and fluid mechanics applications. The code has been widely used in hydrology and hydrogeology and benchmarked with many other codes (e.g., Bath and Jackson, 2002; Michel et al., 2004; Rutqvist et al., 2005).

Nevertheless, the water depth gradient is not a natural variable in MHFE. The different attempts to linearize Eq. (40) in the Picard scheme have shown to fail, the results were quite unstable. Our experience in MHFE is that a Dirichlet BC is preferable in highly non linear problems. Therefore, the zero water depth gradient BC was approximated by an iterative water depth BC. At the centre of an outlet face (see Fig. 2), the imposed BC at time step (n + 1), iteration (i + 1) is

(

nþ1;iþ1

hs;o

nþ1;iþ1 Q s;o

nþ1;i

¼ hs;c ¼0

nþ1;i

>0

nþ1;i hs;c

<0

if hs;c if

ð41Þ

where Q s;o is the outlet flowrate defined at the centre of the outlet face, hs;o is the imposed water depth at the outlet and hs;c is the mean water depth of the cell, defined at the cell centre. It can be shown that when the Picard iterative algorithm converges this BC becomes equivalent to a zero water depth gradient boundary condition (Weill, 2007). Nevertheless, as the Picard algorithm is stopped when a convergence criterion is reached a non zero water depth gradient, proportional to the convergence criterion, still remains after convergence. Consequently we expect the outlet flux to be overestimated during rising limbs and underestimated during receding limbs. Indeed, during rising limbs, a small but positive water depth gradient is created and a flux greater than expected escapes from the runoff layer outlet. The inverse effect is expected to happen during receding limbs. Apart from this small correction this iterative BC allows the water depth to vary from one step to another as if the outlet was ‘‘transparent”, i.e. water flux guided only by gravity.

Boundary conditions Convergence and time step strategy The flow BC may be Dirichlet type (imposed head), Neumann type (imposed flux), or both. For Dirichlet BC the head is imposed at the centres of the cell faces describing the concerned boundary domain. For the Neumann BC it is the flux which is imposed. The same holds for transport BC. For Dirichlet type concentration is imposed at the cell face centres and for Neumann type the diffusive dispersive flux is imposed. At the outlet of the runoff layer, we impose the zero water depth gradient boundary condition (BC) given by

rn ðhs Þ ¼ 0

ð40Þ

where index n stands for normal to the outlet and rn is the normal derivative. This BC is equivalent to the one defined in Eq. (22).

The studied system is highly nonlinear; it shows some similarity with free boundary problems in fluid mechanics (Crank, 1988). For instance a cell of the runoff layer may commute during a time step from a state with no water (capillary capacity and hydraulic conductivity equal to zero) to a state with water (capillary capacity equal to one and non zero hydraulic conductivity). For this type of problem, experience shows that underrelaxation may be needed to make the system numerically stable, see for instance (Huyakorn et al., 1996) for a discussion of this technique. Underrelaxation means that for a current iteration i the capacity and the conductivity terms are not only calculated from the pressure obtained at the previous iteration ði  1Þ but also from that of iteration ði  2Þ:

15

S. Weill et al. / Journal of Hydrology 366 (2009) 9–20

ð42Þ ð43Þ

with a1 þ a2 ¼ 1. Note that not only the water pressure may be underrelaxed but also the coefficients themselves. They are evaluated with the water pressure defined at the cell centre, which is a natural choice in MHFE, but other choices are possible (Belfort and Lehman, 2005). The time stepping strategy adopted in Cast3M is classical: time step grows, up to a maximum, as long as the number of iterations required for convergence is low, typically smaller than ten. When the number of iterations exceed this number time step is divided by two. Experience shows that convergence is easier (greater time steps can be used) in subsurface (soil alone) or surface (impermeable soil) systems than in coupled surface–subsurface systems. In surface–subsurface systems when the extension of the saturated part of the runoff layer varies with time, i.e. number of saturated cells increases or decreases, the time step may be reduced, up to a factor ten, as compared to the time step used when the extension is permanent. More precisely if the number of saturated cells is increasing (or decreasing), this process must be captured cell by cell. When one cell is, for instance, about to be saturated, that is n when hs is commuting from a negative value hs < 0 to a positive nþ1 value hs > 0 the transition must be captured in a way that this cell only is concerned by saturation during iterations. This imposes a small time step. When the time step is too large the number of saturated cells becomes oscillatory with iterations and convergence is not reached. Depending on the studied system different convergence criteria may be used. The classical one is defined as

P

nþ1;iþ1

cells jh

nþ1;i

nþ1;i

h

n

h j

j

<

Comparison with 1D analytical solution of the kinematic wave approximation Two approximations of the full St. Venant equations are often used to model overland flow (Beven, 2001): the kinematic and the diffusive wave approximations. Nevertheless, only the kinematic wave equation has an analytical solution (Eagleson, 1970). The overland flow model implemented in Cast3M is then compared to the analytical solution of the kinematic wave equation on a simple 1D test case (Kazezyilmaz-Alhan et al., 2004) to assess the performance of the overland flow model developed. A 30 min and 1; 4  105 ms1 intensity rainfall is simulated over a 183 m long parking lot. The Manning’s coefficient value is 0.025 sm1/3 and the slope value 0.0016. The sensibility to spatial and time discretizations is first examined: a mesh of either 10 or 100 cells is considered with constant time steps of 1 s or 100 s. Fig. 3 shows the simulated and the analytical kinematic wave discharge. The agreement between Cast3M results and the analytical solution for the rising and receding limbs is good. As expected, the hydrograph peak is poorly represented. The diffusive wave curve approaches the peak more smoothly, due to the diffusive term (i.e. the head pressure gradient) added to the kinematic wave equation to obtain the diffusive wave equation. This agreement shows that the introduction of a high vertical conductivity K zz does

3.

ð44Þ

where the sum is performed on all the cells of the discretized domain. Following the preceding discussion the convergence criterion imposed on hs is much more drastic in surface–subsurface systems, typically  ¼ 105 , than in subsurface or surface systems,  ¼ 103 . At last, when the time step decreases underrelaxation makes convergence easier and is often necessary. As already emphasized by different authors (e.g., Panday and Huyakorn, 2004), convergence is very sensitive to soil water retention and relative permeability curves. There again experience shows that for a given model, van Genuchten for instance, a little change in soil parameters values may completely change convergence.

Analytical solution of the kinematic wave equation Simulated outflow with 10 meshes: dt=100s Simulated outflow with 10 meshes: dt=1s

2.5 3 -1

cells jh

Overland flow model validation

2.

-3

P

area at the toe of a 3D virtual hillslope during a rainfall event is presented.

Flowrate (10 m s )

rnþ1;i ðhÞ ¼ rða1 hnþ1;i1 þ a2 hnþ1;i2 Þ nþ1;i nþ1;i1 nþ1;i2 ðhÞ ¼ kða1 h þ a2 h Þ k

1.5 1. 0.5 0

0

10

20

Validation, verification, and application examples

40

50

60

70

Analytical solution of the kinematic wave equation Simulated outflow with dt = 1 s and 10 meshes Simulated outflow with dt = 1s and 100 meshes

2.5

-3

3 -1

Flowrate (10 m s )

3.

The numerical validation of integrated surface/subsurface models, and in particular surface/subsurface coupling, is a difficult task. Indeed, none of the classical hydrological situations illustrating the coupling have an analytical solution. Only elementary processes like runoff or infiltration in simple geometries may be described analytically. The only integrated surface/subsurface examples available in the literature are either numerical or experimental. This turns any tentative of model validation into a model intercomparison or a model verification. The approach adopted here is: first validate the overland flow model, then verify the integrated model on an experiment, and finally present some coupled application examples. To validate the overland flow model, the results are compared to an analytical solution of the kinematic wave equation and to results from other overland flow models of the literature. The verification of the integrated flow and transport model is performed with the Abdul and Gillham laboratory sandbox experiment (Abdul and Gillham, 1984). Finally, an illustration of the model ability to simulate the dynamics of a contributed saturated

30

Time (min)

2.

1.5

1. 15

20

25

30

35

40

Time (min) Fig. 3. Comparison between Cast3M results and the analytical solution of the kinematic wave equation for different temporal and spatial discretizations.

16

S. Weill et al. / Journal of Hydrology 366 (2009) 9–20

not perturb the runoff dynamics along the surface. Concerning the impact of spatial and time discretizations on the simulated hydrograph, it seems that the spatial one does not greatly influence results. Nevertheless, when time steps increase, the precision of the approach decreases. This is linked to the iterative boundary condition which needs small time steps to follow correctly the water depth evolution at the runoff layer outlet. The 2D tilted V-catchment test case The problem considered here is overland flow generated by a rainfall event on a 2D tilted V-catchment (see Fig. 4). This test case was first published by Giammarco et al. (1996) and stands today as a reference to validate runoff models (e.g., VanderKwaak, 1999; Kollet and Maxwell, 2006). The system is composed of two 1000 m  800 m slopes connected in the middle to a 1000 m length channel 20 m wide. Only one half of the domain is considered due to symmetry. Surface slopes are respectively 0.05 perpendicular to the channel and 0.02 parallel to the channel. Manning’s coefficients values are 0:015 sm1=3 for the slope and 0:15 sm1=3 for the channel. A rainfall rate of 3:106 ms1 is applied on the whole surface for 90 min with a total simulation time of 180 min. A constant spatial discretization with Dx ¼ Dy ¼ 20 m and constant time steps of 10 s are used. Fig. 5 compares the simulated discharge to predictions of Di Giammarco (Giammarco et al., 1996) and Panday and Huyakorn with their model MDMH (Panday and Huyakorn, 2004). Cast3M results are obtained with the itera-

l ne an

1000 m

Ch

Sl

op

e=

0.0

2

Symmetry axis

Slope = 0.05

20 m

800 m

Fig. 4. Three-dimensional view of the tilted V-catchment area (Giammarco et al., 1996).

6

3 -1

Verification of the integrated flow and transport models In most of previous works on integrated surface/subsurface modelling (VanderKwaak, 1999; Kollet and Maxwell, 2006), the experimental system presented by Abdul and Gillham (1984) is the only system used for validation. Indeed, few experimental systems are available in the literature to validate integrated surface/ subsurface models. The Abdul and Gillham system was designed to study and demonstrate the rapid capillary zone response first described by Sklash and Farvolden (1979). The experiment consists in a Plexiglas sandbox 140 cm long, 120 cm high and 8 cm wide (see Fig. 6). The outlet of the system is located at a height of 76 cm. The box was filled with medium-fine sand. The surface slope is uniform equal to 12° and Manning parameter equals 0.185 sm1/3. In one of their experiments, the initial condition is hydrostatic and water table is located at the toe of the slope. A 1:2  105 ms1 constant rainfall rate is applied for 20 min on the whole surface domain. Recharge water contains bromide so that to identify the contribution of rainfall water to the total discharge. Discharge volume and concentrations were measured for 25 min. Abdul and Gillham also modelled their experimental system with a simple 2D unsaturated–saturated model with a matching boundary algorithm to be able to follow the evolution of the seepage face. The primary drainage and wetting curves of the sand used in this experiment are presented in Abdul and Gillham’s article (Abdul and Gillham, 1984). The saturated hydraulic conductivity is supposed to be 3:5  105 ms1 and the total porosity 0.34. The saturation–capillary pressure relationship was approximated using

lator

ll simu

Simulated hydrograph with our modelling approach Simulated hydrograph from Digiammarco et al, 1996 Simulated hydrograph using MDMH with CD Simulated hydrograph using MDMH with ZDG

5

Flowrate (m s )

tive BC (Eq. 41), those of Di Giammarco with a critical depth boundary and Panday and Huyakorn could simulate both a critical depth or a zero depth gradient boundary. The hydrograph predicted using the present model shows a general good agreement for the entire hydrograph with the one predicted by Di Giammarco. The model fits well the rising and receding limbs and the time to peak discharge and the peak discharge value are correct. The simulated hydrograph is also bounded by the results obtained by Panday and Huyakorn with their MDMH model using either critical depth or zero depth gradient BC. By comparing our result with those obtained with MDMH and a zero depth gradient BC, we observe that, as expected, the outlet flux is overestimated during the rising limb and underestimated during the receding one. Nevertheless, these results show that the 3D overland flow model implemented in Cast3M with its iterative boundary condition and its highly anisotropic hydraulic conductivity tensor represent runoff dynamics properly.

Rainfa

4

)

º slope

ace (12

urf Land s

Initial water table

3 2 0.76 m

1 0 0

No flow boundaries

20

40

60

80

100

120

140

160

180

Time (min) 1.4 m Fig. 5. Comparison between simulated hydrograph and results from other hydrological models (CD – critical depth and ZDG – zero depth gradient).

Fig. 6. The Abdul and Gillham system.

17

S. Weill et al. / Journal of Hydrology 366 (2009) 9–20

Total simulated Event water simulated Pre-event water simulated

Total measured Event water measured Pre-event measured

1 Normalized flux (-)

the van Genuchten model (see Eqs. (4) and (5)) and the following n and a values were obtained for the wetting curve: n = 5.5 and a = 2.3 m1. This sand has a capillary fringe approximately 30 cm thick so that initially the whole domain is nearly fully saturated. The domain was discretized into 50 layers and 50 rows which defined a centimetric cell. To determine relative contributions of preevent and event water, rainfall water was tracked. Simulation rainfall was assumed to contain a tracer concentration equal to unity. The diffusion coefficient was set to 109 m2 s1 and the dispersivity to 0.01 m. No flow boundaries were applied at the bottom and vertical sides of the domain so that the only exchange surface of the domain is ground surface (see Fig. 6). The time step was 0.1 s at the beginning of simulation and increased up to a value of one second at the end. All the fluxes presented in the following results are normalized by the rainfall flux imposed at the top of the runoff layer. The evolution of the normalized flux along the land surface (entering fluxes are negative by convention), displayed on Fig. 7, shows that the model implemented in Cast3M is able to describe the three surface regimes observed and simulated by Abdul and Gillham: at the top of the slope, all rainfall water infiltrates (normalized flux equals 1); in the middle of the slope, there is a partition of rainfall water between a part which infiltrates and another which produces overland (normalized flux between 1 and 0); at the toe of the slope, water exfiltrates (normalized flux positive). Beaugendre et al. (2006), VanderKwaak (1999) and Kollet and Maxwell (2006) have already simulated this particular flow dynamics and presented such results. In Beaugendre’s work (Beaugendre et al., 2006), the influence of the type of soil (sand, silty-loam or clay) on the flow regime, and especially the influence of the capillary fringe length, is investigated. They show that the relative hydraulic conductivity and water content parameters have a great influence on the flow dynamics of the system (equilibrium time, steady state saturated length, etc.). The following part focuses on the hydrograph separation that was only discussed by VanderKwaak (1999) and Cloke et al. (2006). Fig. 8 displays the evolution with time of different simulated and measured water fluxes. It shows that, for the total outflow, the simulated time of concentration is shorter than the experimental one. The lack of air phase in the modelling approach could be responsible for that. Indeed, air phase compression increases water storage capacity at early time and the presence of air could significantly slow down infiltration process (Touma and Vauclin, 1986). This may explain the delay in the development of the seepage face. Concerning the relative contributions of preevent and event water, the measured and simulated curves do not fit well. As these contributions are closely related to the global

0.8 0.6 0.4 0.2 0 0

5

10

15

20

25

Time (minutes) Fig. 8. Simulated vs. measured normalized total outflow, pre-event water outflow and event water outflow.

flow dynamics, the direct comparisons are hampered due to inaccuracies in the early simulated total outflow. Nevertheless, our transport model can properly describe the global trends. Concerning the evolution of pre-event water, Fig. 8 shows that at early time its contribution to the total hydrograph is large, which is in good agreement with the experimental measurements. Then, this contribution decreases with time as rainfall replaces pre-event water through infiltration process. The rainfall concentration field presented Fig. 9 shows that processes affecting transport take place in a thin region near land surface. In comparison, the total head field is affected by rainfall in the entire system. We believe that the difference between flow and transport dynamics could explain the high proportion of pre-event water in the total hydrograph for early time. Indeed, flow kinetics is much faster than transport one. As rainfall begins and event water infiltrates, total head gradients appear nearly instantaneously in the whole domain as it is initially nearly fully saturated. With the two no-flow lateral boundaries considered, it rapidly makes pre-event water exfiltrate at the toe of the slope. Mixing in the thin layer near land surface makes event water participate to total outflow but the main part of the water exiting the system is exfiltrating pre-event water. The influence of dispersivity on pre-event and event water contributions to total hydrograph was also studied. The relative contributions of pre-event and event water to total hydrograph presented Fig. 8 is obtained with a dispersivity of 0.01 m that gives the best fit concerning pre-event contribution peak. During the fitting procedure, it was shown, as discussed by Jones et al. (2006), that dispersivity has a great influence on relative pre-event and

Normalized land surface flux (-)

3. 2.5 2.

Exfiltration

Infiltration and runoff

1.5 1. 0.5 0. −0.5 −1.

Infiltration only

−1.5

0.

1.2 1. 0.2 0.4 0.6 0.8 Distance from the toe of the slope (m) Fig. 7. Normalized flux along the land surface.

1.4

0. 0.05 0.10 0.15 0.20 0.24 0.29 0.34 0.38 0.43 0.48 0.52 0.57 0.62 0.66 0.71 0.76 0.80 0.85 0.90 0.94 1. Fig. 9. Concentration field at the end of the rainfall event (t ¼ 1200 s).

18

S. Weill et al. / Journal of Hydrology 366 (2009) 9–20

event water contributions. The higher the dispersivity, the deeper the extension of the mixing zone near surface. As a consequence, the event water concentration near surface and the event water contribution to total hydrograph decrease when dispersivity increases. Flow and transport fluxes at the system outlet appeared to be quite insensitive to the choice of runoff laws, Manning or DarcyWeisbach, or runoff parameters, Manning coefficient for instance. It seems that, at the scale considered here, which is a small scale, those laws and parameters do not influence the hydrological response of the system, and especially the flux distribution along the land surface. The Abdul and Gillham system seems then to be a subsurface controlled system. Indeed, from a mathematical point of view the subsurface flow regime is governed by a Richards equation with a moving BC along land surface. It is clear that water pressure distribution in subsurface does not depend on runoff height, this one may be set equal to zero, which is equivalent to set Manning coefficient to a very high value, simulating instantaneous transfer of runoff water. The extension of seepage face with time is governed by a Dunne mechanism, i.e. runoff occurs when land surface becomes completely saturated, governed by infiltration only. This conclusion is closely linked to the small size of the system. It can not be generalized to other hydrological systems. Finally a sensitivity analysis to the spatial discretization was also performed with a loose mesh, 10 layers and 10 rows, and a more refined one, 100 layers and 100 rows, and adapted time steppings (Weill, 2007). Flow results show that the flux distributions along the soil surface and the hydrographs are very similar.

cretized with 40 meshes in the x and y directions and with 10 meshes in the vertical direction. The iterative boundary condition presented previously is imposed at the outlet of the channel. No flow boundaries are applied on the rest of the domain. Channel flow is modelled using the same equation as for overland flow. The aim of this section is to show that the model developed can deal with 3D situations. Fig. 11 displays the evolution with time of the normalized outflow and water depth at the outlet of the system. The maximum water depth reached during the simulation is around 4 cm and the maximum outflow is 5  104 m3 s1 . The shape of the hydrograph is classical and similar to those presented by Govindaraju and Kavvas (1991). As expected, water depth evolution at the outlet follows the hydrograph. The evolution with time of the normalized saturated area (saturated area divided by the total area of the surface domain) is plotted on Fig. 12. It is interesting to note that the rising limb of the normalized saturated is convex and shows, at a specific time, a sharp transition to a plateau (Beaugendre et al., 2006; Weill, 2007). This sharp transition is linked to the dimensions of our system and to the fact that the whole land surface becomes saturated. For a bigger system and as shown by Govindaraju and Kavvas (1991), the evolution of the saturated area tends smoothly to the permanent state. The runoff velocity field at the end of the rainfall is presented Fig. 13. As expected, this velocity field follows the greatest slope. One can see that the flow is parallel to the channel in the very

1

3D saturated contributing area dynamics

10 m

0.25m

Normalized value (-)

0.8

0.6 Outflow evolution

0.4

0 0

10 m

Fig. 10. Schematic of the three-dimensional virtual hillslope considered.

10

20

30

40 50 Time (min)

60

70

80

90

Fig. 11. Normalized outflow hydrograph and water depth at the outlet as a function of time.

100

80

60

40

20

0

0.1m

Water depth evolution

0.2

Normalised saturated area (%)

Slope = 0.05

Ch

an

ne

lS

lop

e=

0.0

05

The saturated contributing areas evolution near streams and rivers is of particular interest when studying runoff generation processes. Overland flow on these areas could be responsible for a major part of the hydrograph when the existence of a shallow water table is exhibited (Beven, 2001). A 3D virtual hillslope, that has the same geometry as the one studied by Govindaraju and Kavvas (1991), is simulated. Only the toe of the slope, where water table reaches land surface and saturated area appears, is considered. The system is 10 m long, 10 m wide and less than 1 m deep (see Fig. 10). Slope values along the channel direction and perpendicular to the channel are respectively 0.005 and 0.05. The soil’s relative hydraulic conductivity and water content functions were estimated from those published by Govindaraju. Saturated hydraulic conductivity value is set to 5  105 ms1 and specific storage is 104 m1 . The Manning coefficients are respectively 0.1 and 0.01 sm1/3 in the channel and in the slope. A 60 min and 5  106 ms1 intensity rainfall with a 30 min recession period is simulated. Initially a constant head is fixed and the water table is located at the toe of the slope. The computational domain is dis-

0

10

20

30

40 50 Time (min)

60

70

80

Fig. 12. Normalized saturated area evolution with time.

90

S. Weill et al. / Journal of Hydrology 366 (2009) 9–20

Channel Fig. 13. Normalized velocity field in the surface domain for t = 60 min.

low part of the domain. In his article, Govindaraju assumed that flows along the channel direction can be neglected and consequently models 3D hillslopes as a succession of coupled 2D hillslopes. Regarding the velocity field presented here, this approach seems simplistic to describe accurately the evolution with time of saturated contributing areas, at least at the spatial scale investigated here. The study performed by Govindaraju and Kavvas (1991) on the influence of different parameters (geometry, soil type,...) on saturated area evolutions is very instructive but fully coupled 3D simulations, as presented in this section, may be needed to further extend our understanding of runoff generation processes. Conclusion A new modelling approach to deal with surface/subsurface flow and transport has been presented. The implemented model couples in a fully integrated way a 3D variably saturated flow equation with a 2D diffusive wave approximation for overland flow. The surface equation is treated as a Richards equation. A 3D layer, called runoff layer, is then introduced at the top of the computational domain to simulate surface processes. The whole dynamics is described in a Darcian continuum, from saturated zone to land surface, where all the water fluxes are proportional to the hydraulic gradient. Coupling between surface and subsurface is done via a tensorial hydraulic conductivity which provides fluxes and pressure continuity at the land surface, without using the two traditional methods that are the first order law coupling and the matching boundary condition. An advection–dispersion equation has also been implemented to deal with hydrograph separation and to be able to determine the relative contributions of pre-event and event water. The governing equations are solved within the finite element code Cast3M using numerical techniques designed for highly nonlinear problems. The overland flow model was validated with an analytical solution and previously published results. The integrated surface/subsurface flow and transport model was then evaluated using a sandbox experiment and showed a relative good agreement. As a matter of fact, if the general characteristics (shape and values) of the experimental water fluxes (total, pre-event, event) are well reproduced, simulated kinetics is much faster than the experimental one. Dispersivity appears to be a key parameter when modelling hydrograph separation with an advection–dispersion model. The dynamics of the toe of a 3D virtual hillslope was finally simulated to illustrate the applicability of the model implemented in Cast3M. It was shown that fully coupled 3D sim-

19

ulations are needed to properly describe the complex coupled surface/subsurface behaviour. The model implemented in Cast3M seems rather competitive. As for many other codes, the time step has to be small (sometimes less than one second) when the hydraulic response of the system is fast and abrupt: at the beginning or at the end of a rainfall event when the area of the contributive zone increases or decreases. Once the dynamics is settled, larger time steps (up to 1 min or more, depending on the cell size) can be used. Therefore, some hydrological simulations may require up to 105 iterations. In conclusion, these results show that the model presented in the paper works well. Nevertheless, there is a need for more benchmarking between codes and models. This could help for a better qualification and a better understanding of the different approaches that couple in an integrated way surface and subsurface flows. This type of integrated approach may be useful to study the hydrological processes at the plot scale (1–100 m2) or at the scale of a small watershed (VanderKwaak and Loague, 2001). The application to the catchment scale (more than 1 km2) demands important computer resources and constitute a challenge for the future. References Aavstmark, I., Boe, O., Mannseth, T., 1998. Discretization on unstrucutred grids for inhomogeneous anisotropic media. Part I: derivation of the methods. SIAM J. Sci. Comput. 19 (5), 1700–1716. Abbott, M.B., Bathurst, J.C., Cunge, J.A., O’Connell, P.E., Rasmussen, J., 1986a. An introduction to the European hydrological system – Système Hydrologique Européen, ‘‘SHE”, 1: history and philosophy of a physically based, distributed modelling system. J. Hydrol. 87, 45–59. Abbott, M.B., Bathurst, J.C., Cunge, J.A., O’Connell, P.E., Rasmussen, J., 1986b. An introduction to the European Hydrological System – Système Hydrologique Européen, ‘‘SHE”, 2: structure of a physically based, distributed modelling system. J. Hydrol. 87, 61–77. Abdul, A.S., Gillham, R.W., 1984. Laboratory studies of the effects of the capillary fringe on streamflow generation. Water Resour. Res. 20 (6), 691–698. Ambroise, B., 1995. Topography and the water cycle in a temperate middle mountain environment: the need for interdisciplinary experiments. Agr. Forest Meteorol. 73, 217–235. Anderson, M.G., Burt, T.P., 1978. The role of topography in controlling throughflow generation. Earth Surf. Process. 3, 331–344. Barnel, N., Le Potier, C., Maugis, P., Mouche, E., 2002. Impact of a thermal radioactive waste on the thermal-hydraulic behaviour of a clay engineered barrier system. Poro-mechanics, vol. 2. Balkema Publishers. Barros, A.P., Knapton, D., Wang, M.C., Kuo, C.Y., 1999. Runoff in shallow soils under laboratory conditions. J. Hydraul. Eng. 4 (1), 28–37. Bath, A.H., Jackson, C.P., 2002. Task force on modelling of groundwater flow and € Hard Rock Laboratory, Review of Task 5, SKB transport of solutes – Aspo International Progress Report (IPR-03-10), SKB Stockholm. Bear, J., Bachmat, Y., 1991. Introduction to Modeling of Transport Phenomena in Porous Media. Kluwer Academic Publishers. Beaugendre, H., Ern, A., Esclaffer, T., Gaume, E., Ginzburg, I., Kao, C., 2006. A seepage face model for the interaction of shallow water table with ground surface: Application of the obstacle-type method. Adv. Water Resour. 329, 258–273. Belfort, B., Lehman, F., 2005. Comparison of equivalent conductivities for numerical simulation of one-dimensional unsaturated flow. Vadose Zone J. 4, 1191–1200. Bernard Michel, G., Le Potier, C., Beccantini, A., Gounand, S., Chraibi, M., 2004. The ANDRA Couplex 1 test case: comparison between finite element, mixed hybrid finite element and finite volume discretizations. Comput. Geosci. 8, 187–201. Betson, R.P., 1964. What is watershed runoff? J. Geophys. Res. 69, 1541–1551. Beven, J.B., 2001. Rainfall–Runoff Modelling. Wiley, New York. Brezzi, F., Fortin, M., 1991. Mixed and hybrid finite methods. Springer-Verlag, New York. Brooks, R.H., Corey, A.T., 1964. Hydraulic properties of porous media. Hydrol. Pap., Colo. State Univ., Fort Collins 3, 27. Cappus, P., 1960. Bassin experimental d’Alrance: étude des lois de l’écoulement Application au calcul et à la prévision des débits. La Houille Blanche A, 493–514. Celia, M.A., Bouloutas, E.T., Zarba, E.L., 1990. A general mass-conservative numerical solution for the unsaturated flow equation. Water Resour. Res. 26, 1483–1496. Chavent, G., Jaffré, J., 1986. Mathematical models and finite elements for reservoir simulation. Studies in Mathematics and its Applications, vol. 17. Elsevier. Cloke, H., Anderson, M.G., McDonnell, J.J., Renaud, J., 2006. Using numerical modelling to evaluate the capillary fringe groundwater ridging hypothesis of streamflow generation. J. Hydrol. 316, 141–162. Crank, J., 1988. Free and Moving Boundary Problems. Clarendon Press. Dabbene F., Mixed hybrid finite elements for transport of polluants by underground water. In: Proceedings of the 10th International Conference on Finite Element in Fluids, Tucson, Arizona, USA, January 5–8, 1998.

20

S. Weill et al. / Journal of Hydrology 366 (2009) 9–20

De Marsily, G., 1986. Quantitative Hydrology. Academic Press, Harcourt Brace Jovanovich, New York. Di Giammarco, P., Todini, E., Lamberti, P., 1996. A conservative finite element approach to overland flow: the control volume finite element formulation. J. Hydrol. 175, 267–291. Dunne, T., Black, R.D., 1970. Partial area contributions to storm runoff in a small new england watershed. Water Resour. Res. 6, 1296–1311. Dunne, T., Zhang, W., Aubry, B.F., 1991. Effects of rainfall, vegetation, and microtopography on infiltration and runoff. Water Resour. Res. 27, 2271–2285. Eagleson, P.S., 1970. Dynamic Hydrology. McGraw-Hill, New York. Esteves, M., Faucher, X., Galle, S., Vauclin, M., 2000. Overland flow and infiltration modelling for small plots during unsteady rain: numerical results versus observed values. J. Hydrol. 228, 265–282. Freeze, R.A., 1972. Role of subsurface flow in generating surface runoff 1. Base flow contributions to channel flow. Water Resour. Res. 8 (3), 609–623. Freeze, R.A., 1972. Role of subsurface flow in generating surface runoff 2. Upstream source areas. Water Resour. Res. 8, 1272–1283. Freeze, R.A., Harlan, R.L., 1969. Blue-print for a physically-based digitally simulated hydrologic response model. J. Hydrol. 9, 237–258. Genty, A., Le Potier, C., Renard, P., 2000. Two-phase flow upscaling with heterogeneous tensorial relative permeability. Computational Methods in Water Resources XIII. Computational Mechanics Publications. Gimenez, R., Govers, G., 2001. Interactions between roughness and flow hydraulics in eroding rills. Water Resour. Res. 37, 791–799. Gottardi, G., Venutelli, M., 1993. A control-volume finite element model for twodimensional overland flow. Adv. Water Resour. 16, 227–284. Govindaraju, R.S., Kavvas, M.L., 1991. Dynamics of moving boundary overland flows over infiltrating surfaces at hillslopes. Water Resour. Res. 27, 1885–1898. Horton, R.E., 1993. The role of infiltration in the hydrologic cycle. Eos Trans. AGU 14, 446–460. Hromadka II, T.V., Lai, Chintu, 1985. Solving the two-dimensional diffusion flow model. In: Proceedings: ASCE Hydraulics Division Specialty Conference, Orlando, Florida, pp. 555–562. Huyakorn, P.S., Springer, E.P., Guvanasen, V., Wadsworth, T.D., 1996. A threedimensional finite-element model for simulating water flow in variably saturated porous media. Water Resour. Res. 22, 1790–1808. Jones, J.P., Sudicky, E.A., Brookfield, A.E., Park, Y.-J., 2006. An assessment of the tracer-based approach to quantifying groundwater contributions to streamflow. Water Resour. Res. 42, W02407. doi:10.1029/2005WR004130. Kazezyilmaz-Alhan, C.M., Medina Jr, M.A., Rao, P., 2004. On numerical modeling of overland flow. Appl. Math. Comput.. doi:10.1016/j.amc.2004.06.063. Kirchner, J.W., 2003. A double paradox in catchment hydrology and geochemistry. Hydrol. Process. 17, 871–874. Kirkby, M.J., 1978. Hillslope Hydrology. Wiley, New York. Kollet, S.J., Maxwell, R.M., 2006. Integrated surface–groundwater flow modelling: a free-surface overland flow boundary condition in a parallel groundwater model. Adv. Water Resour. 29, 945–958. Le Potier, C., 2005. Schémas volumes finis pour des opérateurs de diffusion fortement anisotropes sur des maillages non structurés. C.R. Acad. Paris, Ser. I 340, 921–926. Le Potier, C., Mouche, E., Genty, A., Benet, L.V., Plas, F., 1998. Mixed hybrid finite element formulation for water flow in unsaturated porous media. Computational Methods in Water Ressources XII. Computational Mechanics Publications. McGlynn, B.L., McDonnell, J.J., 2003. Quantifying the relative contributions of riparian and hillslope zones to catchment runoff. Water Resour. Res. 39, 1310. Mosé, R., Siegel, P., Ackerer, P., Chavent, G., 1994. Application of the mixed finite hybrid elements approximation in a groundwater flow model: luxury or necessity? Water Resour. Res. 30, 3001–3012.

Nearing, M.A., Jetten, V., Baffaut, C., Cerdan, O., Couturier, A., Hernandez, M., Le Bissonais, Y., Nichols, M.H., Nunes, J.P., Renschler, C.S., Souchre, V., van Oost, K., 2005. modelling response of soil erosion and runoff to changes in precipitation and cover. Catena 61, 131–154. Ogden, F.L., Watts, B.A., 2000. Saturated area formation on nonconvergent hillslope topography with shallow soils: a numerical investigation. Water Resour. Res. 36, 1795–1804. Panday, S., Huyakorn, P.S., 2004. A fully coupled physically based spatiallydistributed model for evaluating surface/subsurface flow. Adv. Water Resour. 27, 61–382. Paniconi, C., Putti, M., 1994. A comparison of Picard and Newton iteration in the numerical solution of multidimensional variably saturated flow problems. Water Resour. Res. 30, 3357–3374. Ribolzi, O., Andrieux, P., Valles, V., Bouzigues, R., Bariac, T., Voltz, M., 2000. Contribution of groundwater and overland flows to storm flow generation in a cultivated Mediterranean catchement. Quantification by natural chemical tracing. J. Hydrol. 233, 241–257. Richards, L.A., 1931. Capillary conduction of liquids through porous medium. Physics 1, 318–333. Rubin, J., 1969. Numerical analysis of ponded rainfall infiltration. In: Proceedings of the Wageningen Symposium, International Association of Scientific Hydrology. Rutqvist, J., Barr, D., Datta, R., Gens, A., Millard, A., Olivella, S., Tsang, C.-F., Tsang, Y., 2005. Coupled thermal–hydrological–mechanical analyses of the Yucca mountain drift scale test – comparison of field measurements to predictions of four numerical models. Int. J. Rock Mech. Min. Sci. 42 (5–6), 680–697. Singh, V., Bhallamudi, S.M., 1998. Conjunctive surface–subsurface modelling of overland flow. Adv. Water Resour. 21, 567–579. Sklash, M.G., Farvolden, R.N., 1979. The role of groundwater in storm runoff. J. Hydrol. 43, 45–65. Smith, R.E., Woolhiser, D.A., 1971. Overland flow on an infiltrating surface. Water Resour. Res. 7, 899–913. Streeter, V.L., Bedford, K.W., Wylie, E.W., 2003. Fluid Mechanics, ninth ed. Mc GrawHill Company. Therrien, R., Sudicky, E.A., 1996. Three-dimensional analysis of variably-saturated flow and solute transport in discretely-fractured porous media. J. Contam. Hydrol. 23, 1–44. Touma, J., Vauclin, M., 1986. Experimental and numerical analysis of two-phase infiltration in a partially saturated soil. Transport Porous Media 1, 27–55. VanderKwaak J.E., 1999. Numerical simulation of flow and chemical transport in surface–subsurface hydrologic systems. Ph.D. Thesis, Department of earth sciences, University of Waterloo, Ontario, Canada. VanderKwaak, J.E., Loague, K., 2001. Hydrologic-response simulations for the R-5 catchment with a comprehensive physics-based model. Water Resour. Res. 37, 999–1013. van Genuchten, M. Th., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. 44, 892–898. Wasantha Lal, A.M., 1998. Weighted implicit finite-volume model for overland flow. ASCE J. Hydraul. Eng. 124, 941–950. Weiler, M., McDonnell, J., 2004. Virtual experiments: a new approach for improving process conceptualisation in hillslope hydrology. J. Hydrol. 285, 3–18. Weill, S., 2007. Modélisation des échanges surface/subsurface à l’échelle de la parcelle par une approche darcéenne multidomaine, Thèse Ecole des Mines de Paris. Woolhiser, D.A., Smith, R.E., Giraldez, J.V., 1996. Effects of spatial variability of saturated hydraulic conductivity on Hortonian overland flow. Water Resour. Res. 32, 671–678. Younes, A., Fontaine, V., 2008. Efficiency of mixed hybrid finite element and multipoint flux approximation methods on quadrangular grids and highly anisotropic media. Int. J. Numer. Methods Eng.. doi:10.1002/nme.2327.