Vadose Zone Journal 2:193-200 (2003)
© 2003 Soil Science Society of America
Numerical Simulations of One-Dimensional Infiltration into Layered Soils with the Richards Equation Using Different Estimates of the Interlayer Conductivity
B. Brunonea,
M. Ferrantea,
N. Romano*,b and
A. Santinib
a Department of Civil and Environmental Engineering, University of Perugia, Perugia, Italy
b Department of Agricultural Engineering, University of Naples Federico II, Portici (Naples), Italy
* Corresponding author (nunzio.romano{at}unina.it)
Received 30 August 2002.
 |
ABSTRACT
|
|---|
Transient water flow processes in unsaturated soils are usually modeled using the Richards equation. This paper compares several numerical approximations to this equation for vertical infiltration in layered soil profiles. Three formulations of the governing flow equation (i.e., the h-based, 
-based, and mixed forms) are compared for the critical test problem of infiltration into a layered soil profile with initially relatively low soil water contents. An efficient, yet relatively simple weighting algorithm is employed that improves the estimation of the interlayer hydraulic conductivity in an h-based finite-difference formulation. Results highlight improvements in the mass conservative properties of this model, which is termed the H-IL model. A comparison is then performed between the finite-difference H-IL model and a finite-element model for infiltration toward a water table. The H-IL model was found to be computationally very efficient for this test problem. For both illustrative flow examples, the different numerical models were evaluated in terms of their ability to reproduce an exact solution developed by Srivastava and Yeh (1991) in their work on one-dimensional, transient infiltration toward the water table in homogeneous and layered soils.
Abbreviations: FD, finite difference • FE, finite element • H-GM model, same as FD h-based model H-IL, but with geometric means (Eq. [4]) for computing the interlayer conductivities • H-IL, FD h-based model with interlayer hydraulic conductivity computed by algorithm of Romano et al. (1998) (see Eq. [1] and [5]) • MIX-GM, FD mixed model with geometric mean conductivities for all interior nodes (see Eq. [3] and [4]) • THT-GM, FD -based model with geometric mean conductivities for all interior nodes (see Eq. [2] and [4])
|
INTRODUCTION
|
|---|
WATER FLOW IN THE VADOSE ZONE has been the subject of numerous analytical, numerical, and experimental studies (Kutílek and Nielsen, 1994). Many have studied the role that numerical solutions of the governing flow equation exert on quantitative simulations of practical problems, such as estimation of groundwater recharge or the transfer of water in the soil–plant–atmosphere system.
Soil water flow processes are commonly assumed to evolve in a rigid, isothermal porous medium, with flow occurring one-dimensionally in the vertical direction. Combining Darcy's equation and the principle of mass conservation yields the following equation:
 | [1] |
which is known as the Richards equation. In this equation, soil depth, z, and time, t, represent the independent variables, whereas the pressure head, h, is the dependent variable. The system parameters are the soil hydraulic conductivity function, K(h), and the specific soil water capacity, C(h), given by C(h) = d/dh and readily computed from the soil water retention function, (h), where is the volumetric water content. Equation [1] is usually referred to in the literature as the h-based formulation of the Richards equation and has led to valuable results when interpreting experimental data collected at both laboratory or small plot scales (van Dam and Feddes, 2000). The equation has also been employed in process-based distributed models to predict the hydrologic response of river basins (e.g., Paniconi and Wood, 1993).
Exact analytical solutions of Eq. [1] are available, but a basic limitation is that they are usually obtained for simplified flow regimes or for soils with hydraulic properties being realistic only within a relatively limited range of soil water pressure heads. To facilitate the derivation of analytical solutions, or to enhance the numerical tractability of the governing flow equation, the concept of soil water diffusivity D() = K() dh/d was introduced to produce the following diffusive-type flow equation (Klute, 1952):
 | [2] |
in which represents the dependent variable.
A third formulation of the governing unsaturated flow equation, the mixed form of the Richards equation, was used by Celia et al. (1990):
 | [3] |
