Note: This is documentation for version 5.0 of Source. For a different version of Source, select the relevant space by using the Spaces menu in the toolbar above
Groundwater Numerical Model for 1-Dimensional Flow (GN1D) - SRG
Contents
Background
The Groundwater Numerical model for 1-Dimensional flow (GN1D) simulates the stream – aquifer flux and groundwater heads for a saturated connection between a river and connected unconfined aquifer. The model has been developed to meet three primary needs:
- Improved calibration of rainfall-runoff models via better representation of baseflow
- Representation of floodplain interactions
- Evaluation of conjunctive water use
It operates at the scale of a river reach and represents the following types of processes:
- Stream – aquifer flux
- Groundwater levels
- Pumping and injection
- Diffuse infiltration
- Ponded infiltration
- Evapotranspiration
The GN1D model is not intended to capture the details of regional groundwater flow or small-scale variations in groundwater levels across a catchment. These types of applications require the use of a more complex 3D groundwater flow model, such as MODFLOW (Harbaugh, 2005).
Scale
The model is applied at the link scale. The model assumes a daily time-step.
Principal Developer
eWater.
Scientific Provenance
The GN1D model contains two primary components:
- A finite difference scheme for groundwater flow modelling.
- The Mein-Larson equations (Mein and Larson, 1973) for modelling infiltration.
Both of these methods are well established. A project demonstrating their implementation in the Source GN1D model is underway.
Version
Source version 3.8.12 onward.
Availability/conditions
The GN1D model is automatically included in Source Public Beta. It is not available for Source Public.
Dependencies
The GN1D model is created on a river link and it requires that the link has a rating curve to translate river flow to stage. A rating curve can be created from the GN1D model cross section, discussed below.
Theory
Introduction
The model conceptualises a river reach as a straight line and assumes that the subsurface soil properties are homogeneous across the model domain. The geometry of the riverbed and the surrounding land surface are defined by a cross-section that lies perpendicular to the reach.
The regions either side of the river reach are divided into a number of subsections that run parallel to the reach (see Figure 1). Within each subsection, the subsurface is conceptualised as an unsaturated zone overlying a saturated zone:
- The unsaturated zone simulates the impacts of infiltration and evapotranspiration on soil moisture.
- The saturated zone simulates groundwater levels in an unconfined aquifer, including the impacts of stream – aquifer fluxes and other losses and gains.
The geometrical conceptualisation and the unsaturated and saturated zone models are described in more detail below.
Figure 1. Conceptual illustration of the GN1D model structure, showing one bank of the river reach only.
Geometrical Conceptualisation
River Link Topology: Cross-sections
The geometry of the riverbed and the surrounding region is defined by specifying a cross-section that lies perpendicular to the river link (see Figure 2). This cross-section represents the average link topology, and it is applied uniformly along the link.
It is important to note that the cross-section defines the coordinate system used by other components of the groundwater model, such as the locations of pumps and observations bores.
At the same time as defining the cross section, the user must also define:
- The position of the left hand and right hand riverbanks. The groundwater model uses the positions of the left and right banks to determine when overbank flow events occur.
- Manning's N and the link slope. Manning's N and the slope parameters can be used, if desired, to generate a streamflow rating curve from the cross section. Alternatively, the user can calculate and input their own rating curve.
Figure 2. Cross-section conceptualisation.
Groundwater Model Geometry
The groundwater model divides the cross-section into a number of user-defined subsections (illustrated in Figure 3). The subsections run parallel to the river link. They determine the spatial resolution of the groundwater model and form the basis for the 1-D finite difference scheme (saturated zone model).
- The subsections start either side of the riverbanks, as the groundwater model does not extend underneath the river.
- The left and right extents of the subsections are implemented as no flow boundaries.
- The depth of the subsections (z-coordinate) is defined by the user and is implemented as a no-flow boundary.
- The initial groundwater table height is constant across the subsections and is set by the user.
Figure 3. Example illustrating the representation of a cross-section within the groundwater model using 4 subsections per bank (grey shading).
Subsection height calculation
The height of the ground surface within each subsection is calculated by taking the average of the cross-section elevations within the subsection. The horizontal extent of the groundwater model domain may extend beyond the end-points of the cross-section. If so, the gradient of the last two points in the cross-section is linearly extrapolated to enable the average ground height to be calculated (as illustrated in Figure 3).
Unsaturated Zone
Introduction
Each subsection in the groundwater model is divided into an unsaturated zone and a saturated zone, where the unsaturated zone is the region between the ground surface and the groundwater table. The unsaturated zone models the following processes:
- Diffuse and ponded infiltration
- Evapotranspiration
Infiltration
The infiltration module estimates infiltration and runoff for two types of infiltrative processes:
- Diffuse infiltration, which refers to cases where the depth of ponded water at the surface is negligible, e.g. rainfall.
- Ponded infiltration, e.g. surface flooding.
The infiltration module is based on the Mein and Larson (1973) formulation of the Green-Ampt equations. It has four parameters that must be entered by the user. These are:
- Kvsat – vertical saturated hydraulic conductivity
- – the initial soil moisture content (volumetric fraction)
- – the saturated soil moisture content (volumetric fraction)
- Save – average capillary suction at the wetting front
The parameters of the Mein-Larson equations can be related to measurable soil characteristics and the equations are simple and easy to solve. Two constraints of the approach are that:
- The Mein-Larson equations do not describe the redistribution of soil moisture at the cessation of rainfall/ponding.
- As with most infiltration models, it is difficult to account for spatial and temporal parameter variability and vegetative characteristics.
Diffuse Infiltration
Diffuse infiltration takes as input a time series of rainfall (or other diffuse source) intensities. The Mein and Larson equations model diffuse infiltration in two stages as illustrated in Figure 4 (Mein and Larson, 1973). The first stage predicts the volume of water that will infiltrate before the soil surface becomes saturated. From this point, surface runoff begins and the infiltration capacity is predicted by the Green-Ampt equation.
Figure 4. Infiltration rate as a function of time for steady rainfall.
Additional assumptions:
- The depth of ponded water at the surface is negligible
- The rainfall is of constant intensity until run-off begins.
- The soil is a homogeneous medium.
Ponded Infiltration
Ponded infiltration is modelled similarly to diffuse infiltration (Elhanafy et al., 2007), except that:
- The soil surface is assumed to be instantly saturated.
- The infiltration rate is driven by the additional force of the ponding depth at the ground surface.
The method assumes that the ponding depth does not change during a time step, which is reasonable if the infiltrated volume is a small in proportion to the ponded volume. As with diffuse infiltration, the soil is assumed a homogeneous medium.
Evapotranspiration
The evapotranspiration (ET) process model calculates the amount of water that is lost from the unsaturated and saturated zones due to evapotranspiration (the combination of plant transpiration and bare soil evaporation). The ET model takes as input a time series of potential evapotranspiration (PET). PET is first subtracted from the volume of water stored in the unsaturated zone above an extinction depth. The extinction depth, which is defined by the user, is the distance below the ground surface at which ET falls to zero. If some evaporative demand remains and the groundwater table is above the extinction depth, the residual PET is then applied to the saturated zone. Groundwater evapotranspiration is assumed to decrease with increasing depth to the water table and is calculated as:
Soil Moisture Re-distribution
After applying infiltration and evapotranspiration, the soil moisture content of the unsaturated volume is calculated. If the soil moisture is greater than field capacity, then the volume of soil water that is above field capacity is allowed to drain to the saturated zone at a rate limited by the vertical saturated hydraulic conductivity.
Saturated Zone
The saturated zone represents the groundwater aquifer. The 1-dimensional equation for horizontal groundwater flow is given by:
(1) |
where:
- S is the specific storage of a confined aquifer (or specific yield in an unconfined aquifer),
- T Transmissivity
- F is the next flux from vertical losses and gains
- h is the groundwater head
- t is time
- x is the horizontal distance
Equation 1 is derived from Darcy's law using conservation of mass and assuming that transmissivity is homogeneous (constant) and isotropic (equal in all directions).
Figure 5. Illustration of the finite difference scheme.
The GN1D model solves equation (1) for groundwater head using an explicit finite difference scheme (illustrated in Figure 5):
(2) |
Where:
- is the groundwater head in subsection n at time i
- Ksat saturated hydraulic conductivity
- is the aquifer thickness
- is the specific storage of a confined aquifer (or specific yield in an unconfined aquifer)
- Net vertical losses and gains to subsection n at time (i+1)
- is the time step
- is the distance
The initial groundwater head
is set by the user and is uniform for all subdivisions, i.e. for where is the number of divisions.Equation (2) is valid for lateral flow in a confined aquifer and can also be applied to an unconfined aquifer providing that the spatial variation in head is small compared to the aquifer thickness. In the context of estimating stream – aquifer interactions, this requirement is generally met if the stream is not deeply incised compared to the unconfined aquifer thickness. The phreatic surface (i.e. the groundwater table) can be highly variable close to things such as extraction bores. However, at the spatial scale at which the model is designed to operate, the localised impacts of groundwater extractions are likely to be less significant.
Equation (2) becomes unstable if the following criteria is not met:
(3) |
Accounting for Soil Moisture
Equation (2) describes the head fluctuations in the saturated zone assuming that the aquifer material is at field capacity. If the unsaturated zone is not at field capacity, the expected groundwater head rise will differ. To ensure that mass balance is maintained, the groundwater model employs a post-processing step to adjust the head fluctuations calculated by the finite difference scheme to account for the soil moisture content of the unsaturated zone.
Parameters
The GN1D model settings and parameters are described in Table 1. Typical values of the parameters for the Mein-Larson infiltration equations are given in Table 2.
Table 1. GN1D model settings and parameters.
Component | Variable | Symbol | Typical Units | Description |
---|---|---|---|---|
Geometry | River to no-flow boundary | None | m | The distance from the left and right hand river banks to the groundwater no-flow boundary |
Geometry | Divisions per bank | None | The number of subdivisions per bank | |
Geometry | Aquifer floor | None | m | The Y-coordinate of the aquifer basement (no-flow boundary). Uses the coordinate system defined by the cross section. |
Saturated Zone Model | Transmissivity | m2/day | The rate at which groundwater flows horizontally through an aquifer. The transmissivity for horizontal flow | |
Saturated Zone Model | Specific yield | fraction (0…1) | The volumetric fraction that an aquifer will yield when all water is allowed to drain out of it under the force of gravity: | |
Saturated Zone Model | Initial groundwater table height | m | The initial groundwater table height at time step 0. | |
Unsaturated Zone Model | Saturated soil moisture content | fraction (0…1) | The volumetric fraction of water that an aquifer will hold when it is saturated, equivalent to porosity. | |
Unsaturated Zone Model | Initial soil moisture content | fraction (0…1) | The initial (time step 0) soil moisture content of the unsaturated zone as volumetric fraction of the total volume. | |
Unsaturated Zone Model | Saturated vertical hydraulic conductivity | m/day | The rate at which water moves vertically through saturated media. | |
Unsaturated Zone Model | Average capillary suction at the wetting front | m | The average capillary suction (head) at the wetting front as defined by Mein and Larson (1973). | |
Unsaturated Zone Model | Extinction depth | None | m | The maximum depth below the ground surface at which evapotranspiration can occur. |
Table 2. Average values of Green-Ampt parameters by USDA soil texture class (adapted from Rawls et al., 1983).
USDA Soil-Texture Class | Hydraulic Conductivity (cm/day) | Wetting-Front Suction Head (cm) | Porosity (Volumetric fraction) | Water Retained at Field Capacity (Volumetric fraction) | Water Retained at Wilting Point (Volumetric fraction) |
---|---|---|---|---|---|
Sand | 144.48 | 4.9 | 0.437 | 0.062 | 0.024 |
Loamy Sand | 35.97 | 6.1 | 0.437 | 0.105 | 0.047 |
Sandy Loam | 13.11 | 11 | 0.453 | 0.19 | 0.085 |
Loam | 3.96 | 8.89 | 0.463 | 0.232 | 0.116 |
Silt Loam | 7.92 | 16.99 | 0.501 | 0.284 | 0.135 |
Sandy Clay Loam | 1.83 | 22 | 0.398 | 0.244 | 0.136 |
Clay Loam | 1.22 | 21.01 | 0.464 | 0.31 | 0.187 |
Silty Clay Loam | 1.22 | 27 | 0.471 | 0.342 | 0.21 |
Sandy Clay | 0.61 | 24 | 0.43 | 0.321 | 0.221 |
Silty Clay | 0.61 | 29.01 | 0.479 | 0.371 | 0.251 |
Clay | 0.3 | 32 | 0.475 | 0.378 | 0.265 |
Input Data
Each time step, the GN1D model requires the river stage, which is calculated by the link using a rating curve. Other inputs are configured by the user and are summarised in Table 3.
Table 3. GN1D model data requirements, including typical units.
Component | Data Requirements | Comments |
---|---|---|
Cross section | Cross section X-Y coordinates (m) | Required. |
Potential Evapotranspiration (PET) | Time series of values (mm/day) | Optional |
Extraction Bore | Time series of values (ML/day) | Optional |
Injection Bore | Time series of values (ML/day) | Optional |
Diffuse Recharge | Time series of values (mm/day) | Optional |
Direct Diffuse Recharge | Time series of values (mm/day) | Optional |
Observation Bore | X-coordinat |