WEC-C: a distributed, deterministic catchment model — theory, formulation and testing

Environmental Modelling & Software 16 (2001) 583–599 www.elsevier.com/locate/envsoft

WEC-C: a distributed, deterministic catchment model — theory, formulation and testing J.T. Croton b

a,*

, D.A. Barry

1,b

a Water & Environmental Consultants, 28 Orpheus Street, Robertson 4109, Robertson, Australia School of Civil and Environmental Engineering and Contaminated Land Assessment and Remediation Research Centre, The University of Edinburgh, Edinburgh EH9 3JN, UK

Received 31 August 2000; received in revised form 12 April 2001; accepted 1 May 2001

Abstract WEC-C is a distributed, deterministic, catchment-scale water flow and solute transport model containing a number of innovations not present in comparable models. For example, it allows for the imposition of catchment topological changes resulting from land uses such as surface mining. In addition, it features preferential vertical flow in the vadose zone modelled using a dual continuum approach. Numerically, it differs from many existing models in that it employs a direct linkage between the vertical and lateral solvers of its split-solver scheme and, due to its use of explicit solvers, is stable regardless of the form of the soil moisture characteristic curves. The WEC-C model framework is a rectangular grid of uniform cell size in the lateral plane combined with a system of soil layers, of variable thickness, in the vertical direction. Within this structure, the governing equations for flow and transport are discretised and solved. All parameters are defined locally in each computational cell so that all available data on catchment variability can be incorporated directly into the model. This paper describes the formulation of WEC-C, which is based on operator splitting with first-order accurate solvers for both the vertical and lateral flow and transport models. WEC-C was subjected to four stringent, synthetic tests designed to evaluate its accuracy by comparison with available analytical and numerical solutions. These showed that there were scale issues associated with the model, and that induced numeric dispersion of solutes could be significant. Nonetheless, it is suggested that WEC-C is still useful as a distributed catchment model.  2001 Elsevier Science Ltd. All rights reserved. Keywords: Water flow; Solute transport; Catchment water balance; Distributed model; Hydrological processes; Dual-continuum profile; Split operator

Software availability Name of software: WEC-C Developer: James T. Croton of Water & Environmental Consultants, 28 Orpheus St, Robertson, 4109 Australia Comment: The copyright of the WEC-C model is held by James T. Croton. The model is freely available to all researchers for non-commercial use at no charge, as source code for all except the solver engines and as a compiled executable, * Corresponding author. Tel.: +61-7-3345-6756; fax: +61-7-33456756. E-mail addresses: [email protected] (J.T. Croton), [email protected] (D.A. Barry). 1 Tel.:+44-131-650-7204; fax: +44-131-650-7276

with the only limits being that it can not be copied, or used for any commercial purpose, without the author’s permission. The model is written in FORTRAN 77 style with use of 90 features, and is 21,000 lines of code including comments. It is well structured and easy to follow despite its length. A set of user notes and sample applications are currently available and a manual is in preparation. E-mail: [email protected]

Software availability Name of software: MAGIC Developer: Geoff Mauger, Water & Rivers Commission of Western Australia Comment: The MAGIC system, which is a GIS facility

1364-8152/01/$ - see front matter  2001 Elsevier Science Ltd. All rights reserved. PII: S 1 3 6 4 - 8 1 5 2 ( 0 1 ) 0 0 0 4 4 - 5

584

J.T. Croton, D.A. Barry / Environmental Modelling & Software 16 (2001) 583–599

Nomenclature General Units of measurement are given for each definition A B c c0 C d D D D0 Dd DL Dm DN DT En g h h0 Is k kr K K Ks Kss Kx Ky Kz L Ld Ls Lv n p q Q qr qs ra S Sc Se

⌬t 2S ⌬x2 ⌬t (Ti+1,j +Ti,j ) n 2 2S i,j ⌬x Solute concentration, M/L3 Concentration of influent solution, M/L3 ⌬t (Ti+1,j +Ti,j ) n 2 2S i,j ⌬y Depth of water in a cell, L ⌬t (Ti,j+1+Ti,j ) n 2 2S i,j ⌬y Soil–water diffusivity, L2/T Soil water diffusivity at the residual moisture content, L2/T Dispersion coefficient, L2/T Longitudinal dispersion coefficient, L2/T Molecular dispersion coefficient, L2/T Coefficient of numerical dispersion, L2/T Transverse dispersion coefficient, L2/T Net daily evaporation, L/T Magnitude of gravitational acceleration, L/T 2 Pressure head, L Hydraulic head at x=0, L Invert level of the stream, L Unit vector in the z direction Water phase relative permeability Saturated hydraulic conductivity tensor, L/T Hydraulic conductivity, L/T Saturated hydraulic conductivity, L/T Lateral saturated hydraulic conductivity of the upper layer, L/T Saturated hydraulic conductivity in the x direction, L/T Saturated hydraulic conductivity in the y direction, L/T Saturated hydraulic conductivity in the z direction, L/T layer thickness, L Domain length, L Stream length for an element, L Root length per unit volume of soil, L/L3 Time level in the numerical solution Soil water pressure offset from atmospheric pressure, M/T 2L Darcy flux, L/T Surface flux, L/T Potential water flux to the roots per layer, L/T Flux due to interflow, L/T Average root radius for a layer, L Storativity Slope term for lateral flow to a stream channel, L/L Normalised effective saturation (Ti,j +Ti−1,j )

n i,j

J.T. Croton, D.A. Barry / Environmental Modelling & Software 16 (2001) 583–599

SL Ss t T TA TB Tlin Tu Tx Ty v Ws x y z

585

Product of specific storativity and layer thickness Specific storativity, 1/L Time, T Transmissivity, L2/T User input variable in transpiration relation User input variable in transpiration relation User input variable, L/T Transpiration per unit leaf area, L/T Transmissivity in the x direction, L2/T Transmissivity in the y direction, L2/T Magnitude of advection velocity, L/T Stream width, L Horizontal longitudinal spatial co-ordinate, L Horizontal transverse spatial co-ordinate, L Vertical spatial co-ordinate, L

Greek and operators a ⌬· f fs n q ⌰ qi qr qs r n y ye yr ⵜ ·∗ ·ni,j

Parameter in van Genuchten characteristic curve, 1/L Discretisation of · Hydraulic head, L Head in the upper-most soil layer, L Parameter in van Genuchten characteristic curve Volumetric moisture content Normalised moisture content Initial moisture content Residual moisture content Volumetric moisture content at natural saturation Soil water density, M/L3 Parameter used to define soil type Soil water pressure head, p/rg, L Air entry head, L Potential driving flow to roots, L Del operator, 1/L Intermediate value of · Value of · at x=i⌬x, y=j⌬y, t=n⌬t

used to handle distributed input and output data, is available in a WEC-C compatible form for no cost. MAGIC is a combination of FORTRAN code and a Windows Application that allows graphical display, manipulation and management of map data. Documentation in the form of on-line help and examples are available for MAGIC and a user manual is planned. A full copy of MAGIC with source code, example files, etc. is 43 MB. WEC-C can be operated as a stand-alone program without MAGIC; ASCII files are used to handle all data, but data handling can be messy unless the model domain is small.

1. Introduction Catchment models, often called watershed models, can be classified according to the hydrological processes described, the space and time scales used and the method of solution (e.g., Fleming, 1975; Singh, 1995). Process descriptions are categorised as stochastic — based on statistical relations between two or more variables; parametric — empirical equations describing water movement and storages; and deterministic, where the fundamental flow equations are used to simulate actual processes. Spatially, models are normally divided into lumped and distributed types. Methods of solution can include analogue, analytical and numerical. In general, model complexity increases from stochastic through parametric to deterministic; from lumped to distributed;

586

J.T. Croton, D.A. Barry / Environmental Modelling & Software 16 (2001) 583–599