in which and h are the dependent variables.
One feature that is common to Eq. [1], [2], and [3] is the difficulty of their solutions because of the strong nonlinearity of the functions C(h), K(h), and D(). Only for particular types of soils (e.g., linear soils) and specific boundary conditions can these equations be solved using analytical techniques. Milly (1988) showed that numerical evaluation of the capacity term C(h)
h/t in Eq. [1] can generate major non-mass-conserving solutions of the h-based formulation. Rathfelder and Abriola (1994) examined factors that may affect mass conservation for different numerical approximations of the unsaturated flow equation and argued that efficient mass-conservative solutions of the h-based Richards equation can be obtained by appropriately discretizing the capacity term. A major benefit of employing Eq. [2] is due to the fact that soil hydraulic functions generally show less nonlinearity when they are expressed vs. instead of h (Hills et al., 1989). The -based form of the Richards equation can be handled easier from a numerical point of view and appears very suitable when solving specific problems, such as infiltration into an initially dry soil. On the other hand, Eq. [1] and [3] have the advantage of being applicable to the flow region as a whole, including saturated and unsaturated flow, while also being very effective in the case of layered soil profiles, since h remains continuous across interfaces as opposed to .
Equation [1], and equivalently Eq. [2] or [3], is a parabolic partial differential equation that must be solved numerically for most practical unsaturated flow problems. Two methods are usually employed when solving the Richards equation: finite differences (FD) and finite elements (FE). Both methods yield similar systems of nonlinear algebraic equations.
|
Estimation of the Finite-Difference Internodal Hydraulic Conductivity
|
|---|
When FD methods are used, various strategies are available to handle the geometry of the flow domain. Specific arrangements are required when complex geometries are present. For one-dimensional solutions implementation of a single-structured grid is very popular, especially with the present availability of powerful and affordable computers that enable the use of relatively fine steps in both space and time (van Dam and Feddes, 2000). Once a uniform spatial grid has been set up, FD solutions are commonly obtained using a mesh-centered or a cell-centered grid, although other schemes are also available (Patankar, 1980; Mansell et al., 2002). For mesh-centered schemes, nodes are located at the nodal points of the grid. Such schemes efficiently treat point sinks or sources (e.g., pumping or recharging wells) and appear more suitable for Dirichlecht type (i.e., h or specified) boundary conditions (Gureghian et al., 1979). However, mesh-centered schemes are prone to mass balance errors. The importance of setting up mass-conservative methods and the fact that most practical subsurface flow problems evolve under imposed flux boundary conditions (i.e., Neumann type boundary conditions) has spurred the developments of cell-centered (or, block-centered) numerical schemes (Hills et al., 1989). Such schemes assume that the flow variables are evaluated at the center of computational cells, and that the cell boundaries coincide with external or internal domain boundaries (see Fig. 1). Using central FD and a cell-centered grid scheme, the flux is computed with second-order accuracy at the midpoint between nodes where the domain boundary is located.
View larger version (26K):
[in this window]
[in a new window]
|
Fig. 1. Computational grid near the interface between two soil layers, with schematic distributions of the pressure heads at a certain time during the transient flow.
|
|
Because of its nonlinearity, spatial derivatives of the Richards equation are often evaluated by using internodal hydraulic conductivities, Kin; that is, values of K are estimated from the known hydraulic conductivities at two neighboring nodes of the spatial grid system. Previous studies showed that the accuracy of FD models strongly depends on the weighting technique used for computing Kin (Haverkamp and Vauclin, 1979; Zaidel and Russo, 1992; Singh and Bhallamudi, 1998; Baker, 2000). The use of simple arithmetic means has been criticized often, partly because it can cause physical inconsistencies (Haverkamp and Vauclin, 1979; Warrick, 1991a; Baker, 2000). Geometric means as initially proposed by Haverkamp and Vauclin (1979) are now a widely used criterion for evaluating Kin = Ki+1/2 between two consecutive nodes i and i + 1 of the computational spatial grid:
 | [4] |
