Techniques and Methods 6A19
U.S. GEOLOGICAL SURVEY
Techniques and Methods 6A19
Vertical flow through a homogeneous unsaturated zone can be approximated with kinematic waves (Colbeck, 1972; Smith, 1983; Smith and Hebbert, 1983; Charbeneau, 1984). This approximation is made by simplifying Richards’ equation, which can be written in the vertical dimension as:
, (1)
where
is 
the volumetric water content (volume of water per volume of rock); 

is 
the water flux (volume of water per time per unit area); 

is 
the elevation in the vertical direction (length); 

is 
the hydraulic diffusivity (length squared per time); 

is 
the unsaturated hydraulic conductivity as a function of water content (length per time); 

is 
the ET rate per unit depth (length per time per length); and 

is 
time. 
If equation 1 is simplified to remove the diffusive term, assuming that the vertical flux is only driven by gravitational forces, the vertical flux owing to gravity is:
, (2)
where
is 
shown as negative downward. 
Substituting equation 2 into equation 1, neglecting
(the diffusive term), and assuming evaporation is removed from the soil profile instantaneously and all flow is downward vertical yields:
. (3)
The method of characteristics solution to equation 3 is obtained by rewriting equation 3 so that is taken partially in terms of and :
, (4)
and by introducing the equation of variation (Abbott, 1966):
. (5)
Equations 4 and 5 written in matrix form yields:
. (6)
and must be indeterminate along characteristic lines (Abbott, 1966). This occurs only when the determinant of the coefficient matrix in equation 6 and the determinants of the matrices obtained by substituting the right hand side vector for each column in the coefficient matrix are all zero, as described in equation 7:
, (7a)
, (7b)
. (7c)
Expanding the determinants yields:
, (8a)
, (8b)
, (8c)
where
is 
the characteristic velocity restricted to the downward (positive ) direction (length per time) and decreases according to the ET rate along a characteristic path. 
The characteristic equations 8a, b, and c define the velocity of a wave, the change in water content of the wave with time, and the change in water content with depth behind the wave, respectively. This wave represents a wetting front in the unsaturated zone that was generated by a sudden increase in infiltration at land surface and that is being affected by ET. Without ET, the water content of the wave is constant over the wave profile. However, the inclusion of ET creates a
nonlinear slope in water content over the wave profile
if a nonlinear function for is used, such as the Brooks‑Corey function (Brooks and Corey, 1966).
The derivative is discontinuous over a sharp
wetting front that results from neglecting the diffusive term in equation 1, such that the spatial derivative in equation 3 is not defined in the absence of hydraulic diffusion. An analytic solution for the velocity may still be derived by considering the effects of diffusion and substituting an equivalent sharp wetting front of equal mass (Smith, 1983; Charbeneau, 1984). The solution to equation 1 for a wetting front that considers hydraulic diffusion can be found by integrating over a control volume containing a single wetting front (Charbeneau, 1984):
, (9)
where
and are points above and below the wetting front at distances far enough that , respectively. 
can be neglected because at and such that:
, (10)
where
and are the values of at depths and . 
Integrating over a profile containing a sharp front with equivalent mass gives:
, (11)
where
is 
the depth of the sharp front. 
Combining equations 10 and 11 gives:
, (12)
where
is 
the velocity of a sharp wetting front (length per time); 

is 
the volumetric water content above a depth ; and 