and from analytical to numerical. The Water & Environmental Consultants — Catchment (WEC-C) model belongs to the most complex category in that it is a distributed, deterministic model of numerical form, that simulates both water and solute movement within a catchment. Its form is especially useful where direct surface runoff is not the only streamflow generation process; when below ground processes, that is interflow, and groundwater discharge are significant contributors. Distributed, deterministic catchment models, and their practical application as hydrological tools, have a considerable history — examples include SHE (Abbott et al., 1986; Jain et al., 1992), TOPOG (O’Loughlin, 1986) and IHDM (Calver and Wood, 1995). By way of comparison, WEC-C incorporates the following innovations: 앫 It permits the imposition of spatially and temporally varying topological changes, such as excavation, to allow representation of surface mining; 앫 It employs a direct linkage between the vertical and lateral solvers of its split-solver numerical scheme; 앫 It is computationally efficient, allowing distributed catchment modelling on a high-end PC platform; 앫 It is based on a modular structure, allowing relatively easy changes to the various sub-models; 앫 Due to its use of explicit solvers, it is stable regardless of the form of the soil moisture characteristics used; 앫 The solvers are formulated to always conserve both water and solute so there are no mass balance errors; 앫 The form of the lateral solvers means there are no issues associated with partially saturated domains, or transition of any part of the domain between “wet” and “dry” states. 앫 Solver accuracy is at least first order for all components including the split-solver linkage; 앫 It directly models vertical preferential flow, through the use of a dual continuum soil profile; 앫 If required, the vertical preferential flow pathways can be directly linked to a lateral preferential flow system to create a three dimensional preferred pathway system; 앫 A simple, but powerful, GIS-like system handles its arrays of input and output data. 앫 The code is available to researchers at no cost and is well structured and easy to read. Distributed, deterministic catchment modelling has stimulated a considerable body of literature on the appropriateness of models, parameterisation issues, data needs, numerical issues, catchment representations, etc. For example, Beven (1996) highlights the limitations of distributed, deterministic modelling with regard to process conceptualisation, scaling, parameterisation, calibration data, model verification and availability of data for running the model. While the authors acknowledge these limitations, their discussion is beyond the scope of this

present article. Croton and Bari (this issue) give perspective on applying WEC-C to an actual catchment, and present an application methodology which may allow successful use of distributed models on data-rich catchments.

2. Model theory WEC-C employs a rectangular grid of uniform cell size in the lateral plane combined with a system of soil layers in the vertical, to represent the regolith of a catchment Fig. 1. A uniform grid was chosen for the lateral plane to simplify data input and output. For the same reason, the number of soil layers has been kept uniform across the model domain (this allows a layer to be connected laterally across all cells), although the thickness of any layer, or its elevation compared with the model datum, is flexible. This structure permits any soil layering and surface topography to be modelled. The top of the soil profile is defined by the soil surface while the bottom is defined by an impermeable basement surface (usually taken as the top of the parent rock if its permeability is low). The catchment is delineated by defining as active only those cells within the catchment divide; inactive cells are impermeable and act as solid boundaries to lateral flow within the model. All parameters are defined locally in each model cell so that all available data on catchment variability can be directly and fully incorporated into the model, though such flexibility can create its own problems if poorly used. The WEC-C model uses operator splitting such that fluxes are alternately computed for the vertical and lateral models; rather than a true three-dimensional solver which would solve for the whole domain simultaneously Fig. 2. In terms of computation, operator splitting has the advantage of being orders of magnitude faster than a full simultaneous solution; it is this feature that enables WEC-C to simulate real catchments on high end PCs. The disadvantage is that, in reality, all hydrological processes are occurring simultaneously and some miscalculation of fluxes must occur if they are calculated sequentially. There has yet to be a full theoretical evaluation of the splitting error for the flow and transport models used in WEC-C. However, we show below, by comparison with a “non-split” solution, that the splitting error is likely to be acceptable. In addition, since the splitting error is proportional to the discretisation used, in principle it can be reduced (using finer discretisation) at the expense of an increase in computational burden. The vertical and lateral models share the same spatial framework and the heads within the lateral models are defined directly from the soil water potentials of the vertical models. The pressure–water content relation for the two models is also common for positive pressures. This commonality is made possible through the simple

J.T. Croton, D.A. Barry / Environmental Modelling & Software 16 (2001) 583–599

Fig. 1.

587

Schematic of model layout and processes modelled.

explicit form of WEC-C’s solutions to the flow equations and allows a more direct linkage between the vertical and lateral components than is possible within models such as TOPOG which have a spatial separation between the vertical models and the lateral model for the permanent aquifer. Sub-surface lateral flows are simulated in each of the layers by an explicit finite-difference formulation that can have constant-head, constant gradient or zero-flux boundary elements in the bottom-most layer of the profile and bank seepage to stream channels in the top-most layer. Constant-head and constant gradient boundary cells permit simulation of groundwater fluxes across a catchment divide, particularly down-valley flow beneath a stream gauging station. All lateral flow uses the saturated hydraulic conductivity for a cell layer, combined with a depth of water based on the matrix potential. It is modelled as two-dimensional flow in an aquifer (Bear

and Verruijt, 1987) provided that the sum of the matrix potential and elevation head of a cell layer exceeds the elevation of the bottom boundary of the cell layer. Maximum depth of water for a layer is limited to the layer thickness. Vertical fluxes between layers within the model are handled by a dual continuum moisture model. This model includes, as boundary conditions, rainfall interception by a vegetation canopy and evaporation from the soil surface. Within the soil profile it includes flow to plant roots and plant transpiration as source/sink terms (Fig. 1). The solution is based on an explicit scheme applied to Richards’ equation (Richards, 1931) within each continuum. There are a number of significant differences between this explicit solver, and the implicit schemes normally associated with soil water flow models. These will be discussed in detail below. The dual continua are conceptualised as a system of preferred pathways, cylindrical in form, containing highly

588

J.T. Croton, D.A. Barry / Environmental Modelling & Software 16 (2001) 583–599

∂q q ∂y ⫹ S ⫽ⵜ·[Kkr(ⵜy⫺k)], ∂t qs s ∂t

(1)

where symbol meanings are defined in the notation list. The assumptions that underlie Eq. (1) are: 앫 앫 앫 앫 앫 앫 앫

elastic porous medium and liquid; z axis points downwards; isothermal conditions; negligible pressure gradients in the gas phase; no sources or sinks of water; qs is constant; relative permeability is homogenous and isotropic; and 앫 soil water is homogeneous.

Fig. 2.

Flow chart for the WEC-C program.

permeable material, within a matrix of lesser permeability. The linkage between the matrix and preferredpathway profiles is based on radial flow under both saturated and unsaturated conditions. WEC-C considers overland flow via a nine-point finite difference solver that considers both a Manning’s equation-like algorithm and the kinematic wave equation. Channel flow can be either by simple reach summation; or by a channel routing algorithm based on both Manning’s equation and the kinematic wave equation, which is computationally intensive due to the very short time steps required to maintain stability within the channel routing algorithm. The unit of time within the model is a day, and input data for evaporation and the salinity of rainfall are input at this time-step. However, rainfall data can be input at any fraction of a day although, once set, the fraction must remain the same. There are both user-defined controls and automatic optimisers incorporated within WEC-C, to control the computational time-step length so that numerical accuracy and stability are preserved regardless of the data input time-step.

3. Model formulation 3.1. Single-phase flow The major portion of the WEC-C model is simulation of water flow through the soil profile. The governing equation for single-phase liquid flow in porous media is (e.g., Pinder and Gray, 1977):

If a source adds water to the soil profile, then a term quantifying its behaviour is added to the right-hand side of Eq. (1). For simplicity, such terms will be ignored below, although they are included, as appropriate, in the model. Eq. (1) is written in a Cartesian (x,y,z) co-ordinate system. If the principal directions of K are assumed to be aligned with these axes, then Eq. (1) becomes:

冋

册 冋

