Published online 23 January 2008
Published in Vadose Zone J 7:112-114 (2008)
DOI: 10.2136/vzj2006.0182
© 2008 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
TECHNICAL NOTES
Estimation of a Critical Spatial Discretization Limit for Solving Richards' Equation at Large Scales
H.-J. Vogela,* and
O. Ippischb
a UFZ Helmholtz Center for Environmental Research, Dep. of Soil Physics, Theodor-Lieser-Strasse 4, Halle 06120, Germany
b Institute for Parallel and Distributed Systems, Univ. of Stuttgart, Universitätsstrasse 38, 70569 Stuttgart, Germany
* Corresponding author (hans-joerg.vogel{at}ufz.de).
All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
Received 19 December 2006.
 |
ABSTRACT
|
|---|
Water dynamics in soil at spatial scales larger than the representative elementary volume (REV) of the porous structure are typically described by Richards' equation, which relates the flux law of Buckingham–Darcy to the mass balance of soil water. It is based on the soil water retention characteristics and the hydraulic conductivity function as constitutive material properties. In hydrological modeling, Richards' equation is also used at large scales up to hundreds of meters. Increasing the scale is typically accompanied by increasing the spatial discretization scale for the numerical solution of the problem. However, due to the underlying assumption of local equilibrium between water content and water potential, there is an upper limit of spatial discretization above which the solution is expected to be biased. We present a simple approach to estimate this limit, which depends on the shape of the soil hydraulic functions and the local gradient of total water potential. It is in the range between millimeters and decimeters.
Abbreviations: REV, representative elementary volume
|
INTRODUCTION
|
|---|
Application of Richards' equation at large scales is complicated by heterogeneity and nonlinearity, which prohibit upscaling from smaller to larger scales through simple averaging, as shown by Hopmans et al. (2002) and Beven (2001). The heterogeneity of the soil material can be taken into account in distributed modeling of water dynamics at larger scales using methods like pedotransfer functions (Wagner et al., 2001) or geostatistical tools (Yeh and Zhang, 1996). However, because of computational limitations, an increase of the model area is typically accompanied by a coarsening of the spatial discretization used in the numerical solution of Richards' equation. Downer and Ogden (2003) showed that this may actually lead to problems with the convergence of the solver and that the spatial discretization must be in the range of 1 cm or less in the case of infiltration fronts close to the surface due to numerical stability issues. Similar results were found by van Dam and Feddes (2000). Ross (1990) showed that the spatial discretization is especially critical in sandy material and proposed a dynamic discretization scheme at the front, which increased the efficiency considerably.
The aim of this note is to recall that an upper limit of the spatial discretization for solving Richards' equation exists. This is not relevant for stationary flow fields, but it is especially true for infiltration and drainage processes, which are considered here. If the spatial discretization is above a critical limit, this influences not only the convergence of the solver, which may be improved by more elaborate numerical techniques, but also the accuracy of the solution. We provide a simple estimation of the critical limit of spatial discretization for infiltration or drainage processes into homogeneous soil at a given initial water potential, 
m. In this case, the critical limit depends on the shape of the hydraulic functions and the steepness of the local gradient of total water potential.
|
Theory
|
|---|
If the water retention characteristic 
(m), relating water content and matric potential m, and the hydraulic conductivity function K() are known, water flow can be calculated from Richards' equation. The soil hydraulic functions are usually given by a set of parameters for one of the common parametrizations, such as the van Genuchten model (van Genuchten, 1980) or the Brooks–Corey model (Brooks and Corey, 1966). Since analytical solutions of Richards' equation are scarce and limited to very special cases, the partial differential equation is usually solved numerically based on discretizations in space and time. While the time step is often adapted automatically, an adaptive discretization in space is less common, and the spatial discretization is often fixed at a certain grid resolution, in the following referred to as discretization scale 
(L) describing the linear extension of a grid cell. Depending on the discretization and, if applicable, the type of basis functions used, the effects of a fixed differ. For higher-order schemes, the size of necessary to obtain a solution with reasonable accuracy will usually be larger but still limited.
In finite-difference schemes or finite-element and finite-volume schemes with linear basis functions, the hydraulic state variables or m are implicitly assumed to be piecewise constant or piecewise linear functions on each spatial element. Evidently, for large and high nonlinearities, this assumption is no longer valid. Across dynamic fronts of infiltration or drainage, the state variables are actually propagated across the discretization scale in a nonlinear manner. This effect increases with the amount of water required to equilibrate a potential gradient (i.e., with low values for dm/d). We want to derive an upper limit for based on the idea that for the validity of the linearity assumption, the hydraulic potential has to be close to local equilibrium. To estimate the critical limit for local equilibrium, the characteristic time required for equilibration can be compared to the time scale of convective water flow. If the latter is larger than the former, the soil hydraulic functions can be considered to be close to local equilibrium.
The soil water diffusivity
 | [1] |
