Subsurface Flow → Darcy, Surface-parallel Method

The input variables used for modeling horizontal subsurface flow in the saturated zone via Darcy's Law are defined as follows:

Ks = saturated hydraulic conductivity [m / s]
Sy = specific yield or drainable porosity [unitless]
   (less than or equal to the porosity, φ, see Notes)
thickness = soil layer thickness [meters]

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 have 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.

It is possible to have six different soil layers, each with different parameter values. The soil layer number is indicated at the top of the dialog and the Next and Back buttons in the top section can be used to select a layer. In order to use fewer than six layers, a layer can be made impermeable by either setting Ks = 0 or θs = 0 for that layer. All layers below an impermeable layer are then ignored.

Note that the thickness of a layer can vary spatially; this requires choosing "Grid" from the droplist that specifies the thickness "Type".

Important Note: The Depth to bedrock settings are not yet used in this beta version.


Equations Used by the Darcy, Surface-parallel Method

Q(k) = vhor(k) * y(k) * b   = horizontal discharge out of grid cell in kth layer [m3 / s]
vhor(k)   = Ks(k) * SH = horizontal flow rate [m / s]   (from Darcy's Law)
SH = (Hx2 + Hy2)1/2 = local slope of the water table [unitless] (see Notes below)
ΔH / Δt = [ R - Σ Q(k) / (Δx Δy) ] / Sy = integral form of mass balance [m / s] (see Notes below)

In the first equation, y(k), is the saturated thickness of the kth layer and b is the flow width perpedicular to the direction of flow. Note that b is the same for all layers, by assumption, and depends on the D8 flow direction. Both y(k) and b have units of meters.

In the last equation, R is the recharge rate, or the vertical infiltration rate just above the water table, and Δx and Δy are th DEM grid cell dimensions. The notation Σ Q(k) denotes the sum of all flows that are either into a grid cell from neighboring cells or out of the cell. This corresponds to a boundary integral of the flow field around the grid cell boundary.


Notes on the Equations

This method makes a number of simplifying assumptions that are important to understand. It is intended for modeling unconfined aquifers , which have a water table, H(x,y), as their upper boundary. The Dupuit assumptions are often used for unconfined aquifers and are used here, namely:

  1. flow lines are horizontal and equipotential lines are vertical and
  2. the horizontal hydraulic gradient is equal to the slope of the free surface (the water table) and is invariant with depth.
In addition, the saturated hydraulic conductivity is assumed to be isotropic throughout any given layer, so that a single value, Ks, can be used in both the x and y coordinate directions. Note that this does not prevent Ks from varying spatially as a function of x and y.

Another key assumption of this method is that the subsurface flow direction is always the same as the surface flow direction to within an angle of 45 degrees; that is, the D8 flow direction found from Z(x,y) (provided as a DEM) is always used for the subsurface flow. Put another way, the gradient vectors for H(x,y) and Z(x,y) are assumed to differ by less than 45 degrees. (Notice, however, that SH, which is the magnitude of the gradient vector for H(x,y), is used in Darcy's Law.) This assumption allows an integral vs. differential form of mass balance to be applied to the soil column that lies below any given DEM grid cell. That is, instead of numerically solving a second-order PDE known as the Boussinesq equation by the finite-difference method, the amounts flowing across the soil column's six faces and the amount stored within the soil column and the requirement of mass balance are used to compute changes to the height of the water table. Since the subsurface flow direction is always taken to be the same as the surface flow direction, this method does not allow subsurface flow to cross beneath a drainage divide.

The water table, H(x,y), is a dynamic surface (a moving boundary), and may actually rise up to the land surface in some places. Both the land surface, Z(x,y), (given by a DEM) and the water table are represented as grids (or 2D arrays). Wherever H(x,y) is greater than Z(x,y), mass balance is used to compute a seepage rate from the subsurface that contributes to the surface runoff.

Control volumes just above the water table have horizontal flow into or out of their 4 vertical faces as well as vertical flow into (or out of) their top face. The top face has the dimensions of a DEM grid cell. Horizontal flow rates (in the saturated layers below the water table) are determined by Darcy's Law, while the vertical flow rate (in the unsaturated zone just above the water table) is determined by the method that the user has selected for modeling the infiltration process. If the net flow into the control volume is positive, then the water content increases (up to a maximum of θs), the saturated thickness increases (up to a maximum of the layer thickness) and the water table rises locally as required to conserve mass. Similarly, if it is negative, the water table falls.

The specific yield or drainable porosity refers to the water content, as a volumetric fraction, that can be made to flow out of a control volume in response to a water table gradient. Note that the porosity, φ, is a volumetric fraction that defines the amount of void space in a given porous medium that water could potentially occupy. However, there are some air spaces that remain unfilled even when a porous medium is saturated, so the saturated water content , θs, is always somewhat smaller than the porosity. Similary, it is generally not possible to remove all of the water in a given porous medium through the application of a head gradient; there is generally some residual amount that remains within the pores, unless the soil is desicated by some other means. This amount is known as the residual water content, θr. It follows that Sy = (θs - θr).

Note: In the current implementation of this method, discharges from all layers are summed in order to determine the net horizontal flow into soil columns (cylindrical control volumes) that extend up to the land surface. The water table height is then adjusted accordingly to conserve mass. If the net flow into a deep, saturated layer is positive (perhaps due to convergent horizontal flow lines) then this approach causes the excess water to be transferred up to the water table in a single time step. Depending on how Ks varies with distance below the ground surface, z, and how thick the layers are, this may not be physically realistic.

The values of Sy, θs and θr are all dimensionless and must be between 0 and 1. Freeze & Cherry (1979, Table 2.4) give the following table with typical values of θs for various soil types.

Gravel 0.25 - 0.40
Sand 0.25 - 0.50
Silt 0.35 - 0.50
Clay 0.40 - 0.70
Fractured basalt 0.05 - 0.50
Karst limestone 0.05 - 0.50
Sandstone 0.05 - 0.30
Limestone, dolomite 0.00 - 0.20
Shale 0.00 - 0.10
Fractured cryst. rock 0.00 - 0.10
Dense crystalline rock 0.00 - 0.05


References

Anderson, M.P. and W.W. Woessner (2002) Applied Groundwater Modeling: Simulation of Flow and Advective Transport, Academic Press, An Imprint of Elsevier, New York, 381 pp. (see chapter 2 and pages 15 and 42)

Freeze, R.A. and J.A. Cherry (1979) Groundwater, Prentice-Hall, Englewood Cliffs, N.J., 604 pp.