∂ ∂y ∂y ∂q q ∂y ∂ ⫹ S ⫽ Kk ⫹ Kk ∂t qs s ∂t ∂x x r ∂x ∂y y r ∂y

冋 冉 冊册

册

(2)

∂y ∂ ⫺1 , ⫹ Kzkr ∂z ∂z where 0ⱕkrⱕ1. The solution to Eq. (2) requires specification of the boundary and initial conditions. The saturated hydraulic conductivities, Kx, Ky and Kz, may vary with location, while kr can be specified as a function of either q or y. In the saturated zone (y⬎0, ∂q/∂t=0), Eq. (2) must be solved using y as the independent variable. In the unsaturated zone (y⬍0), either y or q may be used. One converts between each formulation using the soil moisture characteristic curve which relates y and q. 3.2. Split-operator solution of the saturated– unsaturated flow model WEC-C implements a split-operator solution to Eq. (2); below we describe the approximations made in arriving at the WEC-C model formulation. Eq. (2) is split into vertical and lateral steps. Assume that at some time step n (where t=n⌬t), q(x,y,z,n⌬t) and y(x,y,z,n⌬t) are known. The split-operator solution over a time step ⌬t proceeds by solving the following equations: Step 1 — Solution for vertical flow using q(x,y,z,n⌬t) and y(x,y,z,n⌬t) as initial conditions

冋 冉

冊册

∂y∗ ∂q q∗ ∂y∗ ∂ ⫹ Ss ⫽ Kzkr ⫺1 , ∂t qs ∂t ∂z ∂z

(3)

J.T. Croton, D.A. Barry / Environmental Modelling & Software 16 (2001) 583–599

where q∗ and y∗ are intermediate results, and are used as the initial condition in the next step. Step 2 — Solution for lateral flow using q∗(x,y,z,n⌬t) and y∗(x,y,z,n⌬t) as initial conditions

∂q∗ ∂q∗ ⫽ . ∂t ∂z

∂ ∂y ∂y ∂q q ∂y ∂ ⫹ Kykr . ⫹ Ss ⫽ Kxkr ∂t qs ∂t ∂x ∂x ∂y ∂y

q∗(t+⌬t)−q∗(t) 1 ∂ ∗ ∗ ⫽ [q (q (t))⫹q∗(q∗(t⫹⌬t))], ⌬t 2 ∂z

冋

册 冋

册

(4)

The solution to Eq. (4) gives q and y at t=(n+1)⌬t. Whatever the accuracy of numerical schemes used to solve Eqs. (3) and (4), the sequential scheme employed here introduces a splitting error. Different splittings were considered, and that of Strang (1968) which consists of carrying out Step 1 over ⌬t/2, Step 2 over ⌬t, and Step 1 again over ⌬t/2 was adopted as a compromise between simplicity and computational efficiency, and accuracy. For linear partial differential operators, Strang splitting maintains second-order accuracy; but not for non-linear operators such as those here. Miller and Rabideau (1993) and Barry et al. (1996a) have made these observations in the context of temporal, as opposed to spatial, splitting. Thus, it is reasonable to conclude that the WEC-C split-operator approach is first-order accurate in the spatial discretisation (due to the scheme used for the lateral flow model it is also first-order accurate in the time step, as discussed below). 3.3. The WEC-C solution of Eq. (3) For each vertical “stack” of cells Fig. 1, Eq. (3) is solved by WEC-C, except that WEC-C takes Ss=0 for the unsaturated case. This simplification is standard when considering unsaturated flow, and amounts to the assumption of a rigid soil-matrix through which an incompressible liquid (water) is flowing. In the transition to saturated flow, however, when the soil moisture capacity approaches zero, the second term on the lefthand side of Eq. (3) becomes relatively more important. WEC-C does not account for this in the vertical flow model, although it is included in the lateral flow. Note that because vertical flow takes place in both the saturated and unsaturated zones, vertical saturated flow in WEC-C is that of an incompressible liquid through a rigid soil-matrix. These assumptions are usually considered acceptable approximations in groundwater hydrology. WEC-C assumes a moisture-content based formulation of unsaturated flow in the derivation of its numerical solution. The usual limitation of this approach is that soil–water pressure heads greater than atmospheric (saturated soil) cannot be modelled since, by definition, the moisture content is at saturation and cannot vary. However, in WEC-C the problem of saturated cells does not arise, because in the implementation of the numerical scheme the moisture content is converted to the equivalent depth of water within each cell. Eq. (3) can be written in the form (with Ss=0):

589

(5)

In WEC-C, this equation is discretised as: (6)

which is second-order accurate in ⌬t. On the right-hand side, a first-order expression for the moisture content at the advanced time level is used: q∗(t⫹⌬t)⫽q∗(t)⫹⌬t

∂q∗(q∗(t)) . ∂z

(7)

Use of Eq. (7) renders the scheme in Eq. (6) first-order accurate in ⌬t. The spatial derivatives in Eqs. (6) and (7) are discretised as: ∗1 ∗1 ∂q∗ qi+2−qi−2 ⫽ , ∂z ⌬z

(8)

where the subscript on the right-hand side indicates the node at which the term is evaluated (an average of adjacent nodes is calculated). Eq. (8) is second-order accurate in ⌬z. Thus, overall the numerical scheme for the vertical flow calculations has an error term of O(⌬z2,⌬t). Other solution approaches are possible. As already noted, WEC-C uses Richards’ equation with Ss=0, giving a form that has been subjected to many numerical treatments:

冋

册

∂q∗ ∂q∗ ∂ ⫽ D(q∗) ⫺Kzkr(q∗) , ∂t ∂z ∂z

(9)

where D is the soil–water diffusivity. Eq. (9) can be solved analytically in some cases (e.g., Barry and Sander, 1991; Barry et al., 1993), and can be used to test the WEC-C solution (see below). 3.3.1. Dual continuum soil profile Generally, for soils that contain significant differences in pore sizes, the assumption of a single domain medium yields poor flow predictions (e.g. Grayson et al., 1992). An alternative is to model the flow domain as a pair of overlapping continua. Various approaches based on this idea have appeared (e.g., Duguid and Lee, 1977; Edwards et al., 1979; Davidson, 1985; Pruess et al., 1990; Gerke and van Genuchten, 1993). WEC-C adopts a dual continuum approach based on vertical cylinders of high conductivity within a matrix of lower conductivity. The relationship between the two flow domains needs to be specified, as well as the flow characteristics within each region. In WEC-C, flow occurs simultaneously within each region; along with inter-region exchange, the extent of which depends linearly on the difference in soil–water pressure heads between regions. The governing equations for flow are simply Eq. (3)

590

J.T. Croton, D.A. Barry / Environmental Modelling & Software 16 (2001) 583–599

written for each region, with additional terms describing the inter-region exchange of water. The direct solution method involves solving the coupled equations simultaneously. To prevent rapid oscillations in the interregion flux, a user-defined allowable maximum flux is invoked which, if necessary, provides a criterion by which the time-step is broken into finer steps. The exchange between continua is based on the radial flow equation (Terzaghi and Peck, 1967), which is essentially Darcy’s law factored by a transfer-efficiency coefficient. Gerke and van Genuchten (1993) use a mathematically similar transfer law. 3.3.2. Solute transport modelling Advection is the major solute transport mechanism in WEC-C: physically, solute (in the case of the catchment studies to date, salinity) moves with the water in the soil, in plug-flow fashion. For a given continuum within a computational cell of the model, all inputs are mixed before advection at the next time step. As shown below, this mixing introduces solute dispersion into the model. Between layers and continua, a solute diffusive flux is also included, though its effect is usually small at the scale of a catchment. 3.4. The WEC-C solution of Eq. (4) In WEC-C, lateral flow is treated as being depth-integrated for each layer in the discretised domain. Fluxes from one layer to another usually appear as source terms in the governing equation for each layer, but the vertical flow model, Eq. (3), handles inter-layer fluxes. In arriving at the appropriate governing equation for each layer a number of assumptions are made:

