Infiltration → 1D Richards' Equation (3-layer) Method

The input variables used for modeling infiltration and unsaturated vertical flow with the 1D Richard's equation are defined as follows:

Ks = saturated hydraulic conductivity [m / s]
Ki = initial hydraulic conductivity [m / s] (typically much less than Ks)
θs = soil water content at ψ = 0 [unitless] (often set to the soil porosity, φ)
θi = initial soil water content [unitless]
θr = residual soil water content [unitless] (must be < θi)
ψB = bubbling pressure head [meters] (also called air-entry pressure, ψae)
ψA = pressure head offset parameter [meters]
λ = pore-size distribution parameter [unitless] (alt. notation = 1/b )
η = 2 + (3 * λ) [unitless] (see Notes)
c = transitional Brooks-Corey curvature parameter [unitless] (see Notes)
dznodes = vertical distance between nodes [meters]
nnodes = number of subsurface vertical nodes

For each variable, you may choose from the droplist of data types. For the "Scalar" data type, enter a numeric value with the units indicated in the dialog. For the other data types, enter a filename. Values in files must also use the indicated units.

Single grids and grid sequences are assumed to be stored as RTG and RTS files, respectively. Time series are assumed to be stored as text files, with one value per line. For a time series or grid sequence, the time between values must coincide with the timestep provided.

Note that the names of Greek letters are spelled out in the dialog, so theta = θ, psi = ψ and lambda = λ.

The current version allows up to 3 different soil layers, each with its own soil properties and number of vertical nodes. The soil layer number is indicated at the top of the dialog and you use the Next and Back buttons to view settings for the various layers.


Equations Used by the 1D Richards' Equation Method

v = K * (1 - ψz) = Darcy's Law for vertical flow rate [m / s]
vz = J - θt = conservation of mass,   with source/sink term J
Θe = (θ - θr) / (θs - θr) = effective saturation or scaled water content [unitless]
θr = θs ( |ψB| / 10000)λ = residual water content [unitless]
K = Ks * Θeη/λ = hydraulic conductivity [m / s] (see Notes below)
ψ = ψBe-c/λ - 1]1/c - ψA    = pressure head [meters] (see Notes below)

Notes on the Equations

These equations are used to compute the time evolution of 1D (vertical, subsurface) profiles for (1) soil moisture, θ, (2) pressure head, ψ, (3) hydraulic conductivity, K and (4) vertical flow rate, v. TopoFlow solves these equations separately to get time-evolving profiles for every grid cell in a DEM. The result is a 3D grid for each of these four variables that spans the unsaturated zone. The third equation above just defines a variable that is used in the 4th and 5th equations, so the coupled set constitutes 4 equations to be solved for 4 unknowns. These equations can be combined into one nonlinear, parabolic, second-order PDE (partial differential equation) known as the one-dimensional Richards' equation.

The infiltration rate is simply the vertical flow rate at the ground surface, denoted by v0.

Soil moisture is simply another term for the water content in the case where the porous medium is a soil.

Subscripts in the first two equations indicate partial derivatives with respect to the vertical coordinate, z. Note that z is the vertical distance below the ground surface, in meters.

The last three equations (which reduce to 2 equations) are empirical and are called soil characteristic relations. Changing the parameters in these equations allows different soil types to be modeled. Typical values have been tabulated for many standard soil types, and the dialog updates with these values when you make a selection from the droplist labeled Closest standard soil type .

The soil characteristic relations used by TopoFlow (and shown here) are those of the transitional Brooks-Corey (T-BC) method. This method combines the advantages of the well-known Brooks-Corey and van Genuchten methods, as explained by Smith (2002). Note that the T-BC equations have ψ = -ψA and K = Ks when the soil becomes saturated, that is, when θ = θs.

Soils in nature do not have water contents lower than that corresponding to hygroscopic water (Dingman, 2002, p. 236-238). At this extreme dryness water is absorbed directly from the air. The corresponding pressure head (tension head) is -31,000 cm. TopoFlow uses the T-BC method described in the previous paragraph to compute the minimum attainable water content for a given soil type and parameter set from this minimum attainable pressure head. Any water content below this value is considered unavailable. Because of this, TopoFlow does not allow the initial water content of a given soil type to be set to less than this value. Note that the permanent wilting point is associated with a pressure head of -15,000 cm. Water content associated with values less than this is considered to be unavailable to plants. The field capacity is the soil water content associated with a pressure head of about -340 cm and is an estimate of the water content that can be held against the force of gravity. TopoFlow uses this value as the default setting for the initial water content of any given soil layer. The default setting for residual water content is computed as shown in the fourth equation for a reference pressure head of -1,000,000 cm. See Dingman (2002, p. 237).

The unsaturated zone extends from the soil surface all the way down to the water table, where θ = θs. The water table is tracked by TopoFlow as a moving boundary, and may actually rise up to the land surface in some places. When this occurs, a subsurface seepage term is included in the water balance equations. Typically, ψA is set to zero so that the pressure head is zero at the water table, but it can also be used to model hysteresis effects, as explained by Smith (2002).


References

Dingman, S.L. (2002) Physical Hydrology, 2nd ed., Prentice-Hall, Upper Saddle River, New Jersey, 646 pp.

Smith, R.E. (2002) Infiltration Theory for Hydrologic Applications, Water Resources Monograph 15, AGU, 212 pp.