A Transformed Pressure Head‐Based Approach to Solve Richards' Equation for Variably Saturated Soils

Abstract

This paper presents a new approach to solve Richards' equation. It introduces a nonlinear transformed pressure, Pt, as the dependent variable with the modified Picard method. The new approach was compared to, and contrasted with, two efficient existing methods: the o-based transformation method (Kirkland et at., 1992), and the h-based modified Picard method (Celia et at., 1990). A total of 12 different one-dimensional cases were considered (saturated, unsaturated, layered and uniform soil profiles, with pressure and flux type boundary conditions). The results show that the new method offers excellent CPU et•ciency and, unlike the h-based method, is numerically robust for all cases of variabty saturated, heterogeneous media, and first or second type boundary conditions. The method does not require difficult numerical coding, and its CPU efficiency is not affected by complicated heterogeneous and hysteretic media. The Pt transformation is easy to incorporate into existing h-based codes.where 0 and h are volumetric water content and soil water pressure head, respectively; K is the hydraulic conductivity, which is assumed constant for saturated soils but varies strongly with 0 or h in unsaturated soils; z denotes the vertical dimension; and t is time. To solve (1), the soil water retention function is introduced to eliminate one of the two dependent variables. This results in either an "h-based" or a "O-based" form of (1). Because of the highly nonlinear K-h and O-h relationships, analytical solutions to (1) are impossible, except for special cases. However, there exist many numerical methods.Numerical models to solve the h-based form of (1) seem to be the most common. They can be used for saturated and unsaturated soils, as well as for layered soils. However, these models suffer from poor mass balances for unsaturated soils, and from unacceptable time step limitations, or poor CPU efficiency especially for very dry initial conditions. Milly [1985] presented a mass-conserving solution procedure that uses a modified definition of the capacity term to force global mass balance. Allen and Murphy [1985, 1986] Bouloutas [1989] and Celia et al. [1990] greatly improved the performance of h-based models by using a mixed form of Richards' equation. Their numerical algorithms are denoted as "quasi-Newton" or "modified Picard," but probably should be called "h-based quasi-Newton" and "h-based modified Picard" because pressure head is still the dependent variable. The strategy used by Celia et al. [1990] is to evaluate the water content change over one time step directly from the change of the pressured head. This results in much improved mass balances. However, as shown by Kirkland et al. [1992], CPU efficiency is still a problem with h-based models for very dry initial conditions.

Gridded Surface Subsurface Hydrologic Analysis (GSSHA) User's Manual; Version 1.43 for Watershed Modeling System 6.1

“…With a head-based formulation, rainfall, evaporation and groundwater recharge can all be modeled without a change in variable. Known mass balance problems associated with solution of the head-based formula (Celia et al 1987;Ross 1990;Pan and Wierenga 1995;Celia et al 1990) have largely been eliminated by the development of new solution techniques (Kirkland et al 1992;Rathfelder and Abriola 1994;Pan and Wierenga 1995).…”

New insights into the differences between the dual node approach and the common node approach for coupling surface–subsurface flow

“…A relatively fine uniform vertical discretization also enables the simulation of sharp saturation fronts with the Richards equation (Pan and Wierenga, 1995;Ross, 1990). As such, the simulation results based on the finest vertical discretization can be taken as reference solutions that enable comparisons of the coupling approaches when a coarser vertical discretization is used.…”

Contribution to Numerical Modeling of Water Flow in Variably Saturated, Heterogeneous Porous Media

“…PAN and WIERENGA (1995) presented a method to simulate problems with convergence difficulty attributed to the presence of sharp wetting fronts. The pressure head variable h is then transformed into new dependent variables (3a) and (3b):…”

“…Where: κ is a constant (-0.05 < κ < -0.01 cm -1 ) independent of both the K(h) and C(h) relationship (PAN and WIERENGA, 1995); he is a free parameter, referred as the air entry value [L] (IPPISCH et al, 2006); is the transformed pressure head [L].…”