is 
the volumetric water content below a depth . 
An increase in the infiltration rate will cause a wetting front to form, which is represented by a lead wave. A decrease in the infiltration rate will cause a drying front to occur, which is represented by a trailing wave. Thus, waves are used to represent both wetting and drying fronts. Attenuation of a lead wave occurs as a trailing wave of higher velocity overtakes it. When a trailing wave overtakes a lead wave, the water content of the lead wave becomes equal to the water content of the trailing wave. Consequently, this reduces the velocity of the lead wave. Conversely, when a lead wave overtakes a trailing wave or another lead wave of lower velocity, the overtaken wave is removed, and the water content and flux of the uppermost lead wave are maintained, resulting in rewetting. The process of waves intercepting each other is termed a shock, and this process is important to the method of characteristics because wave properties, such as water content and flux, become discontinuous at the instant a shock occurs (Abbott, 1966).
In contrast to a wetting front that stays sharp, trailing waves elongate with time owing to gravity. Consequently, trailing waves must be divided into a series of incremental waves or be represented by a function that describes internal drainage with time. In either case, internal drainage is determined based on the water content and relative location of points along a trailing wave. For example, points A and B shown in figure 2 represent two points along the trailing wave where the velocities may be determined. The velocity at points A and B is calculated from the characteristic solution (equation 8a) and the ratio of their velocities is equal to the ratio of their depths:
, (13)
where
is 
the depth of the trailing wave increment below land surface represented by point A in figure 2, and 

is 
the depth of trailing wave increment below land surface represented by point B in figure 2. 
In contrast to a lead wave that forms a discontinuity in the water content profile, the derivative is continuous over a trailing front. Accordingly, the BrooksCorey unsaturated hydraulic conductivity function can be used to evaluate . The BrooksCorey function for unsaturated hydraulic conductivity can be expressed as:
, (14)
where
is 
the saturated hydraulic conductivity; 

is 
the residual water content; 

is 
the saturated water content; and 

is 
the BrooksCorey exponent. 
The drainable water from an unconfined or watertable aquifer, referred to as the specific yield or , is equal to the porosity minus the specific retention , where is the volume of water per unit volume of rock that is retained when the rock is drained by gravity (Meinzer, 1923, p. 28). Specific yield is used by MODFLOW2005 to calculate the amount of storage within an unconfined aquifer. Thus, continuity between the unsaturated zone and unconfined aquifers in MODFLOW2005 is maintained through the specific yield by approximating . Residual water content generally is smaller than , such that caution should be used when representing unsaturated flow parameters in the UZF1 Package with measured values.
The velocity of the deepest point along a trailing wave is determined by taking the derivative of equation 14 with respect to :
, (15)
where
and are the velocity and water content, respectively, at the deepest point along a trailing wave. 
A relation between the deepest point along a trailing wave and all other points along a trailing wave is determined by substituting the velocity defined by equation 15 into equation 13 and results in:
, (16)
where
is 
the depth of a point on a trailing wave; 

is 
the water content of a point on a trailing wave; and 

and are the water content and depth of the deepest point along a trailing wave, respectively. 
For example, the depths of points along a trailing wave may be generated by choosing water content values between and shown in figure 2 and by using equation 16 to calculate the depth at a particular time following a decrease in the infiltration rate. Thus, equation 15 may be used to calculate the velocity and the corresponding depth of the front of a trailing wave and then equation 16 can be used to determine the drainage profile above the trailing front.
During a pulse of infiltration, a lead wave’s velocity and water content will decay during a subsequent period of less infiltration as trailing waves overcome the lead wave. The spreading of a trailing wave caused by gravity is defined by equations 15 and 16; however, these equations are not applicable after a trailing wave has intercepted a lead wave. Analytic relations can be derived for routing trailing waves during interception with a lead wave but these equations are difficult to solve. A simpler approach is to divide the trailing wave into steps or increments and calculate velocities of each increment based on a finitedifference approximation (Smith, 1983):
, (17)
where
is 
the change in water content between two adjacent locations along a trailing wave. 
The UZF1 Package relies on equations 15 and 16 to route trailing waves until the trailing waves intercept another lead wave, at which point the trailing wave is discretized into multiple waves that are routed according to equation 17.
The finitedifference approximation of the true trailing wave as depicted in figure 2 does not conserve mass perfectly when a trailing wave intercepts a lead wave and vice versa. However, for problems involving rapidly varying infiltration rates, mass is conserved within an error of less than 0.05 percent when each trailing wave is represented by 15 increments, although 10 increments are sufficient for most problems to achieve a mass balance error of less than 0.05 percent. The user specifies the number of increments used to represent a trailing wave with the UZF1 Package input parameter NSTRAIL.
As described, waves are generated by changes in the infiltration rate and will have differing velocities leading to the possibility of shocks. Determination of the timing of shocks is required for modeling because the wave properties (, ) become discontinuous in time and have to be recalculated in order to continue the solution. Thus, the UZF1 Package routinely calculates the shortest time for any two waves to intersect. Waves are routed continuously between wave intersections, while wave properties remain constant. Consequently, the UZF1 Package requires a unique time discretization that is specific to the problem being solved; the UZF1 Package time steps are different than those used in MODFLOW2005.
Evaporation and uptake by roots in the unsaturated zone can cause negative potential gradients, such that diffusive forces can be important. Evaporation can cause water to move upward in the unsaturated zone by drying out the soil at land surface. These conditions cannot be modeled with kinematic waves because diffusive forces are neglected. However, infiltrating water affected by evaporation and uptake by roots can be approximated with kinematic waves during relatively wet conditions because negative potential gradients are less important. This approach assumes that evaporation and uptake by roots can be grouped together as ET, and that they occur as an instantaneous loss of water over a depth interval equal to the depth of root uptake.
Assuming that is independent of , equations 8b and 8c can be integrated to find:
, (18a)
, (18b)
where
is 
the water content at the head of a lead wave above the ET extinction depth after time , and 

