Model help:TopoFlow-Saturated Zone-Darcy Layers
TopoFlow-Saturated Zone-Darcy Layers
This module is the saturated Zone process component (Darcy's law, multiple soil layers) for a D8-based, spatial hydrologic model
Model introduction
This process component is part of a spatially-distributed hydrologic model called TopoFlow, but it can now be used as a stand-alone model.
Model parameters
Uses ports
• Infil (Infiltration)
• Channels (surface water flow in a network of channels)
Provides ports
• Satzone (Subsurface flow in saturated zone)
• Configure (tabbed dialog GUI to change settings)
• Run (only if used as the Driver)
Main equations
- Horizontal discharge out of grid cell in kth layer
[math]\displaystyle{ Q \left ( k \right)= v_{hor} \left ( k \right) \ast y\left ( k \right) \ast b }[/math] (1)
- Horizontal flow rate (From Darcy's Law)
[math]\displaystyle{ v_{hor}\left ( k \right)= K_{s}\left ( k \right) \ast S_{H} }[/math] (2)
- Local slope of the water table
[math]\displaystyle{ S_{H}= \left ( H_{x}^2 + H_{y}^2\right)^{\frac{1}{2}} }[/math] (3)
- Integral form of mass balance
[math]\displaystyle{ \Delta H / \Delta t = [ R - \Sigma Q \left ( k \right) / \left ( \Delta x \Delta y \right)] / S_{y} }[/math] (4)
| Symbol | Description | Unit |
|---|---|---|
| Ks | saturated hydraulic conductivity | [m / s] |
| Sy | specific yield or drainable porosity | [-] |
| Q | horizontal discharge | m^3 / s |
| vhor | horizontal flow rate | m / s |
| y | saturated thickness | m |
| b | flow width | m |
| SH | local slope of the water table | - |
| Hx, Hy | water table at the upper boundary | m |
| R | recharge rate | |
| ΔH | change in water surface | m |
| Δt | change in time | sec |
| Sy | drainable porosity, equals to θs - θr |
Notes
Notes on Input Parameters
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. Each layer has its own tab at the top of the dialog. 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".
The Depth to bedrock settings are not yet used in this beta version.
Notes on the Equations
All variables and their units can be seen by expanding the Nomenclature section above.
In the first equation, y(k), is the saturated thickness of the kth layer and b is the flow width perpendicular 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 the 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.
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
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. Similarly, 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.
| Soil type | θs |
|---|---|
| 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 |
Examples
An example run with input parameters, BLD files, as well as a figure / movie of the output
Follow the next steps to include images / movies of simulations:
- Upload file: https://csdms.colorado.edu/wiki/Special:Upload
- Create link to the file on your page: [[Image:<file name>]].
See also: Help:Images or Help:Movies
Developer(s)
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.
Links
Related Help Pages
- Model_help:TopoFlow-Infiltration-Green-Ampt
- Model_help:TopoFlow-Infiltration-Smith-Parlange
- Model_help:TopoFlow-Infiltration-Richards 1D
- Model_help:TopoFlow-Soil Properties Page
Model Metadata
