Techniques and Methods 6-A19

Techniques and Methods 6-A19

Back to Table of Contents

Test Simulations

Two hypothetical test simulations are used to illustrate the capabilities of the UZF1 Package. The first is an application of the model to a watershed that encompasses an area of 27 km2. This simulation was designed to demonstrate all the capabilities of UZF1. However, because the first test simulation involved relatively large data input and model output files that are not easily printed, a second simpler test simulation is described that can be used as a guide for data input and model output structure and to verify proper installation of the UZF1 Package in MODFLOW-2005. Detailed results are described for the first test simulation in this document, but the data input and model output are not included as part of the printed document. However, data input and model output for the second test simulation are included in appendix 2 along with a summary of the unsaturated-zone and ground-water budgets in the discussion of the second test simulation. Both test simulations and this report are available at

Test Simulation 1: Sagehen Watershed

The first test simulation is a conceptual yet realistic simulation of the Sagehen watershed (fig. 5), a USGS Hydrologic Benchmark Network Basin located on the eastern slope of the northern Sierra Nevada, near Truckee, California (Mast and Clow, 2000). The basin drains an area of 27 km2 and ranges in altitude from 1,926 to 2,663 m above mean sea level. The areally averaged annual precipitation is about 970 mm and the annual hydrograph is dominated by snowmelt. The Sagehen watershed comprises volcanic rocks overlain by a veneer of alluvium.

The test simulation consisted of a single model layer with 73 rows and 81 columns in which all cells had a constant width and length equal to 90 m (fig. 5). The top elevations of model cells were set equal to estimated land-surface altitudes for each cell; the elevation of the bottom of the aquifer (bottom of model cells) ranged from 1,830 to 2,480 m above mean sea level. The simulation period was 1-year beginning on December 1 and ending November 30. The simulation period was divided into 365 1-day stress periods. The first stress period was steady state, and subsequent stress periods were transient. One time step was simulated for each stress period. No-flow conditions were simulated across the bottom and sides of the model except for three cells beneath and adjacent to the stream at the basin outlet, which were specified as constant head cells.

Several packages and optional files were used in the simulation of the Sagehen watershed (table 2). The Layer-Property Flow (LPF) Package (Harbaugh and others, 2000, p. 59) was used and both horizontal and vertical hydraulic conductivity values were specified. Because one layer was specified in the model, the vertical hydraulic conductivity of the aquifer only affected ground-water seepage to land surface. All cells were specified as convertible (unconfined). The horizontal hydraulic conductivity ranged from 0.005 m/d on the ridges to 0.3 m/d in the valleys (table 3). A lower hydraulic conductivity was specified on the ridges where volcanic rocks are near the surface (fig. 5). The hydraulic conductivity within each cell was assumed isotropic. Specific storage was set to 5 ×10-7 m-1, and the specific yield was specified as 0.05 on the ridges and 0.25 in the valleys near streams.

The vertical hydraulic conductivity of the unsaturated zone had the same values as the horizontal and vertical hydraulic conductivity used in the LPF Package (table 3). A constant Brooks-Corey exponent of 4 and a constant saturated water content of 0.3 were assigned to all cells.

The streams were represented by 15 segments made up of 201 reaches using the Streamflow-Routing (SFR2) Package (Niswonger and Prudic, 2006) and included the simulation of unsaturated flow beneath streams. All stream reaches were set to a constant width of 3 m (ICALC=1) and ranged in length from 30 to 114 m. The streambed and unsaturated zone hydraulic conductivity were equal to 0.3 m/d (table 3). Unsaturated-zone properties were constant among all stream reaches. Saturated water content beneath the stream was 0.30. The Brooks-Corey exponent was 3.5.

The distribution of hydraulic conductivity was created initially on the basis of the surface geology. The distribution of hydraulic conductivity was adjusted during a steady-state calibration that assumed that the spatially varying ground-water recharge was proportional to the distribution of mean annual precipitation. Precipitation was distributed according to 860 mm below an altitude of 2,100 m and increased linearly above 2,100 m to a maximum of 1,140 mm at an altitude of 2,600 m. The range in ground-water recharge for the steady-state simulation was determined by approximating the mean daily discharge at the outlet of the Sagehen watershed in early December. The steady-state simulation resulted in a calculated water table that was as much as 80 m below land surface along the ridges and at, or slightly above, land surface in the lowest parts of the valleys next to streams (fig. 6). Although no observation wells have been drilled on the ridges in the watershed, the maximum depth to ground water was based on depths measured in wells elsewhere in the northern Sierra Nevada.

