Infiltration:Richardsâ€™ Equation

Detailed modeling of the soil water profile in the vadose, or unsaturated zone, is a key addition in the GSSHA model. It is this dynamic area that controls the flux of water between the surface and groundwater and partitions rainfall into infiltration, runoff, groundwater recharge and ET. The most rigorous way to model these complicated and integrated phenomena is to use Richards’ equation in the solution of the problem. While many have described RE as being infeasible for use in hydrologic predictions the equation has been successfully used in field scale and watershed scale models of soil moisture and runoff (Lappala et al., 1987; Hutson and Waggenet, 1989; Dawes and Hatton, 1993; Refsgaard and Storm, 1995; and others). Simpler methods, such as GA and GAR, which are approximations of RE do not provide detailed soil moisture profiles or simulate the movement of water from the groundwater to the unsaturated zone. Accurate representation of layered soils or soils with a water table difficult is also difficult with these approximations (Short et al., 1995). To simulate infiltration with Richards equation use the INF_RICHARDS option in the project file.

Processes in the unsaturated zone important to surface water hydrology, infiltration, ET, and groundwater recharge, are largely oriented in the vertical direction (Refsgaard and Storm, 1995). GSSHA solves the one-dimensional (vertical direction) head-based formulation of Richards’ equation

(14)

where:

C is the specific moisture capacity,

is the soil capillary head (cm),

z is the vertical coordinate (downward positive) (cm),

t is time (hr),

K (

) is the effective hydraulic conductivity (cm),

and W is a flux term added for sources and sinks (cm/hr).

The head-based formulation allows solution of Richards’ equation in both saturated and unsaturated conditions (Haverkamp et al. 1977). The head-based approach has been successfully implemented in the field scale unsaturated flow models VS2D (Lappala et al., 1987), LEACHM (Hutson and Wagenet, 1989), and SWAP (van Dam and Feddes, 2000) and in the surface water hydrology models TOPOG_IRM (Dawes and Hatton, 1993) and MIKE SHE (Refsgaard and Storm, 1995). With a head-based formulation, rainfall, evaporation and groundwater recharge can all be modeled without a change in variable. Known mass balance problems associated with solution of the head-based formula (Celia et al. 1987; Ross, 1990; Pan and Wierenga, 1995; Celia et al. 1990) have largely been eliminated by the development of new solution techniques (Kirkland et al., 1992; Rathfelder and Abriola, 1994; Pan and Wierenga, 1995).

In GSSHA, RE is solved using an implicit finite difference approach, which maps directly to the overland flow finite difference grid. The soil column below each overland flow cell is sub-divided into multiple unsaturated zone cells, Figure 9.

Figure 9 - GSSHA representation of the unsaturated zone

A 1-D approximation of the unsaturated zone imposes some limitations on the application of the model. Lateral flow in the unsaturated zone may be significant for a perched water table that does not extend to the soil surface. Lateral flow of unsaturated water may also occur in the absence of a perched water table (Zaslavsky and Sinai, 1981). Lateral flow in the unsaturated zone under either saturated or unsaturated conditions is most significant for steep slopes. Solution of the multi-dimensional RE is difficult and the additional computational burden is not commonly justified if the purpose of the model is to simulate surface water hydrology. Aspect ratio problems are also avoided because of the 1-D formulation employed. Other conditions in the unsaturated zone that are not explicitly simulated are macro-pores and hysteresis.

7.1.1 Discretization

The soil column must be sub-divided into cells to numerically solve Richards’ equation. The discretization used in the GSSHA model is shown in Figure 9. The soil column is defined in three layers, namely A, B, and C horizons. The soil properties and discretization can vary in each layer. Any size cell can be used in the discretization, however, if very large cells are used the solution may be poor. As discussed in Downer (2002a) and Downer and Ogden (2002b), cells on the order of 1 cm are generally needed in the top 10 cm of soil to accurately model the infiltration process. If larger cells are used in dry soils a large volume of water must be added to the cell before significant infiltration can occur because the hydraulic conductivity will be very low. This can result in significant underestimation of infiltration, as discussed in Downer (2002a) and Downer and Ogden (2002b). The grid size for each soil in the SOIL_TYPE_MAP in each soil layer is specified in the SOIL_LAYER_INPUT_FILE or in the Mapping Table.

The solution scheme employed is implicit, central difference in space and forward difference in time, and is thus second order accurate in space, first order accurate in time. For j being the cell numbering scheme, increasing downward, such that j is the cell of interest, j-1 is the cell above, and j+1 is the cell below, and n represents the time level, Figure 9. The following discretization is used for non-boundary cells

(15)

where:

C = water capacity,

= pressure head in the cell (cm),

K = hydraulic conductivity in the cell (cm/hr),