It is also important to note that improper treatment of the interfaces can affect the stability of a numerical method. Hence, the presence of layers in a soil profile, with related discontinuities in the constitutive soil hydraulic characteristics, requires careful evaluation of the hydraulic conductivity at or close to the interface. This applies to both mesh-centered and cell-centered schemes (Gureghian et al., 1979; Hills et al., 1989). For a layered soil profile, we will use in our numerical schemes the interlayer hydraulic conductivity, Kil. Some authors compute Kil by using the same averaging criterion as defined for a uniform soil profile; namely, they do not distinguish between internodal and interlayer conductivities (Vauclin et al., 1979; Hills et al., 1989; Berg, 1999).
Romano et al. (1998) calculated the internodal hydraulic conductivity for neighboring nodes in the same soil layer using the geometric mean criterion of Eq. [4]. However, they developed a local refinement technique to compute Kil that accounts for discontinuities in the first derivative of h(z,t). The algorithm for computing Kil can be schematically outlined as follows (see Fig. 1):
- Adopting a FD, Crank–Nicolson-type model and invoking a cell-centered scheme such that two nodes are placed symmetrically at equal distances from the boundaries or a horizontal interface between two soil layers—at each time level of the flow process the pressure head, h, and the flux, q, are then continuous across the interface.
- Fictitious values of the pressure head are defined at two adjacent nodes by extrapolating the pressure head profiles h(z,t) of the upper and lower soil layers, respectively.
- Conservation of mass (i.e., the flux entering the interface must be equal to the flux exiting the interface) yields the following expression for the interlayer hydraulic conductivity:
 | [5] |
in which the K terms do not represent any nodal hydraulic conductivity values. Rather, K- and K+ must be calculated as geometric averages of the conductivity values corresponding to both nodal and fictitious values of the pressure head. We refer to a paper by Romano et al. (1998) for additional details. It is worth mentioning that this procedure is simple, computationally efficient, and accurate in minimizing mass balance errors (used here as a global criterion) and discrepancies between numerical vs. analytical solutions for two-layer flow domains (used here as a local criterion). Recently, the algorithm developed by Romano et al. (1998) was employed in a finite-difference approximation of the -based form of the Richards equation (Schaudt and Morrill, 2002).
Unsaturated flow generally evolves in soil profiles with contrasting physical and hydraulic properties. However, only a relatively few studies have performed numerical tests for layered soils and compared results against analytical solutions of the governing flow equation. We performed a systematic investigation to simulate transient flow in layered soil profiles using numerical and analytical techniques. Besides providing additional verification of the algorithm proposed by Romano et al. (1998) to compute interlayer hydraulic conductivities, the specific objectives of this paper are (i) to compare finite difference results of the three formulations of the Richards equation for the critical case of vertical infiltration into initially dry soils and (ii) to answer the question of how, and to what extent, a more accurate evaluation of interlayer hydraulic conductivities can alleviate local and global mass balance errors in the h-based formulation. The performance of the finite difference model of Romano et al. (1998) is also compared with a finite element model for infiltration in a layered profile in the presence of a water table. All of the numerical results are compared against an analytical model developed by Srivastava and Yeh (1991) to predict water content and pressure head distributions for vertical infiltration in layered soils.
|
Simulation Conditions and Evaluation Strategies
|
|---|
The h-based numerical solutions with Eq. [5] for the interlayer conductivities, as proposed by Romano et al. (1998), is referred to hereafter as the H-IL model. This model adopts an implicit, Crank–Nicolson–type FD scheme with explicit linearization of the C and K functions. The same FD h-based model, but with geometric means (Eq. [4]) for computing the interlayer conductivities, will be called the H-GM model. Assuming only the weighting criterion of Eq. [4] for all interior nodes of the flow domain, the FD solutions based on either the -based, Eq. [2], or mixed form, Eq. [3], are referred to as the THT-GM or the MIX-GM model, respectively. For maintaining consistency in the comparisons and highlighting the effects of different ways of computing the parameter Kil, we note that the -based and mixed-based equations were implemented in similar fashions as the h-based equation, except for changes in the relevant dependent variable where appropriate. Additionally, we verified that the time step size 
t was nearly identical for all of the simulation runs.
In this study all simulations were compared directly with the exact analytical solutions developed by Srivastava and Yeh (1991), rather than evaluating results in terms of mass balances, or their self consistency when using a much finer computational grid system. This would make our analysis more effective since the distribution of the mass within the flow domain may not be correct even with an accurate global mass balance. The analytical solution provides distributions of and h in a layered soil profile that is initially at steady state, but subjected to an abrupt change in the flux density at the upper boundary. The lower boundary of the flow domain was subjected to a constant pressure head. Soil hydraulic properties were described by the following exponential-type relationships (Srivastava and Yeh, 1991):
 | [6a] |
 | [6b] |