The initial water content of the unsaturated zone was calculated on the basis of the simulated steady-state recharge rate, which equaled the infiltration rate at land surface. An ET extinction depth of 2.5 m was assigned to all cells (table 3) and ET from the upper 2.5 m was simulated on the basis of the assigned ET demand rate. Water was first removed from the unsaturated zone and then from the saturated zone only when the water table was within 2.5 m of land surface. The water content for ET extinction was set to a constant of 0.1002. Infiltration and ET demand rates were specified for each stress period and were estimated assuming infiltration and ET were at a minimum from December 1 to mid-April, infiltration increased to a maximum during May, and the daily ET demand increased to a maximum during June (fig. 7). The specified infiltration rates varied in the model and were three times greater at the tops of the ridges as compared to the valley lowlands (fig. 8). The range in the applied infiltration rates were estimated on the basis of the average difference between the two precipitation gages located in the upper and lower portions of the watershed.

The Preconditioned Conjugate-Gradient (PCG) Package was used to solve for ground-water head throughout the model domain. The head-closure criterion used in the transient simulation was 0.0009 m and the flow-closure criterion was 0.005 m3/d. The absolute mass-balance error for the unsaturated zone ranged from 0 to 0.08 percent for all stress periods and the saturated zone ranged from 0 to 0.22 percent. The cumulative mass-balance error at the end of the 365-day simulation was 0 percent for the unsaturated zone and 0.10 percent for the saturated zone.

Results from the transient simulation of ET and percolation through the unsaturated zone indicate that recharge patterns are a result of the cumulative effects of ET, spatially varying recharge, hydraulic conductivity, topography, and thickness of the unsaturated zone. Model calculated ground-water recharge summed over the entire basin indicates two distinct patterns: a nearly immediate response to peaks in infiltration in the valley lowlands where the unsaturated zone is thin, and a diffuse response that varies slowly over several months in upland areas where the unsaturated zone is relatively thick (fig. 9).

The simulated volumetric ET rate was correlated to the specified ET demand rate. The cumulative volumetric ET rate for the watershed was limited by water availability in the unsaturated zone and the depth to the water table, except during periods of high infiltration. Because the water table generally was deeper than the ET extinction depth for most of the watershed, a larger proportion of ET was removed from the unsaturated zone (fig. 9). However, the quantity of ET loss is sensitive to the extinction depth and extinction water content, such that attention should be placed on estimating the extinction depth and extinction water content to reduce uncertainty in recharge estimates.

Most of the simulated streamflow during the first six months was from ground-water discharge to the stream (fig.  10). Total overland flow, including ground-water discharge to land surface plus saturation excess, exceeded ground-water discharge to streams beginning mid-May until mid-August. After mid-August, total overland flow was diminished and ground-water discharge to streams dominated the hydrograph (fig. 10).

The outlet streamflow at the end of the simulation period was greater by 5,900 m3/d than that at the beginning; whereas overland flow caused by ground-water discharge to land surface (spring flow) was greater by 1,700 m3/d (fig. 10). The greater ground-water discharge at the end of the simulation is consistent with the greater mean recharge of 46,700 m3/d for the one-year simulation compared with a mean recharge of 10,400 m3/d for the steady-state period. The results of the simulation show how a year of well-above-mean recharge may have a prolonged effect on ground-water discharge to land surface and to streams.

Infiltration and ground-water recharge during the simulations resulted in a general increase of ground-water seepage to land surface between the beginning of the simulation and the end of June when ground-water recharge was at its peak (fig. 10). The increase in the area over which ground water discharged to land surface caused the increase in overland runoff (ground-water discharge to land surface plus saturation excess) shown in figure 11.

Water-content profiles beneath individual models cells indicate that the water content follows a seasonal envelope that becomes narrower with depth because of the seasonal pattern of infiltration and ET (fig. 12). The water-content profiles illustrate the utility of the kinematic-wave approximation in computing vertical flow through the unsaturated zone. The wetting front created by many different infiltration rates took longer than a year to reach the water table in areas where the saturated vertical hydraulic conductivity was low and depth to ground water was greater than 20 m below land surface.