W = source term (cm/hr),

Δt = time step (hr), and

Δz = grid size (cm).

Values for variables represented between cells, or at the cell interface, are represented by j-1/2 and j+1/2 for the interfaces above and below, respectively. Inter-nodal distances, Δz_j-1\2, Δz_j+1\2, are defined as the distance between the center of cells j and j-1 and j and j+1, respectively, computed as:

(16)

Inter-cell hydraulic conductivities may be calculated two ways, an arithmetic weighting of the values or a geometric average. The inter-cell hydraulic conductivity weighting method is selected with the RICHARDS_K_OPTION project card with the argument being either ARITHMETIC or GEOMETRIC; the default is ARITHMETIC. In using the arithmetic weighting, a weighting factor, , is used to determine how much weight is placed on upper and lower cells. The value of is specified with the RICHARDS_WEIGHT project card; the default value is 1.0. If α = 1, then only the value from the j-1 cell is used; this is commonly referred to as backwards difference, or upwinding. If α = 0 only the j+1 value is used. This is commonly referred to as forward difference, or downwinding. If α is 0.5 then equal weight is applied to both the j-1 and j+1 cells. Lappala (1981) recommends that backward difference be used for modeling the infiltration process. The geometric average is best when there are large changes in the hydraulic conductivity between cells and is generally applicable in all situations.

7.1.2 Non-linear Coefficients

In Richards’ equation, hydraulic conductivity and water capacity values are dependent on the water content of the cell. The method to describe these relationships must be specified with the RICHARDS_C_OPTION project card. The possible arguments are BROOKS for the Brooks and Corey method (1964) and HAVERCAMP for the Havercamp method (Havercamp et al., 1977). The Brooks and Corey method (1964) as extended by Hutson and Cass (1987) may be used to estimate relative hydraulic conductivity from the soil pressure head as:

if < _b (17)

and if > _b K_r = 1.0 if

where: K_r is the relative hydraulic conductivity, _b is the air entry or bubbling pressure, and λ is the pore-size distribution index, which is the inverse of the ratio of the length of flow path through the soil matrix to the straight line length.

With the Havercamp method (Havercamp et al., 1977) as modified by Lappala et al. (1987) the relative hydraulic conductivity is calculated as:

(18)

where: A and B are parameters fitted to laboratory determinations of hydraulic conductivity at different soil pressure heads.

The water capacity C is also dependent on the pressure head . Water capacity is defined as the change in moisture with respect to head . For Brooks and Corey if < b the relationship between moisture content and head is expressed as

(19)

where: , θ is the moisture content, θ_r is the residual saturation, and θs is the saturated moisture content, or porosity. Rearranging and taking the derivative with respect to head yields

(20)

In the original formulation of Brooks and Corey, if > b, θ = s and C=0. This formulation permits changes in pressure head without a resulting change in water content for pressures greater than the bubbling pressure. Hutson and Cass (1987) extended the Brooks and Corey equation into the wet region where > b using the following formulation

(21)

where: , and θ_c is the critical moisture content, the moisture content where the curve changes from the standard Brooks and Corey curve to the modified curve of Hutson and Cass. The critical moisture content is defined as with the corresponding critical pressure head:

(22)

The moisture content can then be represented as

(23)

And, the water capacity is found by taking the derivative of the above equation.

(24)

For the Havercamp equations the relationship between the water content and the pressure head is

(25)

where: α and β are fitted parameters. According to Lappala (1981) the form of the equation can be expressed as:

(26)

Differentiating with respect to pressure yields the water capacity

(27)

The relationship between soil moisture and suction head as represented by the different methods is shown in Figure 10.

Figure 10 – Water retention curves (BC – Brooks and Corey).

7.1.3 Evapo-transpiration Source Term

During long-term simulations potential evapo-transpiration (PET) is calculated by either the Penman Monteith (Monteith, 1965; 1981) or Deardorff (1977; 1978) equations and is applied to each soil column below the overland flow plane. Any water ponded on the surface of the overland flow cell is reduced up to amount of the PET. Any remaining PET demand is applied to cells in the unsaturated zone down to the specified root depth (Figure 9). The actual evapo-transpiration (AET) is distributed over the cells in the specified root zone in proportion to the size of each cell. The AET is computed from the PET by adjusting the PET for the soil moisture in each cell. AET depends on the soil moisture, hydraulic properties of the soil and plant characteristics. At water contents at or above the field capacity (θ_fc) there is no stress on plants and AET is equal to PET (Shuttelworth, 1993; Dingman, 1994). The field capacity is the water content at which the suction pressure prevents gravity drainage of the soil. Field capacity is not an actually physically measurable quantity and there are many descriptions of the conditions that correspond to the concept of field capacity. At water contents below the wilting point, θ_wq, plant transpiration ceases and plants wilt (Dingman, 1994), and AET is zero. At intermediate points AET will depend on PET and the water content. Many relationships to relate AET to PET have been suggested (Dyck, 1983). In GSSHA AET is calculated from PET for water contents greater than the wilting point using the relationship

