Template:GUM8

 

8 - Lateral Groundwater Flow Modeling in the Saturated Zone

GSSHA User's Manual

8.1 General

Simulation of saturated 2-D groundwater flow is selected by placing the GW_SIMULATION card in the project file. When groundwater simulations are performed, the starting water surface elevations (m) must be provided in a file specified with the WATER_TABLE card; the file containing bedrock elevations (m) is specified with the AQUIFER_BOTTOM card. The groundwater time step must also be specified (s) using the GW_TIMESTEP card, typically 1 hour, 3600 s. A variety of boundary conditions can be specified with the GW_BOUNDFILE card. Saturated groundwater may be simulated using either the GAR or RE methods to calculate infiltration. In GSSHA v5 and above, multi-layer GA can also be used to simulate infiltration for groundwater models. The saturated groundwater solution was developed in conjunction with the RE unsaturated groundwater solution, and is intended to be used with RE. The GAR method was linked to the saturated groundwater model to provide a simple, quick, and approximate groundwater recharge estimate. This method is a useful screening tool, and can be used to help define parameter values before applying the more rigorous RE method. The multi-layer GA model was extended for use with continuous simulations and as a result, was also extended for use with the saturated groundwater model. Whether GAR or RE or multi-layer GA is used affects how parameter values are assigned, as described below.

8.2 Formulation

Trescott and Larson (1977) described the solution to the two-dimensional free surface groundwater problem, and the efficiency of various solvers. Their methods were largely followed in the development of this portion of the code; an exhaustive coverage need not be presented here. The overall approach, differences in approach, and integration into the GSSHA model are presented.

The controlling equation, as developed by Pinder and Bredehoeft (1968), is:

(40)

where:

T is the transmissivity (m²/s),

h is the hydraulic head (m),

S is the storage term (dimensionless),

and W is the flux term for sources and sinks (m/s).

It is assumed that off diagonal terms are not important and that transmissivity can be expressed as the product of the saturated hydraulic conductivity of the media (K) and the depth of the saturated media (b). For the free surface problem the head is the surface water elevation (E_ws).

(41)

This equation can be represented using a block-centered finite difference five-point implicit scheme as:

(42)

This representation varies from the original representation of Trescott and Larson (1977) in that the transmissivities and the storage terms are both time dependent and calculated implicitly using Picard iteration.

8.3 Solution

The equation is solved by successive overrelaxation by lines (LSOR) (for example Tannehill et al., 1997). LSOR was shown by Trescott and Larson (1977) to be capable of solving a variety of difficult groundwater problems, though not necessarily being the most efficient method. With LSOR, the two-dimensional problem is linearized by solving by rows or by columns. The user specifies solution by rows or by columns with the GW_LSOR_DIR project card, and the choice is made to align the direction of solution with the principal direction of flow, x (argument - 1) or y (argument - 2) (Trescott and Larsen, 1977). LSOR is an iterative method; the user selects the groundwater convergence criteria (m) with the GW_LSOR_CON project card. A typical value, the default, is 10^-5 m. The user also selects an over relaxation coefficient, ω, with the GW_RELAX_COEF card. During each iteration the head in each cell is adjusted by the over relaxation coefficient, ω, such that:

(43)

where k denotes the iteration number. For ω greater than 1.0, the next head is projected out on a line determined from previous iterations. This speeds the convergence process but may also reduce stability. Typically, a value of ω of about 1.2 results the fastest solution. For very difficult problems ω may need to smaller than 1.0, and the solution is said to be underrelaxed.

For each iteration the transmissivities in both the x and y directions are calculated based on the updated saturated depth, b, determined from the heads, so that:

(44)

where K_gw is the lateral hydrualic conductivity of the porous media. This is a variation from Trescott and Larson (1977) who calculate the transmissivity based on the value of b from the last, n^th, time step such that Tⁿ⁺¹ = K_gwbⁿ. The storage term, S, in the equation is used to represent the change in volume with respect to the change in head:

(45)

where V is the volume. In two-dimensional applications it is common practice to represent the storage as the porosity of the saturated groundwater media. When RE is used to define the unsaturated zone the storage term, S, is not a constant but also depends on the head. The storage term is updated each iteration. The storage term is set to the porosity of the unsaturated cell above or below the groundwater surface elevation, depending on whether the groundwater is rising or falling. For time steps when the groundwater moves more than one unsaturated zone cell, the storage term is calculated as a bulk storage term:

(46)

where ΣΔz is the distance that the groundwater surface is anticipated to move during the n+1 time step or k+1 iteration. The anticipated change in the groundwater surface elevation is determined from the source term, W, containing all fluxes in the cell and the storage capacity, S, of surrounding cells. An internal time step limitation attempts to limit the movement of the groundwater surface to a single unsaturated zone cell. This added step helps maintains overall mass conservation. The actual storage available is Θ_s--Θ, but because water content is a discrete value in each cell, calculating storage as Θ_s-Θ over ΣΔz can result in a lack of convergence of the groundwater solution, which is implicitly calculating h, S, b, and T simultaneously. As discussed below in the section on coupling of processes, the difference in water volume between Θ_s and Θ_s-Θ is added to the source term during the next groundwater update, such that mass is conserved. This may result in a small time lag in the groundwater response.