is 
the water content at the head of a lead wave after time subject to the ET rate per unit depth . 
ET is specified as a rate (length per time) within the UZF1 Package (input variable PET) and is converted internally to a rate per unit depth by dividing the specified ET rate by the extinction depth (EXTDP).
Equation 18a is used to calculate the water content at the head of a lead wave through time and equation 18b is used to calculate the depth of points along the leadingwave profile between and , where is the water content at land surface corresponding to the infiltration rate , and and are the water contents at arbitrary points along the leadwave profile with being deeper than (fig. 3). The water content of a trailing wave is decreased according to the ET rate and equation 18a. If the water table elevation is above the extinction depth and the ET demand is not met by the unsaturated zone, then ET is removed directly from ground water. ET is removed from ground water using the same method implemented in the MODFLOW ET Package (McDonald and Harbaugh, 1988, p. 101).
The specified infiltration rate is converted to water content in order to apply the characteristic solutions for unsaturated flow. A relation between water content and infiltration rate can be determined assuming the flux is equal to the unsaturated hydraulic conductivity using equations 2 and 14 (Charbeneau, 1984):
, (19a)
, (19b)
where
is 
the water content of a wave generated from infiltration, and 

is 
the infiltration rate (length per time). 
The water content is set to the saturated water content when the infiltration rate specified in the UZF1 Package input file exceeds the saturated vertical hydraulic conductivity. If the specified infiltration rate is greater than the saturated vertical hydraulic conductivity, then the difference is multiplied by the model cell planview area and this volumetric rate of water can be added to a userspecified stream reach or lake by setting the UZF1 Package variable IRUNFLG to greater than zero. Equation 19 ignores the time to ponding and infiltration rates that exceed the saturated vertical hydraulic conductivity (Smith and Parlange, 1978).
An algorithm was written in FORTRAN using equations 12, 14, 15, 16, 17, 18 and 19 to simulate vertical unsaturated flow and ET. This algorithm keeps track of the location, water content, velocity, and flux of all waves in the unsaturated zone through time and the interaction of waves overtaking one another. Additionally, changes in the thickness of the unsaturated zone owing to changes in the elevation of the water table are considered with respect to waves reaching the water table and contributing to groundwater storage.
Flow from the unsaturated zone as recharge to the water table in an unconfined aquifer is subtracted from the righthand side of the system of equations that are solved by MODFLOW2005:
, (20)
where
is 
a matrix containing the coefficients of the conductance equations (HCOF array) that are solved by MODFLOW2005; 