SL⫽SsL,

(11)

Tx⫽KxkrL,

(12)

and Ty⫽KykrL,

where L is the layer thickness. Since the functional relationship q and y is known, the left-hand side of Eq. (10) can be rewritten with f as the sole dependent variable: S

冋 册 冋 册

∂ ∂f ∂f ∂ ∂f Tx ⫹ Ty , ⫽ ∂t ∂x ∂x ∂y ∂y

where S⫽L

dq q ⫹ S . dy qs L

3.4.1. Explicit finite difference solution of Eq. (14) Forward differencing of the temporal derivative in Eq. (14) is carried out, while the spatial derivatives are evaluated by central differences at the old (t=n⌬t) time level. Thus, the scheme is first-order accurate in time (Noye, 1982). The solution at the new time level t, i.e., (n+1)⌬t, can then be calculated. In the split-operator solution, the “old” time level means the result from Step 1, the vertical flow step, i.e., y∗, gives the initial condition: fnij ⫽y∗(n+1) ⫹zi,j . i,j

(16)

(10)

⫹

冢

⌬t ∂f ⫺f )⫽ Ti+1,j 2 ⌬x ∂x

n+1 i,j

S (f WEC-C treats vadose zone lateral flow in an identical manner to saturated zone lateral flow. For the unsaturated zone, this amounts to assuming that the saturation in each layer is uniform vertically. To arrive at the WECC formulation, we first use the definition of hydraulic head, f=y⫺z in Eq. (4); the resulting equation is identical except that y is replaced by f. Next, we integrate over an arbitrary layer, assuming that all variables are uniform with z (but not x and y) over the depth of the layer. The resulting equation is (assume for simplicity that there is zero water flux through the top and bottom of the layer):

where

(15)

In the WEC-C model, the second term on the right-hand side of Eq. (15) is ignored; this is consistent with the model assuming Ss=0, as mentioned above. The discretised form of Eq. (14) is derived below.

n i,j

冋 册 冋 册

(14)

Eq. (14) is discretised in two stages. The first stage gives:

앫 Kx and Ky are uniform with depth; 앫 Ss is uniform with depth; 앫 Equipotential lines are vertical.

∂ ∂f ∂f ∂q q ∂y ∂ ⫹ T , L ⫹ SL ⫽ Tx ∂t qs ∂t ∂x ∂x ∂y y ∂y

(13)

冢

n i,j

∂f ⌬t T 1 ⌬y i,j+2∂y

|

n

|

⫺T 1 i+ ,j 2

1 i− ,j 2

n

| 冣

∂f ∂x

(17)

1 i− ,j 2

n

| 冣

∂f 2 ∂y

⫺Ti,j⫺1 1 i,j+ 2

n

,

1 i,j− 2

where dq∗ q∗i,j S ni,j ⫽L ∗兩ni,j ⫹ SL. dy qs

(18)

Next, the right-hand side of Eq. (17) is discretised again to move the function evaluations to the calculation nodes, i.e., ⌬t n−1 [(T ⫹T )(fn ⫺fn ) S ni,j (fn+1 i,j ⫺fi,j )⫽ 2⌬x2 i+1,j i,j i+1,j i,j

J.T. Croton, D.A. Barry / Environmental Modelling & Software 16 (2001) 583–599

⫺(Ti,j ⫹Ti−1,j )(fni,j ⫺fni−1,j )]⫹

⌬t [(T 2⌬y2 i,j+1

(19)

⫹Ti,j )(fni,j+1⫺fni,j )⫺(Ti,j ⫹Ti,j−1)(fni,j ⫺fni,j−1)], Eq. (19) can be rewritten as: n n n n n n fn+1 i,j ⫽fi,j ⫹A(fi,j ⫺fi−1,j )⫹B(fi+1,j ⫺fi,j )⫹C(fi,j

(20)

⫺fn1,j−1)⫹D(fni,j+1⫺fni,j ). Eq. (20) is used in the WEC-C model.

591

3.5.2. Transpiration Transpiration is assumed to be controlled by supply, flow to roots, and by available energy to drive the atmospheric-limited process. Extraction of water from the soil by the root system of the vegetation cover is based on the radial flow equation (Whisler et al., 1970). The equation applied within each computational cell is: q r⫽

2pLvL

冋 冉 冊册

ln 1/ ra冑pLv

K(y)⌬yr.

(21)

3.5. Other catchment processes (source/sink terms) The preceding discussion of the WEC-C formulation was made without reference to water source and sink terms. Source and sink terms are the interface between the numeric schemes of the model and the other hydrological processes of the catchment. For each cell of the model domain WEC-C considers, as appropriate: 앫 canopy interception and evaporation, 앫 plant transpiration, 앫 soil surface infiltration and discharge of soil water to the surface, 앫 soil evaporation, 앫 bank seepage to stream channels, 앫 direct surface runoff, and 앫 groundwater discharge from the catchment. Within the model, all evaporation and transpiration fluxes are driven by a potential evaporation term that is presently based on daily evaporation data from a Class A pan. Other methodologies, such as those based on the Penman–Monteith equation, are more conceptually realistic but are applicable only with the requisite meteorological data; these data are lacking in the catchments where the WEC-C model has so far been applied. The evaporation and transpiration functions outlined below have been developed by us and must be considered preliminary. It is intended to outline their development, and comparison with measured data, in future articles. 3.5.1. Interception and canopy evaporation A simple, single store model with rainfall as input and canopy evaporation and throughfall as outputs is used. To account for the decrease in evaporation rate per unit leaf-area as the canopy closes, a maximum effective leaf-area is defined: leaf-areas above this figure, though increasing canopy storage, do not increase total canopy evaporation. By necessity, realistic interception rates are achieved by an interception efficiency term, at present a simple multiplier.

The potential fluxes for each layer are calculated by Eq. (21) and summed to obtain the total root-flow-limited transpiration for the profile. Atmospheric-limited transpiration is calculated as a linear function of net evaporation when the latter is less than an input variable (Tlin), and as a logarithmic function when it is greater. The equations for transpiration per unit leaf-area are given below. The form of these equations has been based on plant water use measurements made by the heat pulse technique (Marshall and Chester, 1992; Marshall, 1993).

冦

TA ln (En)+TB,

En⬎Tlin

Tu⫽ En[TA ln (Tlin)+TB] , otherwise. Tlin

(22)

The lesser of the root-flow-limited and atmosphericlimited transpiration rates is taken to be the transpiration rate. When transpiration is atmospheric-limited, water is first extracted from the layer and continua with the lowest potential referenced to the soil surface, then extracted from the second lowest, and so on until the demand is filled. 3.5.3. Surface infiltration and discharge Surface infiltration and discharge are calculated using Eq. (6) but assuming that ⌬z is equal to half the first layer thickness, that the computation node for surface water is at the soil surface. 3.5.4. Soil evaporation Soil evaporation is assumed to occur from the top layer of the profile, and is modelled via Eq. (6). It is driven by the difference in potentials between the top layer and a user-defined air potential. The total flux is limited to a fraction of potential evaporation, the fraction being a function of canopy cover. Soil evaporation is also reduced by a leaf litter term.

592

J.T. Croton, D.A. Barry / Environmental Modelling & Software 16 (2001) 583–599

3.5.5. Bank seepage to streamlines The lateral model, Eq. (20) allows bank seepage to a streamline within the upper-most soil layer. The assumed form is that of saturated Darcian flow from the soil matrix to a stream channel via the channel’s sidewall. The equation used is: qs⫽2(fs⫺Is)LsScKss.

(23)

The term Sc is the average surface slope between the stream cell and its surrounding upslope cells and the saturated thickness is considered by the (fs⫺Is) term. 3.5.6. Overland flow Overland flow is calculated as sheet flow by an explicit nine point predictor–corrector solver that takes the individual minima of Mannings equation and the kinematic wave equation. Flow gradients between cells are calculated using the soil surface as reference, as this greatly reduces the computational burden; and flow widths are defined as half the cell width when flow is along a principle axis, and as half the cell diagonal for diagonal flows. 3.5.7. Stream routing Stream routing is an option within the WEC-C model and is based on the minima of Mannings equation and the kinematic wave equation in an explicit predictor– corrector solver; this can be computationally intensive if the channel is smooth and/or steep, as the channel is segmented on the same cell structure as the main model. For most small catchments, it has been found that little in terms of realism is sacrificed if channel flow time (short for small catchments) is neglected and all water deemed to discharge on entering the stream. There are also options in which channel flows can be collected by reach, allowing simulation of special cases such as flows along road drains to sediment ponds. 3.5.8. Groundwater discharge from the catchment The lateral model of the deepest layer can have boundary cells that are either no-flow, constant-head or constant-gradient. This means down-valley groundwater flows, or any other subsurface flows across the catchment boundary, can be simulated by defining constanthead or constant-gradient boundary cells there. 3.6. Options in the vertical solver of WEC-C The solution of Eq. (3), as discussed in Section 3.3, requires an estimate of hydraulic conductivity (K) for interlayer fluxing. Obtaining such a value is not simple as water contents, and their related hydraulic conductivities, are known only at the layer centroids and not between. Haverkamp and Vauclin (1979) tested various approximations including those based on the arithmetic,

harmonic and geometric means of the centroid conductivities, as well as an “upstream mobility” concept using the conductivity of the layer contributing water. They also considered schemes involving linear interpolation and extrapolation. Haverkamp and Vauclin concluded the geometric mean provided the least weighting error. The authors have found none of the methods suitable, each having its particular problems. For instance, use of the geometric mean will result in artificially low fluxes if the flux is from a moist clay to a dry sand: the low unsaturated hydraulic conductivity of the dry sand dominates the calculation. Conversely the arithmetic mean will over supply water if the flow in the example is reversed, that is from a saturated sand to a dry clay: the high saturated hydraulic conductivity of the sand dominates the calculation. While the WEC-C model contains an option that allows calculation of conductivities using the arithmetic mean, the authors are of the opinion that the correct solution is the following logic scheme which considers the hydraulic conductivities and matrix potentials of the layers and defines the actual limits to fluxing. It has been implemented as an option in WEC-C. The importance of using such a scheme is outlined in Test 3 below. The ability to incorporate such a scheme is one of the true advantages of the WEC-C explicit solvers over the implicit ones used in most other models. 앫 When the hydraulic conductivity of the receiving layer is greater than that of the contributor, the conductivity of the contributor is used. 앫 When the flux is downwards and the upper layer has a positive head when referenced to its lower boundary, that is y⬎⫺dz/2, then the saturated conductivity of the lower layer is used; unless the conductivity of the upper layer is less, in which case that is used. 앫 When the flux is upwards and the lower layer has a positive head when referenced to its upper boundary, that is y⬎dz/2, then the saturated conductivity of the upper layer is used; unless the conductivity of the lower layer is less, in which case that is used. 앫 If none of the above cases apply, then the conductivity of the receiving layer is used. Another important concept introduced into WEC-C as a solver option is that of air-entry limiting the effectiveness of a given value of y to drive fluxes between continua and layers. It is well known that soil is far from uniform, and that textural discontinuities exist at all length scales. These discontinuities break the capillary column and limit the ability of large negative values of y to draw soil water over large distances. Air-entry is specified in the WEC-C model as a y value for a given soil. Once this value is exceeded, the model uses effective y values to drive the fluxes. These are presently based on the square root of the actual y

J.T. Croton, D.A. Barry / Environmental Modelling & Software 16 (2001) 583–599

values as this functional form creates a smooth transition. Air-entry requires further research, but for now it is an important concept that is included in the WEC-C model and greatly adds to its realism.

4. Model testing Here, we examine a series of WEC-C model test cases to check important components in various operational modes. Although more thorough testing could be undertaken, the tests show how and where the WEC-C solution scheme leads to results that differ from those of a more “exact” analysis. The results of the tests are discussed in the context of the catchment-scale applications for which the WEC-C model was designed. The tests performed, and the reason for their choice, are: 앫 Test 1: Steady Horizontal Flow (Confined or Unconfined Aquifer). This test is a check of the lateral groundwater flow model. For steady, one-dimensional flow, simple analytical solutions are available for both confined and unconfined flow. The analytical results were compared with the WEC-C numerical predictions. 앫 Test 2: Steady Water Flow with Transient Solute Transport. As mentioned above, at the end of each time-step, the solute transport model in WEC-C induces dispersion by mixing water within cells. For steady, uniform water flow (i.e., no flow variations with time or in space), the solute transport can be modelled using an exact analytical solution. The case considered is analogous to steady flow in a laboratory soil column into which solute is introduced at a particular instant. Solute concentrations in the column effluent are then monitored to measure the solute breakthrough curve. This situation can be modelled in WEC-C by considering flow in a uniform layer, with injection at one end and with solute concentrations measured at some distance along the layer. 앫 Test 3: Vertical Infiltration. Tests 1 and 2 are concerned with saturated confined and unconfined flow, and deal with the lateral flow model exclusively. Test 3 is concerned with vertical flow through the unsaturated zone, which is an important component of any catchment model; this part of WEC-C is examined explicitly by comparison with an exact analytical solution. 앫 Test 4: Three-Dimensional Flow and Transport. Test 4 is fully three-dimensional, and involves constantflux infiltration over a surface strip, with pumping in another part of the aquifer. Both flow and solute transport are simulated for saturated and unsaturated conditions.

593

4.1. Test 1: steady horizontal flow (confined or unconfined aquifer) A simple initial test of the WEC-C model is that of steady lateral flow in a confined, one-dimensional, uniform aquifer. Darcy’s Law gives the steady water flux. For an unconfined aquifer, the steady lateral flux of water in one dimension is given by the Dupuit– Forchheimer discharge formula (e.g., Bear and Verruijt, 1987): K q⫽⫺ [h20⫺h2(x)]. 2x

(24)

A number of cases were run with different parameter values in Eq. (24). Typical results are shown in Fig. 3. The abscissa measures the WEC-C model predictions for h2(x). For the example considered, h0=1, as shown in the figure. Rearrangement of Eq. (24) gives the slope of h2(x) as 2xq/K. Fig. 3 is a plot of the exact solution from Eq. (24) and the WEC-C output. A linear leastsquares fit to the WEC-C output, also shown in Fig. 3, gives the slope as 0.06 (units are unimportant), which is exact for the parameter values used. The intercept is 1.0002, which is very close to the exact value of unity. Thus, the WEC-C results are very accurate for this test. 4.2. Test 2: steady water flow with transient solute transport Steady flow in a one-dimensional, uniform porous medium is described simply by Darcy’s Law. Consider, for a moment, the flow situations described in Test 1. For flow in a confined aquifer, we consider the case where the cross-sectional area is uniform, the porous medium is homogeneous, and the head drop imposed at the boundaries is also uniform (i.e., the rate of flow entering the medium is uniform in cross-section). Under steady conditions, the flow rate does not vary across any cross-section within the flow domain, and the problem can be considered one-dimensional.

Fig. 3. WEC-C results for Test 1. The line is a linear regression of the WEC-C output.

594

J.T. Croton, D.A. Barry / Environmental Modelling & Software 16 (2001) 583–599

For steady, unconfined, horizontal flow, the depth of flow varies with position: h is a function of x in Eq. (24). Since total flow through any cross-sectional area is constant, the flow rate per unit depth varies. In other words, unless a large head drop is imposed, the flow predicted for this unconfined aquifer satisfies well the Boussinesq assumption of hydrostatic conditions along any vertical axis Dagan, 1967; Parlange et al., 1984; Barry et al., 1996b). Thus, flow rates at a given point within the aquifer do not vary with vertical position, although they can with lateral position. In Test 2, both cases (confined and unconfined) considered in Test 1 were examined, though only a single case is presented here: that of flow in a confined aquifer, with a solute source added at one boundary at time zero. The unconfined aquifer case gave very similar results. WEC-C, by construction, introduces numerical dispersion because of the mixing-cell approach. That is, WECC models solute transport by advection only, but between the advective steps it mixes the water in the cells and this mixing induces the numerical dispersion. For steady flow conditions, the numerical dispersion induced, DN, is given by (Bajracharya and Barry, 1994): v DN⫽ (⌬x⫺v⌬t). 2