where s (cm3 cm-3) and Ks (cm h-1) are the soil water content and hydraulic conductivity at saturation, respectively, whereas r represents the residual water content when h tends to -
. Parameter 
(cm-1) is related to the particle-size distribution of the porous material and has the same value in both the soil water retention (Eq. [6a]) and hydraulic conductivity (Eq. [6b]) functions.
|
MODEL COMPARISONS
|
|---|
Infiltration into an Initially Dry Layered Soil Profile: Comparisons among Finite-Differences Models
For this example, we assume flow conditions similar to those reported in a study by Hills et al. (1989). These authors analyzed infiltration into an initially dry, layered soil profile made up of a sequence of layers of Berino loamy fine sand (fine-loamy, mixed, superactive, thermic Ustic Calciargid) and Glendale clay loam (fine-silty, mixed, superactive, calcareous, thermic Typic Torrifluvent), with the Berino soil as the uppermost layer. The soil water retention and hydraulic conductivity properties of these soils were described by van Genuchten's parametric relations (Kool and van Genuchten, 1991). A uniform pressure head distribution, h(z,0) = h0, was imposed as initial condition, with h0 for the different cases ranging between -1000 and -10 000 cm.
To assess the effects of initial and boundary conditions on the results, we define the following dimensionless parameters:
 | [7] |
where W0 is the total volume of water in a soil at t = 0 (i.e., when h = h0) and Wmax is the volume of water at complete saturation, and
 | [8] |
where K0 is the hydraulic conductivity for h = h0 and Ks is the value of K at saturation. For each of the h0 values imposed by Hills et al. (1989), Table 1 reports the values of the parameters 
and 
for the two types of soils used in their example.
The current problem considers a vertical infiltration into a 100-cm-deep (L = 100 cm) soil profile. The profile consists of two layers, with the upper soil layer having a thickness of L1 = 20 cm and the lower soil layer having a thickness of L2 = 80 cm. Values of the soil hydraulic parameters in Eq. [6a] and [6b] are Ks1 = 1.0 cm h-1, Ks2 = 10.0 cm h-1, s1 = s2 = 0.40, r1 = r2 = 0.06, and 1 = 2 = 0.10 cm-1, where the subscripts 1 and 2 are for the upper and lower horizons, respectively. We are thus considering a scenario in which the less conductive soil overlies the more conductive one. The lower boundary condition is subjected to a pressure head of -100 cm. The flow equation is solved by assuming as initial condition steady-state profiles corresponding to a prescribed rain intensity qini of 4.54 x 10-4 cm h-1. The imposed boundary flux at the soil surface is subsequently changed abruptly to a constant value, qfin, of 0.95 cm h-1. The nodal spacing of the computational grid was set at z = 0.25 cm, whereas the time step was adjusted such as the maximum change in water content between two consecutive times {i.e., max = max[(z, t + t) - (z, t)]} is <1.0%.
Application of the H-IL model to the above example produced a global value of of 0.15 for the entire soil profile, and values of that ranged from a maximum of 4.05 x 10-4 at the soil surface to a relatively constant value of approximately 4.54 x 10-6 in the lower soil layer. Hence flow occurs in a soil profile that initially has very low water contents.
Figure 2 illustrates pressure head, h(z,t) (Fig. 2a), and water content, (z,t) (Fig. 2b), profiles at selected times during the infiltration process as computed with the Srivastava and Yeh (1991) analytical model. The distributions at time t = 0 represent the imposed initial condition. Notice that the water contents are continuous across the interface, which is consistent with the analytical solution and the fact that we used the same values for s, r, and for both soil layers. Assuming the analytical solution to be the "true" solution of the flow field, comparisons between the various finite difference solutions are performed by calculating the local relative error in the water content as follows:
 | [9] |