describes the diffusive propagation of water due to capillarity. With respect to the discretization scale , an estimation of the characteristic time tdiff required for capillary equilibration can be obtained through
 | [2] |
In contrast, the characteristic time for convective water flow tflow with respect to the distance is given by
 | [3] |
where w = m + g is the total water potential including the component of gravity g.
From Eq. [2] and [3], it is clear that the time required for diffusive equilibration is proportional to 2, while the characteristic time for water flow is proportional to . Using Eq. [1] and [2], we obtain
 | [4] |
and
 | [5] |
The critical discretization scale can be estimated by equating Eq. [4] and Eq. [5], which yields
 | [6] |
An intuitive interpretation of this relation is that with increasing steepness of (m), the resolution of the spatial discretization has to increase. The same is true for an increasing gradient of total water potential. This is consistent with the experimental and numerical results of Downer and Ogden (2003), van Dam and Feddes (2000), and Ross (1990).
|
Results
|
|---|
The estimated minimum discretization length crit is shown in Fig. 1
as a function of matric potential and for different gradients of total water potential for a typical sand and a typical loam. A van Genuchten–Mualem model (van Genuchten, 1980) was used to describe the water retention characteristic and the hydraulic conductivity function. The van Genuchten parameters 
and n, together with the saturated water content s, saturated hydraulic conductivity Ks, and tortuosity 
, are given in Table 1
. The residual saturation r was assumed to be zero.
View larger version (15K):
[in this window]
[in a new window]
|
FIG. 1. Critical discretization length for infiltration as a function of matric potential for a typical sand (left) and a typical loam (right), with three different gradients of total water potential.
|
|
Obviously, in sand the critical spatial discretization length at unit gradient has a minimum of about 5 cm, which decreases quickly to a value below 1 cm as the gradient of total water potential becomes steeper. The location of the minimum corresponds to the steepest point of the water retention characteristic where 
m/ has its minimum. The behavior of the loam is similar. However, as the water retention characteristic is less steep, the critical spatial discretization length is larger than in the sand.
Generally, for each discretization scheme used in a numerical solution of Richards' equation, there is a minimal resolution that is necessary to obtain a solution with reasonable accuracy. This minimal resolution can be realized either by adaptive or by uniform refinement of the grid used for the computations. The approach presented here provides an estimate for the upper limit of this resolution obtained from the soil hydraulic functions and the maximal gradient of total water potential across infiltration or drainage fronts. It should be noted that this limit is thought to be applicable for perfectly homogeneous materials where small-scale heterogeneities can be lumped into effective hydraulic properties in the sense of an REV. In cases where the heterogeneities are in the size range of the discretization scale or larger, they may dictate a much finer resolution than prescribed by Eq. [6].
|
REFERENCES
|
|---|
- Beven, K. 2001. How far can we go in distributed hydrological modelling? Hydrol. Earth Syst. Sci.
5
:1–12.
- Brooks, R.H., and A.T. Corey. 1966. Properties of porous media affecting fluid flow. J. Irrig. Drain. Div. Am. Soc. Civ. Eng.
92
:61–88.
- Downer, C.W., and F.L. Ogden. 2003. Appropriate vertical discretization of Richards' equation for two-dimensional watershed-scale modelling. Hydrol. Earth Syst. Sci.
18
:1–22.
- Hopmans, J.W., J.
im
nek, N. Romano, and W. Durner. 2002. Simultaneous determination of water transmission and retention properties: Inverse methods. p. 963–1008. In J.H. Dane and G.C. Topp (ed.) Methods of soil analysis: Part 4. SSSA, Madison, WI. - Ross, P.J. 1990. Efficient numerical methods for infiltration using Richards' equation. Water Resour. Res.
26
:279–290.[CrossRef]
- van Dam, J.C., and R.A. Feddes. 2000. Numerical simulation of infiltration, evaporation and shallow groundwater levels with the Richards equation. J. Hydrol.
233
:72–85.[CrossRef]
- van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J.
44
:892–898.[Abstract/Free Full Text]
- Wagner, B., V.R. Tarnawski, V. Hennings, U. Müller, G. Wessolek, and R. Plagge. 2001. Evaluation of pedo-transfer functions for unsaturated soil hydraulic conductivity using an independent data set. Geoderma
102
:275–297.[CrossRef][Web of Science]
- Yeh, T.-C., and J. Zhang. 1996. A geostatistical inverse method for variably saturated flow in the vadose zone. Water Resour. Res.
32
:2757–2766.[CrossRef]
TOP
ABSTRACT
INTRODUCTION
used to estimate the critical discretization scale in