(25)

For confined flow, v is calculated as v=q/qs where qs is the porosity of the porous medium. Under steady flow conditions, the analytical solution for solute movement with a first-type entrance condition and assuming an initially solute-free medium is (Lapidus and Amundson, 1952):

Fig. 4. Comparison of WEC-C results with Eq. (26) for solute transport under steadyflow conditions for Test 2.

4.3. Test 3: Vertical infiltration The operator splitting in WEC-C separates the vertical and lateral flow portions of water movement. This test exclusively considers vertical flow in the unsaturated domain. It is designed as a “worst-case” scenario, i.e., infiltration at a high rate into an initially dry soil. Sander et al. (1988a) (cf. Broadbridge and White, 1988a,b; Sander et al., 1988b; White and Broadbridge, 1988) have published an exact solution for constant-flux infiltration into a dry soil. That is, they provide a solution for Richards’ equation (this same equation is solved in the WEC-C model):

冋

册

∂q dK(q) ∂q ∂q ∂ ⫽ D(q) , ⫺ ∂t ∂z ∂z dq ∂z

(27)

with D(q)⫽K(q)

dy , dq

(28)

subject to:

冤冢 冣

冉冊

冢 冣冥

x−vt vx x+vt c(x,t) 1 ⫽ erfc ⫹exp erfc c0 2 D d 冑4Ddt 冑4Ddt