where num(z,t) is numerically calculated water content at depth z and time t and true(z,t) the water content computed with the analytical solution of Srivastava and Yeh (1991).
View larger version (23K):
[in this window]
[in a new window]
|
Fig. 2. (a) Pressure head and (b) water content distributions at selected times as computed with the analytical model of Srivastava and Yeh (1991) for a two-layer soil profile assuming an infiltration scenario with qini = 4.54 x 10-4 cm h-1 and qfin = 0.95 cm h-1. The horizontal bold line indicates the location of the interface between the two soil layers.
|
|
Figure 3 shows the local error 
(z,t) during the transient flow event for the four numerical models tested in this paper. To obtain a better visual comparison, results for only the upper 40 cm of the soil profile are depicted. Inspection of Fig. 3a and 3b confirms that the H-IL solution performs better than the H-GM model, which utilizes the geometric mean for all internodal conductivities of the flow domain. Distributions of vs. z appear similar for both the H-IL and H-GM models only at the beginning of the infiltration process when water contents are still relatively low. For t > 2 h, the error at the layer interface for the H-GM model becomes large and affects also the local error elsewhere in the domain. These features will be discussed in more detail below. As expected, THT-GM (Fig. 3c) and MIX-GM (Fig. 3d) perform well at the beginning of the flow process, but also produce relatively large errors at the layer interface. The errors at the layer interface for THT-GM and MIX-GM eventually become even larger than those for the H-GM model and reach a maximum of about -3% at t = 40 h, as compared with -2% for the H-GM model.
The results depicted in Fig. 3 can also be compared in terms of the relative mass balance error, W(t), as follows:
 | [10] |
in which:
 | [11] |
are the volumes of water in the upper 40 cm of the soil as computed numerically and with the analytical solution, respectively. For the four numerical models considered, Fig. 4 shows mass balance errors, W(t), as a function of time, t. Notice that at early times the mass balance errors for H-IL and H-GM models are nearly identical, showing a peak value of about 0.8% at t = 2.0 h. This can be explained by the fact that at early times the wetting front has not yet reached the interface between the soil layers (i.e., before the differences between H-IL and H-GM are expected to increase). We note here also that the geometric mean criterion generally does not describe infiltration into relatively dry soils very well. The geometric mean (Eq. [4]) is weighted toward the lower value of K and, therefore, often underestimates water fluxes under dry soil conditions. For times greater than approximately 3 h, the H-IL and H-GM curves start diverging. Mass balance errors for the H-GM model rapidly reach a nearly constant value of about -0.22% for t > 15 h. When t exceeds 15 h, mass balance errors for the H-IL model are close to 0.05%, which is negligible from a practical point of view. Relative mass balance errors for the THT-GM and MIX-GM models are nearly identical and show typical oscillations at early times of the simulation run. For t > 15 h, mass balance errors for these two numerical models stabilize at values approximately close to -0.40%. Notice that for times greater than about 20 h flow has essentially reached a new steady-state condition, referred to hereafter as flux density qfin. This is reflected by the mass balance errors stabilizing to a constant value, W, for each particular numerical model.
View larger version (33K):
[in this window]
[in a new window]
|
Fig. 4. Relative mass balance error, W, as a function of time, t, obtained with four finite difference models. The H-IL model is the h-based numerical model with Eq. [5] for computing the interlayer conductivity. The H-GM, THT-GM, and MIX-GM are the h-based, -based, and mixed formulations of the Richards equation, respectively, all using the geometric mean (Eq. [4]) as averaging criterion for computing the interlayer conductivity.
|
|
The present example is for a layered soil in which a less conductive zone (Ks1 = 1.0 cm h-1) overlies a more conductive layer (Ks2 = 10.0 cm h-1). Romano et al. (1998) previously showed that for this case the H-IL numerical model was less efficient computationally as compared with other schemes.
Infiltration toward a Water Table: Comparison between Finite-Differences and Finite-Elements Models
In this section, the performance of the H-IL model is compared with the results obtained with the HYDRUS finite element model of Kool and van Genuchten (1991). HYDRUS employs a fully implicit, Galerkin-type linear finite element solution of the Richards equation and assumes that the soil hydraulic properties, (h) and K(), are described by the parametric functions of van Genuchten (Kool and van Genuchten, 1991). Consequently, the original HYDRUS code was modified to handle the constitutive relationships given by Eq. [6].
This example is for infiltration in a layered soil profile with a constant groundwater level at its base. The flow domain has a total length, L, of 200 cm and consists of two soil layers of equal thickness (100 cm each). Soil water retention (Eq. [6a]) and hydraulic conductivity (Eq. [6b]) parameters for the two soil layers were: s1 = s2 = 0.40, r1 = r2 = 0.06, 1 = 2 = 0.10 cm-1, but Ks1 = 10.0 cm h-1 and Ks2 = 1.0 cm h-1, where the subscripts 1 and 2 are for the upper and lower horizons, respectively. Hence, the more conductive soil overlies the less conductive one.
The initial condition is a steady-state pressure head distribution corresponding to a constant flux qini equal to 0.010 cm h-1, while the surface flux at t = 0 is abruptly changed to qfin = 0.90 cm h-1. Assuming the presence of a water table at z = 200 cm, the lower boundary condition equals h(L,t) = 0. Figure 5 shows calculated pressure head distributions at selected times during the infiltration event (including the initial steady-state profile) as computed with the analytical model. The new steady-state condition corresponding to the flux qfin is reached at approximately 50 h after the beginning of the infiltration process.
View larger version (29K):
[in this window]
[in a new window]
|
Fig. 5. Pressure head distributions, h(z), at selected time, t, as computed with the analytical model of Srivastava and Yeh (1991) for a two-layer soil profile assuming an infiltration scenario with qini = 0.010 cm h-1 and qfin = 0.90 cm h-1. The horizontal bold line represents the location of the interface between the two soil layers.
|
|
Comparisons between the reference analytical solution and the FD H-IL model and FE HYDRUS numerical results are presented in terms of the relative local error in the pressure head; that is,
 | [12] |
