Published in Vadose Zone Journal 4:206-221 (2005)
© 2005 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
A Set of Analytical Benchmarks to Test Numerical Models of Flow and Transport in Soils
J. Vanderborghta,*,
R. Kasteela,
M. Herbsta,
M. Javauxb,
D. Thiéryc,
M. Vancloosterb,
C. Mouvetc and
H. Vereeckena
a Agrosphere Institute, ICG-IV, Forschungszentrum Jülich GmbH, 52425 Jülich, Germany
b Department of Environmental Sciences and Land Use Planning, Université Catholique de Louvain (UCL), Croix du Sud, 2 Bte 2, B-1348 Louvain-la-Neuve, Belgium
c Bureau de Recherches Geologiques et Minieres, BRGM, Avenue Claude Guillemin, F 45060 Orléans Cedex 02, France
* Corresponding author (j.vanderborght{at}fz-juelich.de)
Received 13 May 2004.
 |
ABSTRACT
|
|---|
Accurate model predictions of water flow and solute transport in unsaturated soils require a correct representation of relevant mechanisms in a mathematical model, as well as correct solutions of the mathematical equations. Because of the complexity of boundary conditions and the nonlinearity of the processes considered, general solutions of the mathematical equations rely on numerical approximations. We evaluate a number of numerical models (WAVE, HYDRUS_1D, SWAP, MARTHE, and MACRO) that use different numerical methods to solve the flow and transport equations. Our purpose is to give an overview of analytical solutions that can be found for simple initial and boundary conditions and to define benchmark scenarios to check the accuracy of numerical solutions. Included are analytical solutions for coupled transport equations that describe flow and transport in dual-velocity media. The relevance of deviations observed in the analytical benchmarks for more realistic boundary conditions is illustrated using an intercode comparison for natural boundary conditions. For the water flow scenarios, the largest deviations between numerical models and analytical solutions were observed for the case of soil limited evaporation. The intercode differences could be attributed to the implementation of the evaporation boundary condition: the spatial discretization and the internode averaging of the hydraulic conductivity in the surface grid layer. For solute transport, accurate modeling of solute dispersion poses the most problems. Nonlinear and nonequilibrium sorption and coupled transport in pore domains with different advection velocities are in general accurately simulated.
Abbreviations: BTC, breakthrough curve • CDE, convection dispersion equation
|
INTRODUCTION
|
|---|
MODELS THAT PREDICT flow and transport processes in soils are increasingly being applied to address practical problems. Examples are the assessments of the water balance in engineered covers of waste disposal sites, of leaching of pesticides and fertilizers through the soil to surface near groundwater, and of water and salt balances of irrigated fields. The use of simulation models is triggered by legislation that prescribes the use of pesticide fate simulation models in the registration procedure of new pesticides (e.g., the EU-directive 91/414; Council of the European Union, 1997). The use of models allows extrapolation in time and space of data from leaching experiments and monitoring campaigns. A model is by definition a simplification of complex reality. To have a predictive model, all processes and properties of the system that have a decisive impact on the system's behavior and functioning must be represented. To test the capability of a model to represent reality, model simulations must be compared with experimental data and with simulations by other models. Scanlon et al. (2002) presented an intercode comparison of seven codes to simulate the water balance of nonvegetated covers. Differences in model simulations were attributed to differences in model concepts (Richards' equation vs. storage routing), definition of upper and lower boundary conditions, and parameterization of the water retention curve. Examples of intercode comparisons of pesticide leaching models were presented by Pennell et al. (1990) and Bergström and Jarvis (1994). Besides differences in model concepts, the parameterization of models and their use by different modelers may also lead to important differences in simulation results (Vanclooster et al., 2000). For leaching processes, nonequilibrium phenomena in the soil may lead to a behavior different from what is expected based on equilibrium assumptions. Nonequilibrium can be caused by rate-limited sorption or by an incomplete mixing of percolating with resident waters. A comparison between different concepts to model nonequilibrium flow in soils was given by Simunek et al. (2003). These authors concluded that there is still a need for comparisons of codes and concepts simulating nonequilibrium processes in soils. In particular, experimental data sets that can be used for such comparisons are missing.
Given the nonlinearity of the governing equations, the continuously changing boundary conditions at the soil surface, and the vertical variations of soil properties, simulation of flow and transport relies on a numerical solution of the flow and transport equations. In model and code comparisons that were done so far, models were tested against each other and against experimental data. It was implicitly assumed that the numerical models solved these partial differential equations with sufficient accuracy. The outcome of the study by Vanclooster et al. (2000), showing "user-specific" variability of model simulations, indicates that numerical aspects may play an important role. To test the accuracy of numerical solutions, these solutions can be compared with analytical solutions of flow and transport equations that exist for given simplified initial and boundary conditions. Our objective in this paper is to present a set of analytical solutions and to define a set of benchmark simulation scenarios that can be used to test numerical models. All the numerical models that we considered use the one-dimensional Richards and convection dispersion equations (CDE) to describe flow and transport in soils. To describe nonequilibrium flow and transport, we also considered nonequilibrium sorption and dual permeability models. These equations imply certain assumptions about the processes and make certain simplifications of the complex reality, which we will not discuss in this paper. An overview of flow and transport models and the underlying process assumptions can be found in Feyen et al. (1998). Here we focus on the accuracy or precision with which the numerical flow and transport models solve the flow and transport equations. Therefore, we will make the assumption that these equations represent reality exactly. As a consequence, the notion of "correct prediction" should be interpreted in this paper as a precise numerical solution of an equation.
Because the models are based on the same flow and transport concepts and use the same functions to describe the unsaturated hydraulic soil properties, differences among model simulations can be attributed to differences in the numerical solution. In the Theory section, analytical solutions of the flow and transport equations are presented and listed according to mathematical transformation that is used to obtain the solution. In Benchmarks, the numerical models we consider are presented, and the benchmark scenarios are defined. A set of benchmarks, for which analytical solutions exist, and a benchmark using realistic boundary conditions are considered. The differences among the model simulations for the realistic scenario are then interpreted in terms of the results obtained from the first set of benchmarks.
|
THEORY
|
|---|
Flow and Transport Equations
The vertical water flux, Jw (L T–1) in an unsaturated soil is generally described by the Buckingham–Darcy equation:
 | [1] |