.

q(z,0)⫽0, (26)

The WEC-C model can be checked by replacing Dd in Eq. (26) with DN from Eq. (25). The analytical solution was compared with the WECC results. Again, good agreement was obtained, as shown in Fig. 4, though it should be noted that such a good agreement may not be obtained for another case. An implication of the induced numerical dispersion is that solute transport will be grid-size dependent. In most practical applications of WEC-C, the product v⌬t will be small compared with ⌬x, so the numerical dispersion in the direction of flow will be approximately v⌬x/2. Transverse to the flow, there will be no diffusion/dispersion of solute.

zⱖ0,

q(z,t)⫽0, z→⬁, t⬎0,

(29) (30)

and K(q)⫺D(q)

∂q ⫽Q, z⫽0, t⬎0. ∂z

(31)

Eq. (27) is just another form of Eq. (9). The normalised volumetric moisture content, ⌰, is related to the actual moisture content, q, via: ⌰⫽

q−qi . qs−qi

(32)

The hydraulic functions, D and K, are defined by: D0 , D(⌰)⫽ (1−n⌰)2 and

(33)

J.T. Croton, D.A. Barry / Environmental Modelling & Software 16 (2001) 583–599

Table 1 Parameter values used in the simulations for Test 3 Parameter

Value

ye Q Ks D0 n qi qs

⫺0.2 m 0.04 m/d 0.05 m/d 0.05 m2/d 0.85 0.1 0.55

Ks(1−n)⌰ K(⌰)⫽ , 1−n⌰

(34)

respectively. Observe that, as the soil type parameter n→1, the soil properties become increasingly non-linear. The soil moisture characteristic curve for these properties is: y⫽ye⫹

冋

册

D0 (1−n)⌰ ln . (1−n)Ks 1−n⌰

(35)

Note that K can also be expressed in terms of y as: K(y)⫽Ksexp

冋

册

Ks(1−n)(y−ye) . D0

(36)

Eq. (36) shows the typical exponential dependence on soil water pressure head, y, of the unsaturated hydraulic conductivity, K. That is, the soil hydraulic properties described by Eqs. (33) and (34) are, sometimes, reasonable approximations of real soil properties. The hydraulic properties for the model described above were implemented in the WEC-C model for the case shown in Table 1. In the Sander et al. (1988a) model, the definition of Q needs to be adjusted slightly to take into account the conversion from normalised to actual moisture contents. That is, the Q appearing in Eq. (9) of Sander et al. needs to be factored by (qs⫺qi)−1. The results for the exact and WEC-C solutions for simulations with 10 layers (uniform layers of 2 m thickness) are shown in Fig. 5. There are four WEC-C

Fig. 5. Test 3 moisture content profiles for WEC-C and Sander et al. (1988a).

595

solutions, three for estimates of the between-layer hydraulic conductivity based on the arithmetic, geometric and harmonic means and one for the case-based logic scheme particular to WEC-C (see Section 3.6 for an explanation of this scheme). All three of the solutions based on mean values fit poorly to the exact solution by Sander et al.. In particular, the wetting fronts for the geometric and harmonic solutions have entered only the first layer, and surface water totalling 4280 mm depth is present at the end of these simulations. The reason for the lack of penetration of the wetting front is that the low initial hydraulic conductivity of the soil, only 7.5×10−9 mm/day, is being used in the calculation of the between-layer conductivity at the wetting front, effectively zeroing the calculation and preventing the movement of the wetting front from one layer to the next. For the arithmetic mean, the wetting front is at least close to correct position, but the moisture contents behind it are excessive. The case-based logic scheme of WEC-C produces a much closer match to the exact solution by Sander et al.. Given the coarse spacing used, layers of 2 m thickness, the match is exceptional. This level of correspondence to the exact solution should more than suffice for catchment modelling given: the natural variability in unsaturated flow parameters (e.g., Sharma et al., 1987); the known lack of correspondence between actual infiltration processes and Richards’ equation at the field scale (e.g., Beven, 1996); and that a key component of WEC-C is its dual-continuum approach which is not being tested here. 4.4. Test 4: three-dimensional flow and transport Test 4 is fully three-dimensional, and involves constant-flux infiltration over a surface strip, with pumping of the basement aquifer. Both flow and solute transport are simulated. This test assesses the general performance of WEC-C, particularly the split-operator scheme’s ability to mimic three-dimensional flow. This is the Test 4 set-up: water is infiltrated, at a constant-flux of 0.1 m/d, at a 10 m by 40 m disposal site and passes through the unsaturated zone to an unconfined aquifer. Down gradient, water is pumped from the aquifer via a well (Fig. 6). Initially, the aquifer is assumed to be in vertical equilibrium, i.e., the equipotentials are vertical and the flow is horizontal only. Across the flow domain the phreatic surface (atmospheric pressure) location varies linearly from z=28 m to z=26 m. The initial flow configuration is uniform in the transverse (y⫺z) plane. The screened section of the well is located from the impermeable base of the aquifer, rising to a distance of 20 m. The well is modelled as a fixedhead boundary condition with the head for a given point being equal to the elevation above the base. Solute enters the medium with the infiltrating water at the waste disposal site and has a concentration of 1 (units are

596

J.T. Croton, D.A. Barry / Environmental Modelling & Software 16 (2001) 583–599

Fig. 6.

Details of the flow domain and boundary conditions for Test 4.

irrelevant). Due to symmetry about y=0, only half the flow domain was actually modelled. As a comparison with WEC-C, the three-dimensional numerical solution to the flow and transport equations as given by Simunek et al. (1995), SWMSF3D, was used. The problem defined above was derived from an example given by Simunek et al. (1995). The numerical solution consisted of a grid composed of rectangular prisms with 10,560 elements and 12,144 nodes. Nodal spacings were adjusted in regions where rapid changes of moisture content were likely. A single soil type was used across the whole model domain, with its properties based on the van Genuchten model; the van Genuchtentype soil moisture characteristic curve and hydraulic function are given in Eqs. (37) and (38). Relevant parameters, including those for solute transport, are given in Table 2. Table 2 Parameter values used in the simulations for Test 4 Parameter