The marked difference among infiltration, recharge, and streamflow for the example of the Sagehen watershed suggests that the unsaturated zone could have a strong effect on fluctuations in ground-water levels and streamflow during and after storms. Years in which infiltration rates exceed the mean annual rate may increase ground-water storage even in relatively poorly permeable rocks and result in increased discharge to streams that may extend well beyond the period of infiltration. The UZF1 Package in MODFLOW-2005 provides a method in which daily-averaged values of overland runoff and saturation excess runoff can be simulated in a complex, yet simplified, manner.

Test Simulation 2: Simple Hypothetical Problem

A second smaller test problem was developed for use as a data-input template and testing proper code installation. Abbreviated listings of the data-input and output files for test simulation 2 are included in appendix 2.

Test simulation 2 used the hypothetical stream-aquifer interaction problem described in Prudic and others (2004; p. 13–19), except the UZF1 Package was used to replace the ET and Recharge Packages. The units used in the hypothetical stream-aquifer interaction problem were in English units, and were not converted to SI units for test simulation 2. Thus, the data input and printed output listed in tables 4–8 and in appendix 2 are in English units.

The hypothetical stream-aquifer interaction problem used in this test simulation was developed for an alluvial basin in a semiarid region. The principal aquifer comprises unconsolidated deposits of mostly sand and gravel. Recharge to the aquifer primarily is leakage from the streams that enter the basin from the mountains on the northwest, northeast, and southwest (fig. 13). The main stream in the southern part of the valley is perennial, whereas all other streams in the valley are intermittent with small drainage areas.

Streamflow entering the model domain and diversions from streams were the same as that used by Prudic and others (2004) and remained the same for all stress periods (table 4). Different methods were used to represent the relations between width, depth and flow in the streams (see ICALC in table 4 and ICALC option descriptions in Prudic and others, 2004). Infiltration rates in the UZF1 Package file and pumping rates in the Well Package file were varied during 12 stress periods (table 5) but the distribution of infiltration and pumping among grid cells did not change. The ET demand was specified as constant over the entire modeled area equal to 1.6 ft/yr to a maximum depth of 15 ft below land surface. The LPF Package was used to specify hydraulic properties for the aquifer. The hydraulic conductivity and specific yield were 173 ft/d and 0.2, respectively, in the vicinity of the stream channels and the hydraulic conductivity and specific yield were 35 ft/d and 0.1, respectively, elsewhere in the alluvial basin (fig. 14; table 6).

The model grid was divided into uniformly spaced cells 5,000 ft on each side (fig. 14). The strongly implicit procedure (SIP) was used to solve the flow equation for test simulation 2. The head-closure criterion was 0.0002. The option to add overland runoff to stream segments was specified in the UZF1 Package; however, no ground-water seepage to land surface was simulated in the test. The number of stress periods in the Discretization file was 12. Each stress period was 2.628 ×106 seconds or 30.42 days. The first stress period was steady state and had one time step. This allowed for the calculation of initial water contents for the unsaturated zone on the basis of the steady-state infiltration rate. The remaining 11 stress periods were transient and were divided into 15 time steps that increased sequentially by a factor of 1.1.

Water budgets at the end of selected stress periods for ground water and the unsaturated zone are summarized in tables 7 and 8, respectively. The results indicate that changes in flow and storage in the unsaturated zone affect ground-water recharge and, when combined with stream interaction, ground-water ET, and ground-water pumping, the effect of the unsaturated zone can be important to the overall assessment of ground-water flow. The absolute mass-balance errors were less than 0.05 percent for both the unsaturated-zone and ground-water budgets, and the cumulative mass-balance error was 0 percent for the unsaturated-zone budget and -0.01 percent for the ground-water budget.

Back to Table of Contents

AccessibilityFOIAPrivacyPolicies and Notices

Take Pride in America home page.FirstGov buttonU.S. Department of the Interior | U.S. Geological Survey
Persistent URL:
Page Contact Information: Publications Team
Page Last Modified: Friday, 02-Dec-2016 15:45:19 EST