is 
a onedimensional vector containing the groundwater heads that are solved by MODFLOW2005; 

is 
a onedimensional vector containing all known terms in the conductance equations that are not multiplied by unknown head values (RHS array); and 

is 
the volumetric rate (volume per time) of recharge computed in a given model cell from the UZF1 Package. 
Recharge that is simulated using the UZF1 Package is dependent on the water table elevation (groundwater head). Thus, the UZF1 Package is coupled to MODFLOW2005 within the iteration loop. However, because groundwater recharge is dependent on the amount of storage and the downward flux rate in the unsaturated zone, there is no general equation that defines the relation between recharge and groundwater head, such as the conductance equation in the StreamflowRouting (SFR1) Package (Prudic and others, 2004). Recharge from the UZF1 Package does not affect the HCOF array.
Based on the test simulations presented in this report, the UZF1 Package does not cause additional instabilities or slow convergence of MODFLOW2005 when ground water is not discharging to land surface. Groundwater levels that rise above land surface add an additional head dependency when using the UZF1 Package owing to the calculation of groundwater discharge to land surface, and consequently, convergence may be slowed significantly when ground water is discharging to land surface from many model cells. MODFLOW2005 simulations using the UZF1 Package may increase simulation time compared to using the Recharge Package owing to the additional calculations required to simulate flow through the unsaturated zone.
An infiltration rate at land surface (length per time) is specified for each vertical column in the active model domain (appendix 1). The infiltration rate is multiplied by the planview area of the cell to produce a volumetric flow rate (volume per time). Three options are used in the UZF1 Package to determine where flow in the unsaturated zone recharges ground water (UZF1 variable NUZTOP). The three options are those used in the Recharge Package (McDonald and Harbaugh, 1988, p. 76; Harbaugh and others, 2000, p. 6768). The first option specifies that all recharge is to the top model layer only; the second option allows for recharge to be specified for a particular layer in each vertical column; and the third option allows for recharge to the uppermost active cell in a vertical column. When the last option is used in the UZF1 Package, the program will continue to route water through the unsaturated zone until it reaches the water table of the uppermost active cell. However, the hydraulic properties of the unsaturated zone will remain constant even if the unsaturated zone spans more than one cell in the vertical direction. If there are no active cells in a column beneath land surface, then infiltration and recharge do not occur and a warning statement is printed to the output file. If a cell becomes active during a simulation and the corresponding element in the UZF1 array IUZFBND is nonzero, then infiltration and recharge will start to occur.
As stated previously, mass balance is achieved by approximating the residual water content variable with the specific retention computed as the saturated water content minus the specific yield. Thus, if the water table rises into the unsaturated zone, the aquifer will yield an amount of water equal to the water table rise multiplied by the difference between the water content and the specific retention. Alternately, if the water table falls due to groundwater discharge or pumping at some other location, then the water content in the unsaturated zone between the old and new water table is set equal to the specific retention.
Recharge occurring during a rising water table over a time step is calculated by first routing waves through the unsaturated zone using the base of the unsaturated zone as the watertable elevation at the end of the previous time step. Flow across the base is accumulated for the time step and the volume of water in the unsaturated zone through which the water table rose over the time step is added to the volume that flowed past the base. This sum is divided by the time step to obtain a volumetric recharge rate. Once the flow equation has been solved for the time step, a new base of the unsaturated zone in each active cell is computed. Recharge occurring during a declining water table is computed based on the flux of water from the unsaturated zone during the time step using the base of the unsaturated zone predicted at the previous iteration or time step.
The UZF1 Package also allows for groundwater discharge to the land surface whenever the water table in a cell is higher than the landsurface elevation. The volumetric rate of groundwater discharge to land surface is calculated on the basis of the following equation:
, (21)
where
is 
the watertable elevation; 

is 
the landsurface elevation; 

is 
the planview area of the model cell (equal to the column length times the row length of the model grid); 

is 
the vertical hydraulic conductivity of the model cell; and 