(28)

The appropriate equation depends on vegetation, climate, and soil type. If the exponent, P, is 1 the relationship is linear, above 1 the curve is convex, and below 1 concave. Currently P is set to 1.0. Future versions will allow the power to be specified by the user to reflect the local conditions. The AET for each cell is added to the source term, W, for that cell in the RE solution.

7.1.4 Upper Boundary Condition

The upper boundary condition varies depending on the condition of the top cell: specified flux for no surface ponding, and specified pressure (head) when infiltration excess results in surface ponding. The first cell in the column, j=0 (Figure 9), is located above the ground surface, and a pressure is always specified for this cell. For a flux boundary condition, the pressure in the top cell is zero and the flux is added to cells via the source term, W. For a head boundary the pressure in the top cell is equal to the depth of ponded water. The typical sequence of events is described below.

At the beginning of simulations and at times between rainfall events there is typically no ponded water on the overland flow plane. The upper boundary condition for the unsaturated zone is a negative flux, equal to the AET of the top cell, cell 1 in Figure 9. When rainfall, runoff or some other source of water is added to the overland flow cell a flux boundary condition is initially specified for the RE solution. A flux equal to the depth of ponded water divided by the current time step is added to the source term, W, of the first non-boundary cell, cell 1 in Figure 9. New heads are calculated with this assumed flux. These heads are used to compute the inter-cell fluxes. The computed flux in cell 1 is compared to the source term. If the calculated flux is less than the specified (assumed) flux, then all the water cannot infiltrate into the soil column during the current time step. In this case, water ponds on the soil surface and the upper boundary is changed to a specified head, the head being the depth of ponded water. Heads are re-computed using the head boundary. To save computation time, the upper boundary condition remains a specified head until the overland flow cell becomes dry. At this point the boundary condition changes back to a specified flux, and remains a specified flux until infiltration capacity of the soil column is again exceeded. Water that enters the top cell is infiltration. Any water that does not infiltrate can become runoff or direct ET (DET).

As discussed in Downer (2002a) and Downer and Ogden (2002b), the assignment of the hydraulic conductivity at the ground surface, K_1/2 between j=0 and j=1 (Figure 9), has important implications in determining infiltration, infiltration excess, and ET. Three methods, specified with the RICHARDS_UPPER_OPTION project card and the appropriate argument, can be used to compute the hydraulic conductivity at the ground surface. The default option, NORMAL, is to use the cell-centered value of hydraulic conductivity in cell 1 (Figure 9). The cell-centered value of hydraulic conductivity is always used for a flux boundary condition. If the soil surface boundary condition is a head, hydraulic conductivity at the soil surface may be assumed to be the saturation value, selected by using the GREEN_AMPT argument, or an average of the saturation value and the cell-centered value, selected with the AVERAGE argument. The assumption is that there is always a very thin layer of saturated material at the soil surface any time water ponds. Either method results in increased infiltration compared to using the cell-centered value of hydraulic conductivity. As discussed in Downer (2002a) and Downer and Ogden (2002b), testing at two watersheds indicated either alternative method allows the use of larger cell sizes in the unsaturated zone without seriously affecting calculated hydrologic fluxes.

7.1.5 Lower Boundary Condition

Three different lower boundary conditions can be specified. When the groundwater table is far from the ground surface the lower boundary condition is a zero head gradient. Water entering the N-1 cell in the soil column exits at the incoming rate (Figure 9). This boundary condition is valid when the water table is so deep that its effect on processes in the upper soil column is negligible, and is the default option when a WATER_TABLE is not specified by the user. The lower boundary can also be a fixed water table, specified with the WATER_TABLE project card that specifies the name of a map containing starting groundwater elevations. When a fixed water table is simulated the top of the last cell, N, is just at the surface of the saturated groundwater (Figure 9). The pressure at the top of the cell is zero, and the pressure at cell’s center is positive, calculated as 0.5ΔZ_N (Figure 9). This groundwater boundary fluctuates when saturated groundwater is simulated. For a moving water table, the size of the last non-boundary cell, N-1, and the number of cells, N, changes as the water table rises and falls (Figure 9). A moving water table is specified by using the GW_SIMULATION card in the project file, and then supplying the required saturated groundwater inputs.

