Groundwater: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.