When GAR is used to provide recharge estimates to the groundwater model, the storage term is constant, equal to the value moisture deficit of the soil, effective porosity minus soil moisture. The moisture deficit is updated every rainfall event.

8.4 Assignment of Parameter Values

The value of K_gw (cm/hr) may be specified as a constant value with the GW_UNIF_HYCOND or a map of spatially distributed values specified with the GW_HYCOND_MAP card. The porosity of the groundwater media below the specified unsaturated zone can be a uniform value, specified with the GW_UNIF_POROSITY card, or can be distributed by specifying a map of spatially distributed values using the GW_POROSITY_MAP. If neither uniform or distributed values of these parameters are specified with the above cards, the values will be set to the those of the last soil layer of the soils in the unsaturated zone specified in the SOIL_TABLE_INPUT_FILE or Mapping Table files. When GAR is used to provide estimates of recharge, only hydrualic conductivity need be specified, as the storage is computed using the GAR infiltration parameters.

8.5 Boundary Conditions

Each cell in the active domain is assigned a boundary condition with a map of integer values specified with the GW_BOUNDFILE project card. Boundary conditions around the border of the watershed may be constant head (type 2), constant flux (type 0) or a combination of both. Flux boundaries are represented by setting the transmissivity of flow boundary cells to zero and entering the flux in the source/sink term, W. For all non-head boundary cells (type 1), the unsaturated zone provides a flux to the source/sink term, W. The lower boundary is a flux term (cm/hr) that represents a leaky aquifer, specified with the GW_LEAKAGE_RATE project card. The default is zero. This single user-specified value is used in all cells.

Possible internal boundary conditions include constant head (type 2), constant flux (type 6), variable flux (type 3), river head (type 5) and river flux (type 4), and reservoir cell (type 7). The constant flux and variable flux boundaries are used to represent pumping wells of constant and variable rate, respectively. Water added/subtracted by pumping is added to the source/sink term, W. The possible boundary types are:

0 – no flow

1 – regular infiltration cell (no special boundary condition)

2 – specified head, taken from WATER_TABLE file

3 – dynamic flux (well)

4 – stream cell with calculated flux between stream node and groundwater cell

5 – stream cell with a specified head, value from stream routing solution

6 – static flux (well) value taken from the GW_FLUX_BOUNDTABLE

7 - reservoir cell

8.5.1 Stream Boundary Cells – Type 4 & 5

When 1-D channel flow is simulated, grid cells containing stream channels may be represented as either head (type 5) or flux boundaries (type 4). For a head boundary, the elevation of the water surface in the stream node is used as a specified head in the groundwater solution. For a river flux boundary condition the flux between the groundwater cell and stream node is calculated during the channel routing update. Fluxes accumulated over multiple channel routing updates are accumulated and then added to the groundwater cell in the source/sink term, W. Additional inputs are required to simulate groundwater/stream interactions:

• thickness of the riverbed material, M_rb, (cm), and
• vertical hydraulic conductivity of riverbed material, K_rb (cm/hr).

These may be specified as uniform values using the M_RIVER and K_RIVER project cards, respectively. These may also be distributed along the stream network in the channel input file (cip file), specified with the CHAN_INPUT card, as described in Section 5.1 Channel Routing.

The flux throught the river bed is calculated based on the difference in elevation between the groundwater surface and the water surface elevation in the stream node. McDonald and Harbaugh (1988) developed a simple method to compute flows between the channel network and saturated groundwater based on the Darcy equation. If the groundwater surface elevation is above the riverbed elevation but below the river water surface elevation, the river discharges to the groundwater:

(47)

where:

f = per unit area discharge (m/s),

K_rb = hydraulic conductivity of the river bed material (cm/hr),

M_rb = depth of the river bed material (cm),

E_r = elevation of the river water surface (m), and

E_ws = elevation of the groundwater water surface (m).

The negative flux indicates that the discharge is into the groundwater. If the groundwater surface elevation is below the riverbed elevation the river leaks to the groundwater at the rate

ƒ = -Krb

(48)

When the STREAM_LOSS card is used in the absence of a known or calculated water table location, Equation 48 is used to compute the stream loss.

If the groundwater surface elevation is above the river water surface elevation, the groundwater discharges to the river according to Equation 46. In this case, the flux will be positive, indicating flow is from the groundwater to the river. The unit flux, ƒ, is multiplied by the top width and length of the channel segment. In calculating the top width an effective depth is defined. If the water surface elevation in the channel is higher than the water surface elevation of the groundwater, the channel depth is used in the calculation. When the groundwater level is higher than the water level in the stream, the effective depth used is the groundwater elevation minus the bed elevation.

The channel side slope properties are assumed to be the same as those of the channel bed. If properties vary significantly between the bed and the side slope an average or weighted average values should be input. This approach is highly empirical and the riverbed property values will require calibration to the system being modeled.

8.5.2 Static Pumping Wells – Type 6