where x (L) is the vertical coordinate, which is chosen to be positive in the downward direction, 
(L) is the pressure head, and K() (L T–1) is the hydraulic conductivity. Using Eq. [1] in a mass balance, the one-dimensional Richards equation is obtained:
 | [2] |
where 
is the volumetric water content. Both and K are highly nonlinear functions of . Equation [2] is the so-called mixed form of the Richards equation, with both and as dependent variables. To obtain analytical solutions, the "-based" form of the Richards equation may be used:
 | [3] |
where Dw() = K()(d/d) (L2 T–1) is the water diffusivity. To solve the Richards equation (Eq. [2] and [3]), the hydraulic functions (), K(), and Dw() need to be defined. The van Genuchten–Mualem parameterization (van Genuchten, 1980) is frequently used to describe experimentally determined (), K(,), Dw() functions
 | [4] |
 | [5] |
 | [6] |
where s is the saturated and r the residual volumetric water content; 
= ( – r)/(s – r); Ks (L T–1) the saturated hydraulic conductivity; and 
(L–1), n, and l are empirical constants describing the shape of the curves.
Transport is simulated using the CDE:
 | [7] |
where C (M L–3) is the solute concentration in the water phase, v = Jw/ (L T–1) is the pore water velocity, Ds (L2 T–1) is the hydrodynamic dispersion coefficient, which can be assumed to be independent of the solute concentration, 
b (M L–3) is the soil bulk density, and S (M M–1) is the sorbed solute mass per unit mass of dry soil. The sorption isotherm, F, expresses the relation between the sorbed and dissolved concentration under equilibrium conditions, Seq and Ceq:
 | [8] |
We consider a Freundlich isotherm:
 | [9] |
where kf 
is the Freundlich sorption coefficient, and nf is the Freundlich exponent. When sorption and desorption occur much faster than solute concentration changes in the water phase due to advective and dispersive fluxes, the equilibrium sorption isotherm Eq. [8] is used to write S as a function of C in Eq. [7]. For a uniform soil profile and steady-state flow conditions (, v, Ds, and F are independent of depth and time) and equilibrium sorption, Eq. [7] can be written as
 | [10] |
where the retardation coefficient R(C) is
 | [11] |
When sorption or desorption are rate-limited, a first-order kinetic is commonly considered, and transport in a uniform soil profile and steady-state-flow are described by the following set of equations:
 | [12] |
where s (T–1) is the first-order kinetic sorption rate coefficient. The kinetic sorption model may be extended by considering a set of sorption sites with different sorption rates. For porous media with two separate pore networks, which have different transport properties (e.g., fractured media, soil with macropores or interaggregate pores embedded in a microporous soil matrix or soil aggregates), a set of coupled Richards and CDE equations may be used to describe transport (e.g., Gerke and van Genuchten, 1993). For a vertically uniform soil profile and steady-state flow, the coupled set of CDE equations is
 | [13] |
where 
(T–1) is the first-order exchange coefficient, and the subscripts refer to the two pore regions.
Analytical Solutions
We grouped analytical solutions of the flow and transport equations according to the transformation method that is used to obtain the solution.
Kirchhoff Transform
To predict pressure head profiles for steady-state flow in layered soil profiles, that is, steady-state downward (infiltration) or upward (evaporation) flow to or from a shallow groundwater table, the Kirchhoff transform method is used. Integration of Eq. [1] yields depth as a function of pressure head:
 | [14] |
Laplace Transform
To solve the linear transport equations (Eq. [10], [12] with nf = 1) in semi-infinite soil profiles, the Laplace transform method is used (e.g., Jury and Roth, 1990; Toride et al., 1993). The initial g(x) and boundary h(t) conditions are defined as
 | [15] |
The concentrations that are predicted using the boundary conditions in Eq. [15] are volume-averaged or resident concentrations. However, concentrations that are measured in the effluent from a soil column are flux-averaged concentrations, Cf = Js/Jw (Js = solute flux including advective and dispersive fluxes). Cf is related to the resident concentration, C, as (Parker and van Genuchten, 1984):
 | [16] |