is 
the thickness of the model cell. 
The value that is subtracted from the HCOF array in MODFLOW2005 (corresponds to the A matrix in equation 20) is equal to . The value subtracted from the RHS array in MODFLOW2005 (corresponds to the matrix in equation 20) is equal to .
The UZF1 Package includes an option to specify a twodimensional array of values (UZF1 variable IRUNBND) that identify where groundwater discharge to land surface will be added to a stream or lake if the SFR2 and LAK Packages are active. Additionally, the UZF1 Package variable IRUNDFLG must be greater than 1 for water to be added to streams or lakes. Groundwater discharge to land surface and infiltration in excess of the saturated vertical hydraulic conductivity are added instantaneously to specified streams and lakes. The SFR2 Package is documented by Niswonger and Prudic (2006) and the LAK3 Package is documented by Merritt and Konikow (2000).
A water budget for the unsaturated zone is tracked independently of the groundwater budget within MODFLOW‑2005. The mass balance for the unsaturated zone accounts for all water between the land surface and the water table. The mass balance for the saturated zone is accounted for in the groundwater budget. When the water table rises above the top of a model cell, then the unsaturatedzone mass balance ceases to be calculated for that cell, and unsaturated storage and the change in unsaturatedzone storage for that cell are set to zero.
There are no procedural limitations regarding stress period and time step lengths associated with the UZF1 Package. However, because the UZF1 Package delays recharge to the water table, changes in stress do not necessarily correspond to the beginning of a stress period. Consequently, caution should be used when simulating time steps that increase during a stress period (specifying the variable TSMULT>1.0 in the data input for the Discretization file; Harbaugh, 2005) because water percolating through the unsaturated zone may reach the water table at the end of a stress period, when the time step length is a maximum.
Steady state can be specified for the initial stress period of a simulation when using the UZF1 Package. When steady state is specified for the initial stress period, the groundwater recharge rate computed in UZF1 is set equal to the infiltration rate and ET is not removed from the unsaturated zone. The water content is constant within each model cell during steady state. The specified ET rate is subtracted from the infiltration rate to calculate a steadystate recharge rate. ET can also occur from ground water when the water table is above the ET extinction depth. A uniform water content profile is computed for each active cell, which is used automatically as the initial conditions for a following transient stress period.
Steadystate water content profiles may be printed to a file for specified locations. Initial water content values are not required when the initial stress period is specified as steady state. Thus, initial watercontent values must not be specified within the UZF1 input file when the first stress period is steady state. Another option for developing initial water content profiles is to run the model for several years while repeating a cycle of infiltration rates. This step may be carried out subsequent to a steadystate simulation.
The kinematicwave approximation does not account for capillaryinduced infiltration at the onset of wetting when negative pressure gradients are large relative to gravity potential gradients. However, infiltration approaches the unsaturated hydraulic conductivity when the wetting front has advanced sufficiently far beneath the surface (Mein and Larson, 1973; Childs and Bybordi, 1969). Negative pressure gradients also may result in lateral and vertical redistribution in the unsaturated zone. Because negative pressure gradients are ignored, UZF1 likely under predicts the advancement of the wetting front at early times.
Because capillaryinduced infiltration is neglected, there is no consideration for the period prior to ponding when the infiltration should be represented by a constant flux equal to the rainfall rate (Smith and Parlange, 1978), rather than the constant water content assigned at the boundary (equation 19). Thus for greater rainfall or snowmelt rates, the model will underestimate infiltration during the onset of infiltration. UZF1 does not simulate a capillary fringe above the water table and it assumes that ground water is instantaneously released from or taken into storage.
The initial version of this package (UZF1) does not simulate unsaturated flow through multiple layers of varying hydraulic properties. The model can simulate unsaturated flow through inactive cells to deeper active model cells. However, the hydraulic properties of an unsaturated column between the land surface and the uppermost active cell are uniform during the simulation and are always equal to their values at the onset of the simulation.