Static pumping rate wells may be placed in any cell in the watershed domain. Pumping wells are located by placing a value of 6 in the GW_FLUX_BOUNDTABLE at each desired well location. The GW_FLUX_BOUNDTABLE specifies the name of a file with the pumping rates. The GW_FLUX_BOUNDTABLE has the following format.

I location well 1          J location well 1          Pumping rate (m³/d) well 1
I location well 2          J location well 2          Pumping rate (m³/d) well 2
I location well 3          J location well 3          Pumping rate (m³/d) well 3
…
…
…
I location well N-1    J location well N-1     Pumping rate (m³/d) well N-1 I location well N        J location well N         Pumping rate (m³/d) well N

Values are separated by spaces or tabs, not commas. Groundwater extractions have a positive pumping rate; injections have a negative pumping rate. The number and I, J locations of the wells in table must match the number and location of wells in the GW_BOUNFILE file.

8.5.3 Reservoir Cells - Type 7

Any cell that is within an active reservoir is assigned a boundary condition of 7. As the reservoirs change size within the GSSHA model, type 7 boundary conditions are assigned internally to GSSHA. If a cell has a type 7 boundary a flux between the reservoir and the groundwater will be calculated as the Darcy flux based on the head between the reservoir and the groundwater, the depth of sediments in the bottom of the reservoir, and the hydraulic conductivity of the sediments. The equations are the same as those described for groundwater channel interaction. Values of hydraulic conductivity for each cell are taken as the hydraulic conductivities prescribed for the infiltration model. The thickness of the bed sediments is specified with the MLAKE card in the CHAN_INPUT file. A uniform value, for all lakes, can be specified in the project file with the M_LAKE card. The value is specified in cm. There is no default value. A different value can be entered for any or all of the lakes. Groundwater interaction will not be calculated for reservoir links without the MLAKE card in the CHAN_INPUT file or the M_LAKE card in the project file.

8.5.4 Static and Dynamic Wells using the Mapping Table

In addition to specifying static wells above, both static and dynamic wells can be created by using the WELL_TABLE table of the mapping table file. An index map of wells with unique pumping characteristics is created (0's elsewhere.) This index map and the associated table is then used to modify the groundwater boundary map. This table has two parameters. The first is a flag (1 for yes and 0 for no) to indicate whether the well is dynamic (1) or static (0). If the well is static, the second parameter is the pumping rate (positive values) or injection rate (negative values). If the well is dynamic then the second value refers to a time series index; you will need to also set up the time series index table.

To create a dynamic well, there are four steps that need to be followed.
1) set up a time series of the pumping rates (positive values) or injection rates (negative values) (see Time series file format)
2) reference the time series in the project file (TIME_SERIES_FILE "filename.ts", as indicated in the above link to time series) preferably before the mapping table.
3) set up the time series index mapping table
4) set up the WELL_TABLE table, with the is_dynamic value set to 1 and the value set to the time series index for that particular time series.

8.7 Coupling of the Saturated Zone Model to the GAR Infiltration Model

The GAR routine described in Section 7.4 calculates infiltration and soil moisture based on the movement of multiple sharp wetting fronts that occur with rainfall burst followed by rainfall hiatus. The multi-layer GA model described in Section 7.5 calculates infiltration and soil moisture in a user specified distributed three layer soil model. The infiltration calculated by GAR or multi-layer GA can be used to provide a rough estimate of groundwater recharge to the 2-D saturated groundwater flow simulation. The groundwater recharge for each time step is equal to the infiltration computed by GAR or multi-layer GA. In the model the groundwater recharge is computed each time step as

(52)

where: R is equal to the recharge m³, F is the total infiltration (m), and A is the area of the overland flow cell (m²). The simple two layer soil moisture accounting routine is used to calculate soil moisture. When the groundwater elevations exceed the ground surface elevation, infiltration calculations for the cell cease, and the groundwater surface exchange is calculated as described above. Anytime exfiltration occurs, the infiltration and overland flow processes are initiated if they are not already active. These processes remain active as long as exfiltration occurs and until all water on the overland flow plane stops flowing and infiltrating. Infiltration is not calculated for cells in which exfiltration is occurring.

When using the GA models to couple the surface soils to the groundwater table the properties of the unsaturated groundwater media existing between the specified soil column, SOIL_DEPTH for GAR, and water table are assigned the groundwater porosity. The soil saturation (fraction) for this unsaturated soil may be specified by the user with the SINGLE_UNSAT_SAT card. The default value is 0.75, meaning the soil is 75% saturated. This corresponds, roughly, to field capacity, the value at which the soil can no longer gravity drain. For comparison, the equivalent value in MODFLOW, would be 0.0. To obtain a MODFLOW like simulation the user would specify:

SINGLE_UNSAT_SAT 0.0

When using the RE solution to couple to saturated groundwater, soil moistures are calculated along the entire soil column to the free groundwater surface, regardless of the depth of the groundwater. This provides a variable soil moisture along the entire soil column and represents an advantage of the RE method.

More detail on the coupling of GAR and the saturated groundwater model can be found in these documents:

JD31 TN
Soil Moisture Accounting