For linear transport, the solution for general initial g(x) and boundary h(t) conditions can be expressed in terms of transfer functions: cI(x,x';t), sI(x,x';t) and cB(t;x), sB(t;x), which are solutions of Eq. [10], [12] for, respectively, g(x) = 
(x – x') and h(t) = 0, and for g(x) = 0 and h(t) = (t):
 | [17] |
 | [18] |
Since Cf is linearly related to C, Cf can be expressed using Eq. [17] in terms of transfer functions cfI
and cfB
. These transfer functions may be derived directly from the transfer functions cI(x,x';t) and cB(t;x) through Eq. [16]. Alternatively, cfB may be obtained by solving Eq. [10] or [12] for a first-type boundary condition: C(x = 0,t) = h(t).
The solution cB(t;x), sB(t;x) is obtained by taking Laplace transforms of Eq. [10], [12], and of the boundary condition that leads to an ordinary homogeneous differential equation. The transfer functions cI(x,x';t), sI(x,x';t) for an initial concentration distribution are obtained by taking Laplace transforms with respect to time and depth. A complete set of transfer functions that are obtained for Eq. [10] and [12] were given by Toride et al. (1993). Below, we give the transfer functions that we used in the benchmark study. The transfer functions for the "equilibrium" transport equation, Eq. [10] are (Toride et al., 1993)
 | [19] |
 | [20] |
 | [21] |
 | [22] |
The transfer functions for the "nonequilibrium" transport equation, Eq. [12], can be expressed in terms of the transfer functions of the CDE without sorption (i.e., R = 1 in Eq. [19]–[22]) (Toride et al., 1993):
 | [23] |
 | [24] |
 | [25] |
 | [26] |
where
 | [27] |
and Ij is the modified Bessel function of the first kind and order j. The transfer functions for the flux concentrations are analogous and are obtained by replacing ceq by cfeq.
Analytical solutions of the set of coupled CDE equations (Eq. [13]) were discussed by Walker (1987). For general parameters, Eq. [13] can be solved for an initial concentration distribution in an infinite soil profile:
 | [28] |
where Mi is the mass per unit area in pore domain i. The following solution is obtained (Walker, 1987):
 | [29] |
 | [30] |
where I0 is the modified Bessel function of the first kind and zero order, and w(x,t) is the solution of the CDE for the case that no solute transfer occurs between the two regions:
 | [31] |
and where
 | [32] |
 | [33] |
 | [34] |
 | [35] |
 | [36] |
Boltzmann Transform
Introducing the Boltzmann variable 
= xt–0.5, and neglecting the second term on the right-hand side of Eq. [3] (the gravitation term), the partial differential equation can be transformed to an ordinary differential equation:
 | [37] |
For semi-infinite soil columns with a constant water content at the soil surface, ( = 0) = sur and a uniform initial water content, ( = 
) = i, the Boltzmann transformed flow equation can be used to predict the water flux at the soil surface: Jw(x = 0,t) when gravitational forces can be neglected (at the early stage of an infiltration in an initially dry soil or during evaporation from an initially uniform wet soil). Since the water content profiles can be written as a function of , it follows directly from the mass balance that
 | [38] |
where Sw (L T–1/2) is the sorptivity (infiltration) or desorptivity (evaporation).
For evaporation from an initially wet soil (sur < i), Lockington (1994) proposed the following relation between the desorptivity and water diffusivity Dw():
 | [39] |
where
 | [40] |
 | [41] |
 | [42] |
 | [43] |
To model evaporation from a soil surface, the surface flux is set equal to the potential evaporation rate, Jwpot, until a critical soil pressure head, crit, or water content at the surface is reached crit. Then, the pressure head or soil water content remains constant and the evaporation flux decreases with time when the soils dries out further. The time period of potential evaporation, tpot, can be approximated by the Time Compression Analysis. It is assumed that the water content profile at tpot that results from a constant water flux at the surface matches the one that is obtained for a constant sur = crit when the surface flux reaches Jwpot. This implies that I(t') (i.e., the cumulative evaporation from the soil profile with sur = crit at the time t' when the surface flux equals Jwpot) is equal to Jwpottpot (i.e., the cumulative evaporation from the profile with constant evaporation rate at time tpot when the surface water content reaches crit) Using Eq. [38], t' and tpot are derived as
 | [44] |
 | [45] |
The flux at the soil surface is
 | [46] |
 | [47] |
Traveling Wave Solutions
For certain situations, traveling wave solutions of concentration or water content profiles that keep a constant shape when they propagate through the soil can be obtained:
 | [48] |
For solute transport in a homogeneous soil profile, a traveling wave emerges for the following conditions:
 | [49] |
Because of the decrease of the retardation coefficient with increasing C, higher concentrations propagate faster than lower ones. If there is no solute dispersion, a shock front results:
 | [50] |
where H is the Heaviside function, and
 | [51] |
If Ds > 0, dispersion will smooth out the shock front until equilibrium between smoothing, due to dispersion, and sharpening, due to nonlinear sorption, is reached. Using the following transform:
 | [52] |
the CDE Eq. [10] is transformed to an ordinary differential equation:
 | [53] |
For relatively long travel times, the effect of the surface boundary on the transport can be neglected, and the semi-infinite profile can be approximated by an infinite soil profile with the following boundary conditions:
 | [54] |
For the case of a Freundlich isotherm (Eq. [9]) with nf = i/(i + 1), where i is a positive integer, the following expression is obtained (van der Zee, 1990):
 | [55] |
where s = Ds/v, the dispersivity.
Since [d2K()]/d2 > 0, it follows from a comparison between the flow and transport equations (Eq. [3] and [10]) that a traveling water front will emerge during infiltration in a dry soil. For the following boundary condition:
 | [56] |
and using the following transform:
 | [57] |
the Richards equation (Eq. [3]) is transformed to
 | [58] |
If we consider again an infinite soil profile (i.e., the boundaries are sufficiently far from the infiltration front), the following boundary conditions can be defined:
 | [59] |
so that after integration, the following traveling wave solution is obtained (Philip, 1969):
 | [60] |
where a is a reference water content and 
() the position of the front relative to the position of the reference water content.
|
BENCHMARKS
|
|---|
Numerical Models
Table 1 lists the numerical models that we considered in the benchmark study, some of their main features, and their references. The models were chosen to cover a range of different numerical approximations of the flow and transport equations. The MACRO, SWAP (which is a part of the PEARL model), and WAVE models are used in European (FOCUS, 2000) and national procedures for the registration of plant protection products. The HYDRUS_1D model is a commonly used model in the scientific community. All the models can simulate rate-limited sorption–desorption and, except for the WAVE model, nonlinear sorption. The MARTHE model solves the three-dimensional flow and transport equations in a full continuum and is used to predict flow and transport processes in coupled soil and aquifer systems. Three models, MACRO, MARTHE, and DUALP_1D, can solve coupled flow and transport equations (Eq. [13]) that describe flow and transport in dual-porosity media. The DUALP_1D is based on the HYDRUS_1D code; therefore, we refer to it as HYDRUS_1D. In HYDRUS_1D and MARTHE, flow and transport in micro- and macropore regions are described by the Richards' equation and CDE, respectively. More technical details about the DUALP_1D are given by Tseng et al. (1995). In MACRO, a kinematic wave approach neglecting capillary forces and a purely convective transport neglecting dispersion are considered to model flow and transport in the macropores. Flow and transport in the micropore region are modeled by the Richards equation and CDE. Both models use a first-order mass exchange term to model mass exchange between both regions. For a comparison between both model approaches, we refer to Simunek et al. (2003).
All considered models use an implicit scheme to discretize the one-dimensional Richards equation. In the finite element scheme of the HYDRUS_1D model, mass lumping is implemented to the time derivative term so that the finite element scheme leads to a standard finite difference scheme:
 | [61] |
where the subscript i refers to the spatial position, the superscript j to the time level, and k to the iteration level; tj = tj+1 – tj, x = xi+1 – xi–1, and xi = xi+1 – xi. All models use a first-order expansion of the time derivative coupled with a Picard iteration to assure mass conservation (Celia et al., 1990):
 | [62] |
Using Eq. [62], Eq. [61] is solved for j+1,k+1i until the solution converges. Following Huang et al. (1996), who showed that a -based convergence criterion is more efficient than a -based criterion, all models check the convergence of j+1,k+1i and j+1,ki.
The differences between the models are ascribed to definition of the hydraulic conductivity between the spatial points where the pressure head is calculated, Ki+1/2j+
. First, Ki+1/2j+ may be considered explicitly at the previous time level ( = 0) (MACRO, SWAP) or implicitly at the new time level ( = 1) (HYDRUS_1D, MARTHE, WAVE). Second, Ki+1/2j+ may be the arithmetic: 
/2 (SWAP, HYDRUS_1D, MACRO), or geometric 
0.5 (WAVE) average of the hydraulic conductivity at the nodes. In the MARTHE model, the averaging scheme is user defined. Besides the arithmetic or geometric average, the harmonic average or the upstream conductivity also can be selected. The geometric average was used for the MARTHE simulations.
The differences between the models are more pronounced for the transport equation. The discretized form of the transport equation, Eq. [7], may be written in general as
 | [63] |
where the superscripts refer to the time level and the subscripts to the spatial position. U(C) is a sink–source term that describes rate-limited sorption or exchange of solutes between different flow domains (Eq. [12] and [13]). In the SWAP model, the right-hand side of Eq. [63] is explicitly evaluated at time level j ( = 0). In the MACRO and WAVE models, the average of the right-hand side term at time level j and j + 1 is considered ( = 0.5, Crank-Nicholson scheme). But, the WAVE model treats the sink term explicitly [U(C) is evaluated at time level j]. In the MARTHE model, the total variation diminishing (TVD) scheme is used. The dispersivity term in the right-hand side of Eq. [63] is treated implicitly, but the advection term explicitly. The weight given to the time levels in the HYDRUS_1D model can be defined by the user. We used a Crank–Nicholson scheme for the simulations with the HYDRUS_1D model. To approximate the spatial derivatives, a Galerkin finite element scheme is used in the HYDRUS_1D model and a finite difference scheme in the other models. SWAP and MACRO use a first-order finite difference scheme whereas second-order terms are included in the WAVE model. In the SWAP model, a central finite difference scheme is used to discretize the right-hand side of Eq. [63]. In the MACRO and WAVE models, upstream weighting or backward finite differences are used to discretize the advection term of the right-hand side of Eq. [63]. Upstream weighting can be selected by the user in the HYDRUS_1D model but was not used in the simulations. Because upstream weighting leads to artificial numerical dispersion, an empirically determined correction term is subtracted from the dispersion coefficient in the MACRO and WAVE models. In the MARTHE model a third-order backward interpolation is used to evaluate the advection term. To avoid numerical oscillations, a flux limiter that prescribes a monotonous variation is applied.
Scenarios
A set of flow and transport scenarios were simulated. In the scenarios, we considered three different standard soil materials. The hydraulic parameters of these soils are listed in Table 2 and were based on hydraulic parameters for sand, loam, and clay texture classes in the soil catalog of the HYDRUS_1D software (Simunek et al., 1998). Unless specified differently, numerical simulations were done in a 2-m-deep soil profile using a spatial discretization of 1 cm. In a first set of scenarios, the boundary and initial conditions were defined so that numerical solutions of the flow and transport equations could be compared with analytical solutions. To solve the transport equations numerically, a zero concentration gradient was defined at the bottom of the soil profile in all transport scenarios. Details of the different benchmark scenarios are given in Table 3. The agreement between the numerical and analytical solutions is qualitatively evaluated using a graphical presentation of the results. As a quantitative measure of the agreement, the coefficient of determination, R2 is calculated:
 | [64] |
where pi and mi are the ith elements of the numerical and analytical model output vectors, respectively, and N is number of entries in the vector. The coefficients of determination are provided in the graphs as additional information on the quality of the numerical simulations.
In a second set of scenarios, a mixed boundary condition with realistic water fluxes at the soil surface was used. For –100000 cm = crit < (x = 0) < 0, the flux at the soil surface Jw(x = 0) was defined as the difference between the precipitation and evaporation rates, (Prec – E0). Otherwise, the pressure head at the surface was defined: (x = 0) = crit or (x = 0) = 0, and kept constant as long as (Prec – E0)/Jw(x = 0) > 1. A 2-yr record of daily Prec and E0 was used to define the soil surface boundary conditions (see Fig. 1)
. To determine the initial pressure head profile, we ran the HYDRUS_1D model for the loamy soil using an initially uniform pressure head, (x,t = 0) = –200 cm, and the 2-yr record of Prec and E0. The pressure head profile that was obtained at the end of this simulation was used as the initial pressure head profile (see Fig. 1). The initial concentration distribution was defined as: C(x,t = 0) = 100 for 0 < x < 20 and C(x,t = 0) = 0 otherwise.
View larger version (24K):
[in this window]
[in a new window]
|
Fig. 1. (a) Precipitation and potential evaporation and (b) initial pressure head, , depth profile used as input in the scenario with climatic boundary conditions.
|
|
|
RESULTS AND DISCUSSION
|
|---|
Water Flow Simulation
The pressure head profiles, (x), for steady-state downward water flow in layered soil profiles, and evaporation from a shallow water table, are shown in Fig. 2
and Fig. 3
, respectively. All numerical models approximate the analytical benchmark, Eq. [14], relatively well. For steady-state flow in a layered soil profile, the pressure head, , in the deeper layer is constant with depth. The deviations of the pressure head predictions by HYDRUS_1D from the analytical benchmark in the deeper soil layer (Fig. 2) are the result of the discrete approximation of the hydraulic functions K() and (). When the hydraulic functions are exactly evaluated, the HYDRUS_1D predictions coincide with the predictions by the other models (results not shown). Of note is also the sharp transition of the pressure head at the sand–loam interface (Fig. 2b). Except for WAVE, all other models predict a somewhat more gradual transition of the (x) profile. SWAP does not reproduce the (x) profile in the clay–sand layered soil (Fig. 2c). This may be explained by the -based convergence criterion to solve the nonlinear flow equation iteratively. In the clay soil, small changes in correspond with large changes in . As a consequence, a convergence criterion that is based on may lead to inaccurate profiles when the capacity C() = d/d is small, as in clayey soils or close to water saturation.
View larger version (14K):
[in this window]
[in a new window]
|
Fig. 2. Pressure head, , depth profiles in layered soils for a constant downward flow rate of 0.5 cm d–1: (a) loam–sand, (b) sand–loam, and (c) clay–sand profiles. (Black line is analytical benchmark, Eq. [14]; similarity in output may overlap in the graph.)
|
|
View larger version (23K):
[in this window]
[in a new window]
|
Fig. 3. Pressure head, , depth profile in a loamy soil with a water table at the 54-cm depth and a constant evaporation rate of 0.5 cm d–1. (Black line is analytical benchmark, Eq. [14]; similarity in output may overlap in the graph.)
|
|
The results for the transient water flow scenarios—infiltration in a dry soil and evaporation from an initially wet soil—are shown in Fig. 4
and Fig. 5
, respectively. Figure 4 clearly illustrates that during infiltration, the shape of the wetting front does not change with time. WAVE does not reproduce the infiltration fronts in the sandy soil (Fig. 4a) but tends to predict an increase in water content just before the wetting front. The overestimation of the water contents before the wetting front is compensated by a smaller infiltration depth and a steeper wetting front. Infiltration fronts in the clay soil (Fig. 4c) could be predicted by MACRO, WAVE, and MARTHE (using upstream weighted conductivities) whereas the other models (HYDRUS_1D and SWAP) did not converge for this scenario. Due to the strong nonlinearity of the hydraulic soil functions, the wetting fronts are very steep, but when convergence is obtained, the numerical solutions using a 1-cm spatial discretization reproduce the shape of the wetting fronts fairly well. However, models that use the arithmetic average of the nodal conductivities (HYDRUS_1D, MACRO, SWAP) as an estimate of the grid cell conductivity tend to predict more dispersed wetting fronts. Warrick (1991) illustrated that sharper wetting fronts are simulated when a geometric mean of the conductivity is used. The WAVE model simulations corroborate this. The more dispersed fronts simulated by MARTHE are caused by the upstream weighting. For sandy soils, the underestimation of the grid cell conductivity by the geometric mean may lead to an underestimation of the water flux in the leading edge of the wetting front and a corresponding accumulation of water further above.
View larger version (27K):
[in this window]
[in a new window]
|
Fig. 5. Evaporation rate, Eact, from initially wet (a) sandy, (b, c) loamy, and (d) clayey soil profiles. Dashed lines are simulated Eact using a spatial discretization of 1 cm, full lines using a discretization of 0.25 cm. (Black line is analytical benchmark Eq. [46] and [47]; similarity in output may overlap in the graph; R2 is calculated for simulations with a discretization of 0.25 cm.)
|
|
The deviations between the analytical benchmarks and the numerical simulations are clearly larger for the transient evaporation (Fig. 5) than for the infiltration scenarios (Fig. 4). The nonsmooth behavior of the SWAP predictions is due to the numerical precision of the output. HYDRUS_1D, MARTHE, and SWAP tend to predict a longer time period of potential evaporation than MACRO and WAVE. Since evaporation from a soil surface leads to a sharp increase of pressure head with depth close to the soil surface (see Fig. 3), a high spatial resolution is required to reproduce the pressure head profile at the soil surface. For HYDRUS_1D, MARTHE, and SWAP, an increase in spatial resolution, 0.25- instead of 1-cm cells, leads to a decrease of the simulated time of potential evaporation and a convergence to the analytical benchmark. This contradicts the results of van Dam and Feddes (2000), who found no significant differences between SWAP simulated evapotranspiration rates for 1-cm and 0.1-cm grid cells. For MACRO simulations in general and WAVE simulations in the loamy soil, an increase in spatial resolution results in a longer simulated time of potential evaporation. In the sandy soil (Fig. 5a), in which WAVE predicts almost no evaporation, and the clayey soil (Fig. 5d), the predictions by WAVE did not improve using a finer spatial discretization. The simulations by MACRO in the clayey soil were farther from the analytical benchmark for the finer than for the coarser discretization. But, they converged with the simulations by SWAP and MARTHE. The divergence from the analytical benchmark for a finer spatial discretization may be explained by the fact that the analytical benchmark, which is based on the Time Compression Analysis, is an approximation itself. Therefore, this analytical benchmark should best be compared with numerical simulations by different models using different spatial discretizations. The difference between WAVE and MACRO on one hand and HYDRUS_1D, SWAP, and MARTHE on the other can be attributed to the averaging of the hydraulic conductivity in a grid cell and the implementation of the flux boundary condition. Van Dam and Feddes (2000) showed that longer periods of potential evaporation are simulated when arithmetic rather than geometric averages of the conductivities in a grid cell are used. WAVE and MARTHE use geometric averages of the conductivities at the grid nodes whereas arithmetic averages are used in SWAP, HYDRUS_1D, and MACRO. In contrast to the other models, MARTHE calculates water fluxes in the middle of a grid cell. As a consequence, the surface water flux is applied to the middle of the first grid cell. For relatively thick cells, this approach leads to an overestimation of the duration of the potential evaporation stage. MACRO uses an arithmetic average of the hydraulic conductivity except for the conductivity of the surface compartment when the soil evaporates. Unlike the other models, the MACRO model does not switch to Dirichlet boundary condition when the surface water potential reaches the critical pressure head, crit. For a given pressure head in the surface layer, the surface potential that leads to the maximal evaporation rate is iteratively sought. The actual evaporation is calculated using Darcy's Law assuming a linear decrease of the total head between the center of the first layer and the soil surface and a constant hydraulic conductivity in that layer, which is calculated using the arithmetic mean of the pressure head in the first layer and at the soil surface. It should be noted that if the flux in the surface layer, which remains constant during a time step, were calculated exactly using the Kirchhoff transform, Eq. [14], then the maximum evaporation rate would be obtained for sur = – cm and would be larger than the evaporation rate that is obtained for a preset minimal sur = –10000 cm. Since the K() function is highly nonlinear, the arithmetic average of the conductivities at the surface and in the first layer is larger than the geometric average (WAVE) and larger than the conductivity evaluated at the averaged pressure head (MACRO).
Figure 6
shows the cumulative drainage from a 2-m-deep soil profile that is simulated by the different models for climatic boundary and initial conditions shown by Fig. 1. The differences in simulated total drainage are in line with the results from the transient evaporation benchmark. SWAP, HYDRUS_1D, and MARTHE predict a larger evaporative loss than MACRO and WAVE (Fig. 5). As a consequence, the cumulative net infiltration (dashed lines) predicted by MACRO and WAVE on one hand and SWAP, HYDRUS_1D and MARTHE on the other diverge during the spring and summer periods (from Day 150 to 350 and from Day 550 to 750). For the sandy soil, WAVE underestimates the evaporation losses considerably (Fig. 5a), which is in line with the larger predictions of the cumulative drainage by this model. A comparison between the cumulative net infiltration and the cumulative drainage in the sandy soil suggests mass balance errors in the simulations by the WAVE model. This may be related to the "overshooting" of the simulated water contents by this model at the wetting front during an infiltration event in a sandy soil (Fig. 4a). Because MACRO predicted the transient evaporation scenario better, the net cumulative infiltration and cumulative drainage may be better reproduced by this model. However, the differences between the predictions are relatively small.
View larger version (19K):
[in this window]
[in a new window]
|
Fig. 6. Simulated cumulative drainage at the 2-m depth (full lines) and cumulative net infiltration (precipitation – actual evaporation) (dashed lines) in (a) sandy, (b) loamy, and (c) clayey soil profiles. (Full black line is the cumulative precipitation; similarity in output may overlap in the graph.)
|
|
Solute Transport Simulations
Simulated concentration depth profiles and breakthrough curves for a steady-state downward flow are shown in Fig. 7
. Due to the explicit central finite difference scheme used by SWAP to solve the transport equation, this model does not accept dispersivity values that are smaller than x/2. Therefore, the s = 0.1 cm scenario (Fig. 7a, 7d) could not be simulated by SWAP. WAVE clearly overestimated the dispersion and did so also for relatively large dispersion coefficients. The predictions by all the other models agreed fairly well with the analytical benchmark, even for the small dispersivity case with a relatively large grid Peclet number, x/s. Small differences between HYDRUS_1D and the analytical benchmark are observed. These are due to a smaller estimate of the soil water content for the applied water flux by HYDRUS_1D, which results from the discrete approximation of the soil hydraulic functions. Again, these deviations disappear when the hydraulic functions are exactly evaluated. The effect of the zero concentration gradient bottom boundary on the numerical solution of the CDE and the deviation from the analytical benchmark, which assumes a semi-infinite profile, is apparent in the concentration depth profile for the = 10 cm scenario (Fig. 7c). However, the effect of the bottom boundary on the breakthrough curve (BTC) is small. A zero concentration bottom boundary condition implies that flux and resident concentrations are identical at the bottom boundary. Therefore, the deviation between analytical solutions in semi-infinite profiles and numerical solutions in bounded domains will be more important for smaller column Peclet numbers, L/s (van Genuchten and Parker, 1984). Because a zero concentration bottom boundary has no relevance for a hydrodynamic dispersion process, solutions of the CDE in a semi-infinite profile are more relevant to predict BTCs or to interpret BTCs from leaching experiments in short columns.
View larger version (27K):
[in this window]
[in a new window]
|
Fig. 7. (a, b, c) Simulated resident concentration, Cr, depth profiles and (d, e, f) flux concentration, Cf, breakthrough curves at the 2-m depth for different dispersivities, . (Black lines are analytical benchmarks: Eq. [19] for depth profiles and Eq. [21] for breakthrough curves; similarity in output may overlap in the graph.)
|
|
Figure 8
shows fronts of a nonlinearly sorbing solute that infiltrates in a homogeneous soil profile. Due to the convex shape of the Freundlich isotherm, the shape of the fronts remains constant with time and traveling waves are developed. The shapes of the simulated solute fronts coincide well with the asymptotic analytical solution, Eq. [55]. The small deviation of the numerically simulated front for t = 20 d indicates that the asymptotic regime was not yet reached at that time.
View larger version (17K):
[in this window]
[in a new window]
|
Fig. 8. (a) Simulated concentration depth profiles of a nonlinearly sorbing substance, which is continuously applied at the soil surface and (b) concentration depth profile against the transformed depth coordinate Eq. [52]. (Full black line is the analytical benchmark Eq. [55]; similarity in output may overlap in the graph.)
|
|
Figure 9
shows a transport simulation during a stop-flow experiment of a substance undergoing rate-limited sorption. In a stop-flow experiment, the flow is interrupted for a certain period allowing equilibration of the sorbed and dissolved concentrations. When equilibrium is reached, flow is restarted. The effect of rate-limited sorption is evident from the early arrival of the concentration peak followed by a long tailing. For equilibrium sorption, the peak arrival would be expected at approximately 20 pore volumes. However, the effect of the flow interruption on the effluent concentrations is small. Because of a small increase of the sorbed concentration with increasing distance from the column outlet, the effluent concentrations increase just after restarting the flow. The numerical simulations all agree well with the analytical solutions except for the WAVE model. WAVE predicted sorbed concentrations before flow interruption were higher than predicted with the analytical solution.
View larger version (18K):
[in this window]
[in a new window]
|
Fig. 9. (a) Depth profiles at 200 h after tracer pulse injection of concentration in the soil solution Cr (full line), and sorbed concentration, Sr (dashed line), (Black lines are analytical benchmarks Eq. [23] and [24].) (b) Breakthrough of flux concentrations before and after a flow interruption at 200 h or 20 pore volumes T. (Black symbols are analytical benchmark Eq. [23] before flow interruption and Eq. [25] after interruption.) Similarity in output may overlap in the graph.
|
|
Time series of resident concentrations in the macropore and matrix region of a dual-porosity medium are shown in Fig. 10
for two different solute exchange rates. Simulations by HYDRUS_1D and MARTHE coincided with the analytical solution very well. For the fast exchange rate (Fig. 10b), MACRO simulations also coincided with the analytical solution. MACRO's assumption of purely convective transport in the macropore region did not lead to large differences between the MACRO simulations and the analytical solutions, which were obtained for a convective–dispersive transport in the macropore region. Therefore, the dispersion of solutes in the macropore and micropore regions is predominantly determined by the fast solute exchange between both regions. When the exchange rate is smaller (Fig. 10b), the dispersion in the macropore region defines the peak concentration. Although MACRO assumes purely advective transport in the macropore region, the peak concentrations are underestimated by the MACRO model. This indicates that the numerical solution of the purely advective transport in the macropore region leads to some numerical dispersion. However, the tailing part of the BTC in the macropore region and the concentrations in the matrix region are not sensitive to dispersion processes in the macropore region.
View larger version (26K):
[in this window]
[in a new window]
|
Fig. 10. Time series of resident concentrations in the macropore (full line) and matrix region (dashed line) of a dual-velocity medium with (a) slow and (b) fast exchange of solutes between both flow domains. Time series at 100 cm downstream of the initial concentration pulse are shown. (Black symbols are analytical benchmark for the macropore Eq. [29] and for the matrix region Eq. [30]; similarity in output may overlap in the graph.)
|
|
Simulated solute fluxes from a 2-m-deep soil profile that were simulated with climatic boundary conditions are shown in Fig. 11 . Time series of simulated concentrations at the soil surface and at the 1.0-m depth are shown in Fig. 12
. We considered soil surface concentrations since these concentrations are relevant for modeling volatilization losses. Simulations by WAVE, which suffered from numerical dispersion (Fig. 7), did not coincide with results of the other models. WAVE simulated an earlier solute breakthrough, which can be explained by the artificially larger solute dispersion, but the WAVE model simulates much higher solute flux peaks, which seem to be in contradiction with the artificially larger solute dispersion. The solute fluxes (Fig. 11) and concentration time series (Fig. 12) simulated by HYDRUS_1D, SWAP, MACRO, and MARTHE are similar. The MACRO model, however, tends to simulate an earlier solute breakthrough, higher solute flux peaks (Fig. 11), and somewhat smaller concentration increases at the soil surface due to upward water flow during evaporation (Fig. 12). This is a consequence of the smaller evaporation fluxes and larger deep drainage that are simulated by MACRO (Fig. 6).
View larger version (17K):
[in this window]
[in a new window]
|
Fig. 11. Simulated solute fluxes at 2-m depth in (a) sandy, (b) loamy, and (c) clayey soil profiles for climatic boundary conditions.
|
|
View larger version (26K):
[in this window]
[in a new window]
|
Fig. 12. Time series of concentrations at the soil surface (full line) and at the 100-cm depth (dashed lines) simulated for climatic boundary conditions in (a) sandy, (b) loamy, and (c) clayey soil profiles.
|
|
|
CONCLUSIONS
|
|---|
Analytical benchmarks are helpful to evaluate numerical flow and transport models. Although the benchmarks only exist for simple and unrealistic boundary conditions, they help to interpret the differences among and accuracy of model simulations for realistic boundary scenarios. Testing of numerical models against analytical benchmarks requires the potential implementation of unrealistic initial and boundary conditions. This goes against building in safeguards against improper definition of boundary conditions or material properties or parameters by less experienced users.
Differences among water flow simulations by different models can be attributed to the implementation of the soil surface boundary condition in the models. Both the spatial discretization of the pressure head profile close to the soil surface and the methods of averaging the hydraulic conductivities in the first grid layer influence the numerical solutions. The impact of conductivity averaging in the first grid layer on the soil water balance under realistic boundary conditions is, however, small. Nonetheless, the differences among simulated fluxes at the soil surface may be more relevant when processes that depend on soil surface states (e.g., volatilization of pesticides) have to be described or when soil surface states (e.g., soil water content or soil surface temperature that are derived from remote sensing) are used to calibrate soil water flow models using inverse modeling (e.g., Jhorar et al., 2002). Because evaporation leads to large vertical pressure head gradients close to the soil surface, a fine vertical spatial discretization is required. This implies that for an explicit three-dimensional modeling of water and solute fluxes at larger scales, a much finer discretization in the vertical than in the horizontal direction should be used. Whether the Richards equation and CDE can be used to describe average flow and transport processes in very thin elements with a large horizontal extent is still a point of scientific debate.
For the simulation of solute transport, a correct simulation of the solute dispersion is essential. When spurious numerical dispersion occurs, high solute flu
TOP
ABSTRACT
INTRODUCTION