The unsaturated zone does not include the saturated zone. The flux between the saturated and unsaturated zones is calculated as part of the overall unsaturated solution. This separation of the saturated and unsaturated zones helps in determining mass balance errors and maintaining mass balance in both the saturated and unsaturated zones because water is either in one zone or the other. When water crosses the boundary between the saturated and unsaturated zones it is removed from one zone and placed in the other. Mass balance errors can occur for a variety of reasons including: model formulation, spatial and temporal discretization used, solution technique, and abrupt changes in material properties. In GSSHA the mass balance for each compartment, overland, stream, unsaturated, saturated, are calculated independently and integrated to compute an overall mass balance. Mass balance errors in any compartment may be controlled by reducing time steps or increasing stability criteria, as needed. The methods used to link the saturated and unsaturated zones are further described in later sections.

7.1.6 Solution

Richards’ equation is highly non-linear because both the water capacity and the hydraulic conductivity depend on the pressure, or water content, of the soil. Numerical solution of RE requires some type of linearization. In GSSHA the RE is linearized by making the water capacity and inter-cell hydraulic conductivity constant during each time step. With flux updating of the heads, as described by Kirkland et al. (1992), this provides a stable, accurate, and mass conserving solution for most conditions.

For N cells including the upper and lower boundary cells, N-2 equations are needed. The well-known Thomas algorithm (Thomas, 1949) is efficient for solution of the resulting tri-diagonal matrix. After solving for heads in each cell, flux updating (Kirkland et al., 1992) synchronizes the heads and soil moistures, and improves the mass balance. Fluxes (cm hr-1) across the top face of each cell, f_j-1/2, are computed as:

(29)

These are used to compute the change in water content of each cell as

(30)

The new water content is

(31)

Once the water content in each cell is updated, pressure heads in each cell are calculated based on the head-moisture content relationship and used in the next time step. Kirkland (1991) found it necessary to restrict flux updating to cells that were both unsaturated and not immediately adjacent to saturated cells. Testing by the authors confirms this needed restriction. The reason for this requirement is that in saturated cells, head changes can occur without a change in water content. This deficiency in the formulation can cause errors in the solution and the mass balance. In this case, iterating on K_r and C can produce substantial improvements in accuracy and mass balance.

The maximum number of iterations is specified by the user using the RICHARDS_ITERMAX project card, the default being one. While iterating on the hydraulic conductivity and water capacity will almost always improve the solution and mass balance to some extent, iterating is particularly useful anytime there are saturated cells in the unsaturated zone. Cells can become saturated when large changes in material properties exist between adjacent layers, when the upper boundary condition is a head, and when the saturated groundwater is rapidly rising. Experience indicates that when the top boundary is the head condition, the accuracy and mass balance are always improved by iterating on the non-linear coefficients. When the top boundary condition is a head the maximum number of iterations changes from the default to five, unless the user has specified more iterations.

Picard iterations (Celia et al., 1990) are used anytime the user specifies iterations or the upper boundary condition is a head. Heads at the n+1 time level are first calculated using values of K_r and C from the n time level. Water contents are calculated based on the fluxes, and heads are updated based on the water contents. K_r and C are then calculated based on the updated values of head and water content. The water contents and heads are set back to the n time level values and K_rⁿ⁺¹ and Cⁿ⁺¹ are used to compute values of head at the n+1 time level. The procedure is repeated until the convergence criterion is met or the maximum number of iterations is reached. The convergence criterion is applied independently to each soil column. The convergence criterion is the difference in head between iterations. The error is adjusted for the head in each cell because the error is most important for wet cells, where a small absolute error in head can cause a large error in the solution. The error for each cell is expressed as

(32)

where k expresses the iteration number, which should not be confused with n, which represents the time level. Iterations for a particular soil column cease once the maximum change in head is less than 1 mm.

7.1.7 Time Step Limitation

Because the discretization is implicit there is no inherent stability criteria for the scheme. However, a time step limitation is desirable to keep the scheme accurate and mass conserving. The time step limitation in GSSHA is adapted from Belmans (1983). The time step is limited such that a maximum change in water content, Δθ_allow, is not exceeded. If the maximum change in water content in any cell, Δθ_max, exceeds Δθ_allow the time step is reduced so that an exceedance during the next time step is not likely. The following limitation is used

(33)

The maximum allowable change in water content is specified by the user, with a suggested range of 0.002 to 0.03 (Belmans, 1983). Smaller limitations will result in longer simulation times. In GSSHA each soil column has its own time step limitation computation. This can greatly increase the speed of the model when rapid changes in water content occur in only a few soil columns. In this case, Δt for these few soil columns may be very small while Δt for the bulk of the soil columns in the watershed is still very large in comparison. The time step is limited by specifying the maximum water content change allowable with the RICHARDS_DTHETA_MAX project card; the default value is 0.025.