where hnum(z,t) is the pressure head at soil depth z and time t computed numerically and htrue(z,t) represents as before the analytical solution of Srivastava and Yeh (1991).
Results obtained with the FD H-IL model and the analytical solutions are virtually indistinguishable in Fig. 6. This figure thus shows only the distributions with depth and time of the local error h for the FE HYDRUS model. Notice that close to the layer interface HYDRUS overestimates the analytical values in the upper soil layer, while underestimates the analytical results in the lower soil layer. Moreover, h was found to increase significantly with time in the upper soil layer.
View larger version (22K):
[in this window]
[in a new window]
|
Fig. 6. Distributions of the local error in h, h, with respect to Srivastava and Yeh's analytical solution as obtained with the finite-elements HYDRUS model.
|
|
The results above demonstrate the importance of selecting a suitable numerical technique to solve specific problems and confirm the excellent performance of the H-IL model when simulating water flow in layered soil profiles. We present here one additional example of the H-IL solution showing the potentially important effects of soil layering on groundwater recharge from the vadose zone. Computations for this last example use the same infiltration scenario as above, but with L1
L2.
Figure 7 depicts relative flux rates q/qfin as a function of relative time t/t* for three different arrangements of a two-layer vertical soil profile as identified by the ratio 
= L1/L2. The relative flux rates were computed at both the layer interface (open symbols) and the groundwater table (closed symbols). The groundwater table itself was kept at the same level (200 cm). Following Warrick (1991b), time t* is defined as:
 | [13] |
where L = L1 + L2. The other terms in Eq. [13] pertain to an equivalent uniform soil profile and are calculated as follows:
 | [14] |