Value

qr qs a (1/m) Ks (m/d) n Dm (m2/d) DL (m2/d) DT (m2/d)

0.05 0.45 4.1 5 2 0.01 0.21 0.21

冦

qs−qr , y⬍0 (1+兩ay兩n)1−1/n , q(y)⫽ yⱖ0 qs, qr+

and K(y)⫽

再

Ks冑Se[1−(1−S 1−1/n ]2, y⬍0 e . yⱖ0 Ks

(37)

(38)

The usual transport equation for a tracer solute is solved in the SWMSF3D model; for a discussion of this equation and its derivation see, for example, Bear and Verruijt (1987). As can be inferred from Table 2, the molecular diffusion coefficient (Dm) is too low to have any effect on the results. A single value for both the longitudinal (DL) and transverse (DT) dispersion coefficients was selected to roughly match the numerical dispersion of the WEC-C model, as predicted above by Eq. (25). It is impossible to match the WEC-C results precisely when transient flow on a non-uniform grid is used, as it is here. The WEC-C model was based on a 5 m by 5 m horizontal array with 10 vertical layers. Layer thickness was 2 m for the first 9 layers and 20 m for the bottom layer. There are 5200 active elements in the WEC-C model, which is about half that for the SWMSF3D. The WEC-

J.T. Croton, D.A. Barry / Environmental Modelling & Software 16 (2001) 583–599

C model was deliberately kept relatively coarse so as to be a realistic representation of the type of WEC-C model set-up normally developed for such a case. The pumping well in the WEC-C model is simulated by a constanthead cell in the bottom-most layer. Fig. 7(a)–(c) shows a comparison of the results as vertical slices at y=0. The times shown are 20, 100 and 500 days. The pressure heads are generally similar for both models though there are marked differences about the infiltration area and near the well, this latter being most obvious at 20 days (Fig. 7(a)). The differences in head about the well are probably primarily related to how the well is represented: in SWMSF3D it is a line of constanthead; while in WEC-C it is a constant-head cell 5 m square, which allows a greater flow rate. The differences in head about the infiltration site result from lateral flow in WEC-C being dependent on the generation of saturated conditions while in SWMSF3D it is true unsaturated lateral flow. The solute plots show that the WECC plume is more dispersed: due mainly to solute transport by induced numerical dispersion in the WEC-C model, though some of the differences are related to predictions of the SWMSF3D model. As configured, the SWMSF3D model produced mild oscillations in that plume solute concentrations were slightly above the inflow concentration of 1. The effect is most marked in the 20-day output (Fig. 7 (a)), where there is an unrealistic cavity at the front of the solute plume. Because, in any case, the WEC-C model is not designed to simulate fine details of water and solute dynamics, it was not considered necessary to experiment further with the SWMSF3D simulations to remove the observed oscillations. In conclusion, while the results of this test imply that the WEC-C model can perform to a reasonable level on such an exercise, it is clear that the shape of the plume is highly dependent on the grid length scales in the model. For instance, if the layer thickness is reduced in the upper layers from 2 m to 1 m, the plume becomes narrower and longer. The 0.1 concentration isohyet reaches the 10 m level after 100 days compared with the 19 m level in Fig. 7(b). This is because the thinner layers tend not to develop the saturation necessary to initiate lateral flow in the layers above the water-table. In terms of the WEC-C model’s intended applications, we note that this type of lateral unsaturated flow is at a scale that is much smaller than the typical cell size of a catchment scale model and would therefore have no relevance to a catchment simulation. This test does, however, clearly show some limitations of the WEC-C model and highlights the need to limit its application to the problems for which it is intended.

597

Fig. 7. (a) Comparison of SWMSF3D and WEC-C results for t=20 days. (b) Comparison of SWMSF3D and WEC-C results for t=100 days. (c) Comparison of SWMS-3D and WEC-C results for t=500 days.

598

J.T. Croton, D.A. Barry / Environmental Modelling & Software 16 (2001) 583–599

5. Discussion The guiding principle in the model’s development has been to maintain simplicity while still considering all the major elements of the hydrology and salinity processes of a catchment. This principle of simplicity has led to the use of explicit solvers, which are in the main only first-order accurate, and to the employment of operator splitting between the lateral and vertical components of the model. The model structure is in the form of a flexible framework which permits: the easy inclusion of components such as a dual continuum for the simulation of preferred pathway vertical fluxing; direct linkage between the vertical and lateral solvers; and considerable freedom with regard to definition of external fluxes. Also, the model is stable regardless of the form of the soil moisture characteristic curves used; and the lateral model component tolerates combinations of “wet” cells (cells with positive head) and “dry” cells (cells with zero head), that cannot easily be solved with implicit schemes. Simple explicit solvers do, however, suffer from a lack of accuracy. To determine the significance of this in the WEC-C model, four test problems were developed. Test 1 showed that the lateral solver in WECC was able to provide accurate results in terms of heads, at least for simple flow cases. Test 2 showed that the lateral model was also accurate for simple cases of transport and that the numerical dispersion within the lateral model was consistent with that due to the mixing-cell approach. Test 3 assessed the ability of the vertical model to simulate the progress of a wetting front into dry soil, and it was found that the general features of the infiltration front were captured exceptionally well considering the coarse layering used. In addition, it is worth noting that WEC-C employs a dual continuum approach that allows a closer representation of actual infiltration processes than does a single continuum approach, though parameterisation remains an issue. The fourth test reviewed the model’s ability to simulate both flow and solute transport for a three-dimensional case. It gave an accurate impression of the general performance of WEC-C in toto. While the variations between the WEC-C model results and those for the “exact” solution by SWMSF3D could be minimised to an acceptable level, it was noted that there was length scale dependence in the WEC-C results. For instance, if the layer thickness was reduced in the upper layers from 2 m to 1 m, the plume became narrower and longer. This was primarily related to lateral flow in WEC-C being dependent on the generation of saturated conditions, while in SWMSF3D it is true unsaturated lateral flow. It was noted, however, that this type of lateral unsaturated flow is normally at a scale much smaller than the cell size of a catchment scale model and would therefore have no relevance to a normal WEC-C catchment simul-

ation. This test did, though, clearly show some limitations of the WEC-C model, emphasising that it should be applied only to the type of problems, and at the length scales, for which it is intended.

6. Conclusion The concept of a distributed deterministic catchment model is not new; there are many examples in the literature. However, the WEC-C model contains a number of innovations. For example, it allows for the imposition of catchment topological changes resulting from land uses such as mining or construction of surface drains. It features preferential vertical flow in the vadose zone, modelled using a dual-continuum approach, which although not tested here, has been shown by model application to be very important. Numerically, it differs from many existing models in that it employs a direct linkage between the vertical and lateral solvers of its split-solver scheme and, due to its use of explicit solvers, is stable regardless of the form of the soil moisture characteristics. The WEC-C model framework of a rectangular grid of uniform cell size in the lateral plane, combined with a system of soil layers in the vertical, creates a simple model structure into which catchment variability can be incorporated by defining each parameter locally. When such a structure is combined with explicit solvers for computation of both the vertical and lateral flow and transport components, a conceptually realistic but simple and numerically robust model is produced. To assess the validity of this approach, four tests were undertaken where the model was compared against exact solutions. It was found that while there were limitations within the solvers, in particular the length scale dependence of the generation of unsaturated lateral flows and the level of induced numeric dispersion for solutes in the vertical solver, the model performed well and should be applicable for its intended role as a catchment model. Overall, we conclude that the model is a useful tool for the study of catchment hydrology in situations where groundwater discharge, saturated interflow and direct surface runoff are all components of the streamflow generation process.

Acknowledgements The authors express their gratitude to Alcoa World Alumina Australia who have sponsored the model’s development; and to the Water and Rivers Commission of Western Australia, particularly Mohammed Bari and Geoff Mauger who are working with the authors on application of the model to catchment studies. Geoff Mauger has also developed the interface between the MAGIC system and the WEC-C model. Thanks also to

J.T. Croton, D.A. Barry / Environmental Modelling & Software 16 (2001) 583–599

Jenny Dalton for technical editing of the manuscript and to Ken McIntosh and Craig Beverly for their help and advice during the development of the WEC-C model.

References Abbott, M.B., Bathurst, J.C., Cunge, J.A., O’Connell, P.E., Rasmussen, J., 1986. An introduction to the Europe Hydrological System — Systeme Hydrologique Europeen, “SHE” 2: structure of a physically-based, distributed modelling system. J. Hydrol. 87, 61–77. Bajracharya, K., Barry, D.A., 1994. Note on common mixing cell models. J. Hydrol. 153, 189–214. Barry, D.A., Miller, C.T., Culligan-Hensley, P.J., 1996a. Temporal discretisation errors in non-iterative split-operator approaches to solving chemical reaction/groundwater transport models. J. Contam. Hydrol. 22, 1–17. Barry, D.A., Barry, S.J., Parlange, J.-Y., 1996b. Capillarity correction to periodic solutions of the shallow flow approximation. In: Pattiaratchi, C.B. (Ed.), Mixing Processes in Estuaries and Coastal. Seas, Coastal and Estuarine Studies, vol. 50. American Geophysical Union, Washington, DC, pp. 496–510. Barry, D.A., Parlange, J.-Y., Sander, G.C., Sivapalan, M., 1993. A class of exact solutions for Richards’ equation. J. Hydrol. 142, 29–46. Barry, D.A., Sander, G.C., 1991. Exact solutions for water infiltration with an arbitrary surface flux or nonlinear solute adsorption. Water Resour. Res. 27, 2667–2680. Bear, J., Verruijt, A., 1987. Modelling Groundwater Flow and Pollution. D. Reidel Publ. Co, Dordrecht. Beven, K.J., 1996. A discussion of distributed hydrological modelling. In: Abbott, M.B., Refsgaard, J.C. (Eds.), Distributed Hydrological Modelling., pp. 255–278. Broadbridge, P., White, I., 1988a. Constant rate rainfall infiltration: A versatile nonlinear model, 1. Analytic solution. Water Resour. Res. 24, 145–154. Broadbridge, P., White, I., 1988b. Reply. Water Resour. Res. 24, 2109–2110. Dagan, G., 1967. Second-order theory of shallow free-surface flow in porous media. Q. J. Mech. Appl. Math. 20, 517–526. Calver, A., Wood, W.L., 1995. In: Singh, V.P. (Ed.), Computer Models of Watershed Hydrology. Water Resources Publications, pp. 595–626. Davidson, M.R., 1985. Numerical calculation of saturated–unsaturated infiltration in a cracked soil. Water Resour. Res. 21, 709–714. Duguid, J.O., Lee, P.C.Y., 1977. Flow in fractured porous media. Water Resour. Res. 13, 558–566. Edwards, W.M., van der Ploeg, R.R., Ehlers, W., 1979. A numerical study of the effects of non-capillary size pores upon infiltration. Soil Sci. Soc. Am. J. 43, 851–856. Fleming, G., 1975. Computer Simulation Techniques in Hydrology. Elsevier, New York. Gerke, H.H., van Genuchten, M.T., 1993. A dual-porosity model for simulating the preferential movement of water and solutes in structured porous media. Water Resour. Res. 29, 305–319. Grayson, R.B., Moore, I.D., McMahon, T.A., 1992. Physically based hydrologic modelling 2. Is the concept realistic? Water Resour. Res. 28, 2659–2666. Haverkamp, R., Vauclin, M., 1979. A note on estimating finite differ-

599

ence interblock hydraulic conductivity values for transient unsaturated flow problems. Water Resour. Res. 15, 181–188. Jain, S.K., Storm, B., Bathurst, J.C., Refsgaard, J.C., Singh, R.D., 1992. Application of the SHE to catchments in India 2. Field experiments and simulation studies with the SHE on the Kolar subcatchments of the Narmada River. J. Hydrol. 140, 25–47. Lapidus, L., Amundson, N.R., 1952. Mathematics of adsorption in beds. VI. The effects of longitudinal diffusion on ion exchange and chromatographic columns. J. Phys. Chem. 56, 984–988. Marshall, J.K., 1993. Yearlong water uptake by jarrah (Eucalyptus marginata) in intact jarrah forest. In: Invited paper to International Symposium on Forest Hydrology, Water Issues in Forests Today, Canberra. Marshall, J.K., Chester, G.W., 1992. Effect of forest thinning on Jarrah (Eucalyptus marginata) water uptake. Land and Water Research News 14, 24–27. Miller, C.T., Rabideau, A.J., 1993. Development of split-operator, Petrov–Galerkin methods to simulate transport and diffusion problems. Water Resour. Res. 29, 2227–2240. Noye, J., 1982. Numerical Solutions of Partial Differential Equations. North Holland, Amsterdam. O’Loughlin, E.M., 1986. Prediction of surface saturation zones in natural catchments by topographic analysis. Water Resour. Res. 22, 794–804. Parlange, J.Y., Stagnitti, F., Starr, J.L., Braddock, R.D., 1984. Freesurface flow in porous media and periodic solution of the shallowflow approximation. J. Hydrol. 70, 251–263. Pinder, G.F., Gray, W.G., 1977. Finite Element Simulation in Surface and Subsurface Hydrology. Academic Press, New York. Pruess, K., Wang, J.S.Y., Tsang, Y.W., 1990. On thermohydrologic conditions near high-level nuclear wastes placed in partially saturated fractured tuff. 1. Simulation studies with explicit consideration of fracture effects. Water Resour. Res. 26, 1235–1248. Richards, L.A., 1931. Capillary conduction of liquids through porous mediums. Physics 1, 983–986. Sander, G.C., Parlange, J.-Y., Ku¨ nhel, V., Hogarth, W.L., Lockington, D., O’Kane, J.P.J., 1988a. Exact nonlinear solution for constant flux infiltration. J. Hydrol. 97, 341–346. Sander, G.C., Parlange, J.-Y., Ku¨ nhel, V., Hogarth, W.L., O’Kane, J.P.J., 1988b. Comment on “Constant rate rainfall infiltration: A versatile nonlinear model 1. Analytic solution” by P. Broadbridge and I. White. Water Resour. Res. 24, 2107–2108. Sharma, M.L., Barron, R.J.W., Fernie, M.S., 1987. A real distribution of infiltration parameters and some physical properties in lateritic catchments. J. Hydrol. 94, 109–127. Simunek, J., Huang, K., Van Genuchten, M.Th., 1995. The SWMSF3D code for simulating water flow and solute transport in three-dimensional variably-saturated porous media. Version 1.0. Rep. 139, US Salinity Lab., Riverside, CA. Singh, V.P., 1995. Watershed modelling. In: Singh, V.P. (Ed.), Computer Models of Watershed Hydrology. Water Resources Publications, pp. 1–22. Strang, G., 1968. On the construction and comparison of difference schemes. SIAM J. Numer. Anal. 5, 506–517. Terzaghi, K., Peck, R.B., 1967. In: Soil Mechanics in Engineering Practice. Wiley, New York, p. 729 (reprinted 1967). Whisler, F.D., Klute, A., Millington, R.J., 1970. Analysis of radial, steady-state solution, and solute flow. Soil Sci. Soc. Amer. Proc. 34, 382–387. White, I., Broadbridge, P., 1988. Constant rate rainfall infiltration: A versatile nonlinear model, 2. Applications of solutions. Water Resour. Res. 24, 155–162.