The q/qfin values computed with the H-IL numerical model are in a very good agreement with those obtained with the analytical model by Srivastava and Yeh (1991), and coalesce essentially to the same curves. Notice that for the different cases considered, the relative flux, q/qfin, vs. relative time, t/t*, follows an S-shaped curve, with some tailing evident especially for case = 1/3 and for the fluxes computed at the layer interfaces. This is partly due to the chosen flow domain where a more permeable soil (Ks1 = 10 cm h-1) overlies a less permeable soil (Ks2 = 1 cm h-1). At q/qfin = 1, the curves pertaining to fluxes at the layer interface (open symbols) virtually overlap when t/t* exceeds a value of approximately 4. It is interesting to note that when the layer interface is relatively far from the water table (cases = 1/3 and = 3), the system conveys all water coming from the soil surface, with a delay of similar extent. Conversely, for = 9 the open and closed triangles are fairly close together, showing that in this case the rain drains rapidly toward the groundwater, since the more conductive upper layer now practically governs the flow process.
|
SUMMARY AND CONCLUSIONS
|
|---|
This investigation has provided further insights into the performance of the finite difference H-IL numerical model developed by Romano et al. (1998) for solution of the Richards equation. This model incorporates a convenient algorithm for computing the equivalent hydraulic conductivity at the interface between soil layers having different hydraulic properties. The model was evaluated quantitatively by comparison of analytical and numerical results for a number of simulation scenarios involving infiltration in layered soil profiles.
The following conclusions can be drawn from our numerical experiments:
- While limited to a certain class of soils, validation against the analytical solution enabled rigorous testing of the H-IL model and confirmed earlier statements about its accuracy and computational efficiency (Romano et al., 1998). Results in this paper confirm that computing the interlayer hydraulic conductivity appropriately is an important prerequisite for obtaining reliable unsaturated flow predictions.
- The proposed approach was verified for highly nonlinear flow conditions involving infiltration into an initially very dry layered soil profile. Although for this example some minor mass balance errors occurred at early times, the FD h-based Richards equation with the interlayer conductivity algorithm of Romano et al. (1998) provided overall robust solutions.
- The excellent agreement between our numerical results and the analytical solutions for the cases considered, together with comparisons among alternative FD and FE approximations of the Richards equation, strengthened our confidence in using the Romano et al. (1998) algorithm for Kil, not only for the h-based form, but also in other modeling approaches, such as the scheme developed by Schaudt and Morrill (2002). Inclusion of the algorithm in multidimensional variably saturated flow models may well permit the adoption of relatively coarse spatial grids without compromising accuracy and stability of the numerical results (Romano et al., 1998).
Testing the different formulations of the Richards equations against known analytical results shows the effectiveness of the mixed formulation when simulating relatively high pressure head gradients in the flow domain. However, relatively low local and global errors may be obtained with the h-form only if interlayer hydraulic conductivities are computed using a specifically designed weighting scheme. Coupling the algorithm of Romano et al. (1998) with the numerical scheme of Rathfelder and Abriola (1994) could make the h-based formulation of the unsaturated flow equation a very attractive model.
Although more accurate FD results employing more classic weighting schemes for Kil, such as the arithmetic or geometric averages, might be obtained by adopting smaller nodal spacings z, we believe that such an approach also has theoretical and philosophical drawbacks. For example, for certain layered or heterogeneous systems, as well as for certain types of highly structured soils, values for z should account for the support of the measuring sensors and the fact that the governing equations are averaged over a suitable representative volume.
|
ACKNOWLEDGMENTS
|
|---|
This work was supported in part by the University of Perugia, and by the Italian Ministry of Scientific Research and the University of Naples Federico II as part of the project "New technologies for the qualitative and quantitative monitoring, management and control of groundwater systems". We also appreciate the support provided by the GNDCI special project of the National Research Council of Italy on "Prediction and prevention of extreme hydrologic events and their control—Monitoring and modeling hydrologic processes in soil". We would like to thank T.-C. Jim Yeh, University of Arizona, Tucson, USA, for helpful comments on using the analytical model. The HYDRUS code was kindly provided by Rien van Genuchten of the USDA Salinity Lab., Riverside, USA.
|
REFERENCES
|
|---|
- Baker, D.L. 2000. A Darcian integral approximation to interblock hydraulic conductivity means in vertical infiltration. Comput. Geosci. 26:581–590.
- Berg, P. 1999. Long-term simulation of water movement in soils using mass-conserving procedures. Adv. Water Resour. 22:419–430.
- Celia, M.A., E.T. Bouloutas, and R.L. Zarba. 1990. A general mass-conservative numerical solution for the unsaturated flow equation. Water Resour. Res. 26:1483–1496.
- Gureghian, A.B., D.S. Ward, and R.W. Cleary. 1979. Simultaneously transport of water and reacting solutes through multilayered soils under transient unsaturated flow conditions. J. Hydrol. (Amsterdam) 41:253–278.
- Haverkamp, R., and M. Vauclin. 1979. A note on estimating finite difference interblock hydraulic conductivity values for transient unsaturated flow problems. Water Resour. Res. 15:181–187.
- Hills, R.G., I. Porro, D.B. Hudson, and P.J. Wierenga. 1989. Modeling one-dimensional infiltration into very dry soils. 1. Model development and evaluation. Water Resour. Res. 25:1259–1269.
- Klute, A. 1952. A numerical method for solving the flow equation for water in unsaturated material. Soil Sci. 73:105–116.
- Kool, J.B., and M.Th. van Genuchten. 1991. HYDRUS—One-dimensional variably saturated flow and transport model, including hysteresis and root water uptake. U.S. Salinity Laboratory, USDA, ARS, Riverside, CA.
- Kutílek, M., and D.R. Nielsen. 1994. Soil hydrology. GeoEcology Textbook. Catena Verlag, Cremlingen-Destedt, Germany.
- Mansell, R.S., L. Ma, L.R. Ahuja, and S.A. Bloom. 2002. Adaptive grid refinement in numerical models for water flow and chemical transport in soil: A review. Available at www.vadosezonejournal.org. Vadose Zone J. 1:222–238.[Abstract/Free Full Text]
- Milly, P.C.D. 1988. Advances in modeling of water in the unsaturated zone. Transp. Porous Media 3:491–514.
- Paniconi, C., and E.F. Wood. 1993. A detailed model for simulation of catchment scale subsurface hydrologic processes. Water Resour. Res. 29:1601–1620.
- Patankar, S.V. 1980. Numerical heat transfer and fluid flow. Hemisphere Publishing Corp., Washington, DC.
- Rathfelder, K., and L.M. Abriola. 1994. Mass conservative numerical solutions of the head-based Richards equation. Water Resour. Res. 30:2579–2586.
- Romano, N., B. Brunone, and A. Santini. 1998. Numerical analysis of one-dimensional unsaturated flow in layered soils. Adv. Water Resour. 21:315–324.
- Schaudt, K.J., and J.C. Morrill. 2002. On the treatment of heterogeneous interfaces in soil. Geophys. Res. Lett. 29:DOI: 10.1029/2001GL014376.
- Singh, V., and S.M. Bhallamudi. 1998. Conjunctive surface-subsurface modeling of overland flow. Adv. Water Resour. 21:567–579.
- Srivastava, R., and T.-C.J. Yeh. 1991. Analytical solution for one-dimensional, transient infiltration toward the water table in homogeneous and layered soils. Water Resour. Res. 27:753–762.
- Van Dam, J.C., and R. Feddes. 2000. Numerical simulation of infiltration, evaporation and shallow groundwater levels with the Richards equation. J. Hydrol. (Amsterdam) 233:72–85.
- Vauclin, M., R. Haverkamp, and G. Vachaud. 1979. Résolution numérique d'une équation de diffusion non linéaire. Presses Universitaires de Grenoble, Grenoble, France.
- Warrick, A.W. 1991a. Numerical approximation of Darcian flow through unsaturated soil. Water Resour. Res. 27:1215–1222.
- Warrick, A.W. 1991b. Traveling front approximation for infiltration into stratified soils. J. Hydrol. (Amsterdam) 128:213–222.
- Zaidel, J., and D. Russo. 1992. Estimation of finite difference interblock conductivities for simulation of infiltration into initially dry soils. Water Resour. Res. 28:2285–2295.
This article has been cited by other articles:
|
|
|
|
 
B. Belfort and F. Lehmann
Comparison of Equivalent Conductivities for Numerical Simulation of One-Dimensional Unsaturated Flow
Vadose Zone J.,
November 11, 2005;
4(4):
1191 - 1200.
[Abstract]
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
