Scientific Investigations Report 2005-5227
U.S. GEOLOGICAL SURVEY
Scientific Investigations Report 2005-5227
Several numerical computer models have been constructed to represent the ground-water flow system in at least part of the SVRP aquifer. At least four of these computer models have been published since the early 1980s (Bolke and Vaccaro, 1981; CH2M HILL, 1998; Buchanan, 2000; and Golder Associates, Inc., 2004).
Bolke and Vaccaro (1981) constructed a two-dimensional (plan view) computer flow model of the Washington part of the SVRP aquifer from near Post Falls, Idaho, on the east to near Nine Mile Falls, Washington on the west. This part of the SVRP aquifer was termed the Spokane Aquifer. The finite-element model was developed by USGS and simulated ground-water flow by solving a set of linear equations derived by J.V. Tracy (U.S. Geological Survey, written commun., 1977). A two-dimensional model was deemed adequate since vertical flow or gradients in the aquifer were thought to be insignificant.
The model indicated that pumping at the current rate (1977) had little effect on water levels in the Spokane aquifer. The model was used to forecast increased pumping effects on aquifer heads and streamflow. During a one-year simulation, pumping at twice the 1977 rate of 227 ft3/s resulted in calculated water-level declines of about 3 ft. Spokane River streamflow was calculated to decline 150 ft3/s in the summer and about 50 ft3/s during the rest of the year. Flow to the Little Spokane River was modeled to decrease about 10 ft3/s. The increased pumping rate had a more significant effect on the discharge of the Spokane River than on the change in water levels in the aquifer.
The eastern boundary basically was an approximate north-south line near Post Falls, Idaho, that split the “Spokane Aquifer” from the rest of the SVRP aquifer to the east. This boundary was assigned a specified head value determined from an average of heads in nearby monitoring wells. The specified head value was held constant during simulation. The western part of the model between where the Spokane and Little Spokane Rivers exit the model was classified as specified head boundary to simulate constant ground-water outflow. Ground-water inflow was specified and held constant in areas where surrounding drainage areas entered the modeled area. The ground-water inflow estimates for these areas were obtained from Drost and Seitz (1978b). The remaining external boundaries represented the contact between the aquifer and the surrounding bedrock. These boundaries were specified as no-flow boundaries. The bottom of the aquifer also was specified as a no-flow boundary.
Lateral hydraulic conductivity values were estimated from specific-capacity data. Values ranged from about 0.07 ft/s (6,048 ft/d) in the eastern part of the aquifer to about 0.01 ft/s (864 ft/d) in other parts of the aquifer to the west. Values in the Five Mile Prairie area and other areas where bedrock exists in the modeled domain ranged from 0.001 to 0.00001 ft/s (86.4 to 0.86 ft/d). Available data indicated an absence of vertical stratification in the aquifer lithology. Therefore, the lateral and vertical hydraulic conductivity was interpreted as equal (no anisotropy).
Saturated thickness values were determined from the difference between a water-table map of 1977-78 data and the estimated bottom of the unconsolidated sand and gravel deposits. The basal altitude was based mainly on drillers’ logs and on two seismic profiles by Newcomb and others (1953). Transmissivity values were computed by multiplying estimated saturated thickness by lateral hydraulic conductivity. Specific yield initially was estimated by comparing lithologic information with grain size/specific yield tables from Johnson (1967).
Initially, gains and losses of Spokane and Little Spokane Rivers were obtained from previous work by Broom (1951). Broom (1951) divided the rivers into 7 reaches and computed values for water year 1950. During subsequent model analysis, the Spokane and Little Spokane Rivers were divided into 13 gaining and losing reaches based on the May 1977 to April 1978 modeled period. Previous studies reported that the Spokane River was an entirely gaining reach below the bridge at Barker Road; whereas Bolke and Vaccaro (1981) reported 5 losing reaches and 5 gaining reaches below the bridge at Barker Road (fig. 1).
Within each of the 13 reaches, leakage coefficients were assigned a uniform value based on streambed lithology and whether the river was gaining or losing. Relative differences between head in the river and heads in the adjacent aquifer, leakage coefficients, and the streambed area were used in a form of Darcy’s Law to quantify flow to or from the river. Water level altitudes in reservoirs were treated as specified heads for the model simulation.
Ground-water withdrawals were determined from a water-use inventory of all major water purveyors in the study area. Total ground water withdrawal in the modeled area during 1977 was about 164,000 acre-ft (about 227 ft3/s). Based on the percentage of water pumped each month, the total withdrawal was distributed monthly for the transient simulation.
The percentage of water pumped from the aquifer and applied at land surface or discharged to sewers was estimated based on information from purveyors. The estimates were: 50 percent of the water pumped from seasonally used city of Spokane wells was applied to the land surface and 50 percent was discharged to sewers; 10 percent of the water from continuously used city of Spokane wells was applied at land surface and 90 percent was discharged to sewers; 70 percent of the water from other municipalities was applied at land surface and 30 percent was estimated to directly recharge the aquifer through septic systems; 100 percent of the water pumped for agricultural irrigation was applied at land surface; and none of the water pumped for industry was applied at land surface.
Precipitation was set to be uniformly distributed in the modeled area based on data from weather stations at the Spokane International Airport and at Spokane (fig. 3). An average precipitation rate was set at 1.72 in/mo for the steady-state simulation. Actual monthly values were used for transient simulation.
Depth to water in nearly all parts of the modeled area was thought to be too great to allow for transpiration by plants, except near the Spokane and Little Spokane Rivers. Evaporation also was known to occur from the surface-water bodies. These transpiration and evaporation values were considered negligible. However, estimates of evapotranspiration (ET) were needed to calculate the amount of water that reaches the aquifer from precipitation and water applied at land surface.
Potential evapotranspiration or consumptive use by crops was estimated by the Blaney and Criddle (1962) method. This method assumed total availability of water to the crops and was dependent on air temperature and length of growing season. Estimates of actual ET were made by assuming that some or all precipitation falling during the crop-growing season was available to meet the water requirements (actual ET) of the crops. This precipitation (effective rainfall) was about the same as the actual ET. Actual ET or effective rainfall was estimated by a method developed by the U.S. Department of Agriculture (1967). Effective rainfall was calculated on a monthly basis, which provided an actual ET estimate that was used in the model to determine the amount of precipitation reaching the ground-water table. The amount of water that infiltrated to ground water was calculated as the difference between the sum of the precipitation and applied irrigation water and the amount of ET, where ET was the lesser of either the actual ET or the potential ET. An average potential ET was set at 1.31 in. for the steady-state simulation; monthly values of potential ET and actual ET were used for the transient simulation.
Initial estimated transmissivity values were adjusted upwards by a factor of 1.9 and specific yield values were uniformly halved (ranging from 0.1 to 0.2 in the unconsolidated material and less than 0.05 in the Five Mile Prairie area) during model calibration.
Model sensitivity analysis indicated that the model was most sensitive to changes in transmissivity. Therefore, potential errors in transmissivity would have a greater effect on the model-predicted heads and discharges than would errors in boundary flows or specific yields.
CH2M HILL (1998) constructed a 3-dimensional numerical flow model of the Spokane Aquifer (Washington part of the SVRP aquifer) as part of the city of Spokane Wellhead Protection Program. The modeling was designed to represent two scenarios (conditions in September 1994 and April 1995). The software MicroFem (Hemker and van Elburg, 1986) was used for modeling steady-state conditions. MicroFem is a finite-element model capable of simulating horizontal and vertical anisotropy.
The calibrated model was used to estimate ground-water capture zones using particle tracking procedures. Particle tracking was conducted for eight existing and two planned well fields in the city of Spokane. Particle tracking results then were used to develop wellhead protection areas for each of the well fields.
A specified flux boundary was used at the Washington−Idaho state line (eastern extent of the model). Flux was initially estimated from the hydraulic gradient in the area, the cross-sectional area, and the initial estimate of transmissivity. The flux was adjusted manually whenever transmissivity values were altered at the state line.
Specified flux boundaries also were defined where tributary valleys and creeks intersected the aquifer margins and infiltrated the aquifer. The flux at these sites was allowed to vary during model calibration.
Specified heads were used along the model’s downgradient boundaries at Nine Mile Reservoir and Little Spokane River. Specified head values were based on an autumn 1994, ground-water altitude map.
Areas with very low transmissivity, essentially no flow boundaries, were set around Five Mile Prairie and other bedrock highlands. Low transmissivity values also were set in areas where the aquifer was not present in layer 2 or layer 3, but was present in an overlying layer.
The ground-water system was modeled with 3 layers designed to simulate ground-water withdrawals at different depths: Layer 1 represented the upper 100 ft of the saturated zone; Layer 2 represented depths from 100 to 200 ft below the water table; and Layer 3 represented variable depths from 200 ft below the water table to the bottom of the aquifer.
Hydraulic conductivity was simulated with 20 different zones and was allowed to vary in most of the zones during calibration. Hydraulic conductivities of the zones ranged from <20 to 7,000 ft/d in the calibrated model. The vertical anisotropy was set to 10:1 and was fixed during calibration.
The Spokane River was modeled with specified water-surface elevations measured at 23 stations along the river. Equal elevations were specified for river nodes in Nine Mile Reservoir and Upriver Pool based on stage measurements at each site. Stages in areas between measurement sites were adjusted to incorporate riverbed elevation variations shown in Federal Emergency Management Agency maps. The remaining nodes representing the river were assigned elevations based on interpolation between the stations.
Riverbed leakage rates were specified for 16 river reaches based on the locations of stage measurements. Initial riverbed leakage coefficients of Bolke and Vaccaro (1981) were used and varied during calibration. Rates were varied according to the assumption that coefficients would be higher in areas where ground water discharges to the river than in areas where the river recharges the aquifer. The model simulated leakage rates based on the leakage coefficients and the difference between the water table altitude and river stage. The WADI package of MicroFem (an environment where surface water is separated from ground water) also allowed for leakage simulation from a streambed that is separated from the water table by an unsaturated zone. In those cases, the model recognized that the leakage rate was independent of the depth to the water table.
Withdrawal rates in the model were determined from field observations and from information supplied by purveyors. For many purveyor wells, withdrawal estimates were determined from pumping logs for the September 1994 and April 1995 periods. Withdrawal estimates for other wells were based on the well’s capacity and the estimated pumping duration. Pumping primarily was assigned to Layer 1, with some pumping in Layer 2, and no pumping in Layer 3. Withdrawal rates were not adjusted during calibration.
Recharge from precipitation was based on previous work by Bolke and Vaccaro (1981) and Olness (1993). Bolke and Vaccaro (1981) estimated that a uniform rate of 66 ft3/s percolated to the water table after evapotranspiration. Olness (1993) determined that spatial variations of precipitation existed across the modeled area. An initial rate of 66 ft3/s was set across the model area, but the spatial distribution determined by Olness (1993) also was used. Recharge from precipitation was allowed to vary during calibration.
Recharge from land-applied water, which includes the percolation of water used outdoors in urban and suburban areas and irrigation on agricultural lands, was included in the model. Recharge from land-applied water was not allowed to vary during calibration.
Recharge from outdoor water use in urban and suburban areas within the city of Spokane was specified from pumpage estimates and effluent discharge volumes. The analysis indicated that about 45 percent of Spokane’s annual water use consisted of outdoor use. Recharge to the water table was assumed to be about 30 percent of the total outdoor use, resulting in an average recharge of 2.5 in/yr within Spokane city limits.
Recharge rates from outdoor use outside Spokane city limits were estimated from pumpage rates supplied by purveyors and by estimates of population densities. The pumping rates and population densities were compared with those within the city limits. Recharge rates were assigned 2.5 in/yr in areas with the highest densities (about 6,000 persons per square mile), 1.0 in/yr in areas with moderate population densities (about 1,000 persons per square mile), and 0.25 in/yr in areas with the lowest population densities (about 100 persons per square mile). Recharge from urban and suburban outdoor use was set equal to zero elsewhere, including agricultural and undeveloped fallow lands.
Recharge for agricultural lands was set at 2.0 in/yr. This rate was set for areas outside the city limits where irrigation is known to occur.
Recharge from septic systems was assigned to urban and suburban areas outside the city limits. Recharge rates were based on an assumed discharge rate of 75 gallons per capita per day to septic drain lines and population densities used to quantify recharge from land-applied water. The resulting recharge rate was set at 16 ft3/s and was not allowed to vary during calibration.
The model was calibrated using September 1994 conditions, which were considered to represent short-term steady-state conditions. The following parameters were adjusted during calibration: transmissivity; river leakage coefficients; recharge from precipitation; inflow from tributary valleys; and ground-water flow rates across the state line. The objective of the calibration process was to obtain reasonable simulations of ground-water altitudes, ground-water flow directions, and the aquifer water budget including river gains and losses.
The model was “verified” by simulating conditions measured in April 1995. All values for the April 1995 simulation were the same as the September 1994 simulation except for the pumping rates, river stages, ground-water altitudes at the downgradient fixed-head boundaries, areal recharge rates, and recharge from tributary boundaries. Simulated ground-water altitude contours in most areas of the model were within 5 ft of ground-water altitude contours based on September 1994 measurements. The 5-ft maximum residual value equals about 1 percent of the total observed ground-water altitude difference in the model area. Ground-water flow directions were closely simulated by the calibrated model throughout the model area. Simulated river gains and losses generally were in agreement with historical streamflow records and other independent estimates.
John Buchanan, a professor of geology at Eastern Washington University, constructed the first ground-water flow model of the entire SVRP aquifer. The finite-difference, single-layer, steady-state model was designed as a tool for understanding the overall water balance. The numerical simulation utilized MODFLOW code (McDonald and Harbaugh, 1988). Results of the model indicated that (1) less recharge to the aquifer occurred on the Idaho side of the aquifer than previously estimated and (2) the calculated ground-water flux at the Washington−Idaho state line (390 ft3/ s) was about one-half that of earlier estimates made during the sole source aquifer designation process.
The SVRP aquifer is bounded by lakes and hillslopes around the periphery of the Rathdrum Prairie in Idaho, and downgradient by the Columbia River basalts and the Little Spokane River. Constant head nodes were used to represent peripheral lakes and the Little Spokane River to provide appropriate inflow/outflow to and from the bounding sides of the model. Model data for inflows was from Bolke and Vaccaro (1981) and Painter (1991b). Lateral boundaries were defined geologically with the bedrock contacts and were specified as no-flow boundaries.
Estimates of hydraulic conductivity were derived from previous modeling and are based on field data and calibrated model values. Vertical and horizontal hydraulic conductivity were set to be equal, therefore no anisotropy. Hydraulic conductivity values decreased from east to west ranging from about 10 to 50,000 ft/d. The porosity of the aquifer was set at 20 percent throughout the model domain.
Aquifer thickness was estimated using the seismic reflection profiling by CH2M HILL (1998) in Washington and another seismic line surveyed near Twin Lakes on the Rathdrum Prairie. However, the depth to the aquifer base remains unknown throughout much of the Rathdrum Prairie.
The Spokane River was represented as river nodes using the “Rivers” module in MODFLOW. In each river node, the bed elevation, bed thickness, bed conductance, and river stage were specified.
Since there was no comprehensive source of water usage statistics in the Rathdrum Prairie part of the aquifer, water withdrawal through pumping wells was not included in the model.
Recharge due to precipitation was applied to the top surface of each active cell in the model and was estimated at 25 percent of the rainfall volume. Recharge rates ranged from 0.0014 to 0.0023 ft/d. Additional recharge from hillslopes and small basins adjacent to the aquifer were applied to appropriate cells. These values were from Bolke and Vaccaro (1981) and Painter (1991b).
The model was calibrated with measured heads in the aquifer reported by CH2M HILL (1998) and Brian Painter (Buchanan, 2000). Hydraulic conductivity values were repeatedly changed near the “lake nodes” until the calculated heads approached measured heads. However, it was noted that the calibrated model by no means presents a unique solution.
Golder Associates, Inc., constructed computer models of the Little Spokane River and Middle Spokane River watersheds for Spokane County, Washington. The project’s objective was to simulate all major hydrologic processes in the watersheds and was intended for use in planning and management of watershed hydrologic resources. Components in the model included subsurface flow in terms of ground water and unsaturated flow, surface water in terms of overland and river flow, and coupling between the surface water and ground water. The integrated model included the MIKE 11 HD one-dimensional model to simulate surface-water flow and the MIKE SHE finite-difference, three-dimensional ground-water flow model. The modeling effort represented water years 1994–99.
Surface water discharge measurements from the USGS Post Falls gaging station were used to provide a time-varying flow boundary for the Spokane River on the eastern part of the model at Post Falls. A constant head boundary was assigned along the Long Lake section of the Spokane River. Impermeable or no-flow boundaries were assigned to all lateral boundaries for Layer 1 except on the eastern boundary near the Washington−Idaho state line where a time-varying head boundary (constant head for 3-month intervals) was established. An impermeable boundary was assigned to the model base. Internal surface-water boundary conditions included inflows from Hangman (Latah) Creek to the Spokane River and several wastewater discharge locations. Discharge data from a USGS gaging station were used for the boundary at Hangman (Latah) Creek. Monthly discharge data were used for each wastewater discharge point.
Topography and a spatial roughness coefficient map were used to model overland flow in the watersheds. Topography was obtained from USGS Digital Elevation Models (DEM). Roughness coefficients were assigned to each of the 21 land cover classes in the USGS National Land Cover Data (NLCD) grid coverage for the area.
Unsaturated zone data needed for the model included depth, distribution, and hydrologic characteristics of soils within the model domain and drywells. The soils were divided into one of four Natural Resources Conservation Service hydrologic groups. Vertical saturated hydraulic conductivity and moisture retention curves were assigned to each soil group.
Drywells, or infiltration pits, were simulated with a “bypass” function which recharges a specified percentage of ground water to the saturated zone when water content in the unsaturated zone is greater than a minimum value. Recharge volume resulting from a dry well was decreased linearly when water content fell below the minimum water content value and was discontinued when water content fell to a set “stop” value. Dry wells were only assumed to exist in unsewered areas. The number of dry wells in a cell was estimated by the amount of impervious area like roadways or parking lots in a cell. The percentage of recharge above the estimated natural recharge was assigned to each cell based on the estimated number of dry wells in each cell.
The ground-water flow model included two aquifer layers, a low conductivity lens in the Hillyard Trough area, and an impermeable basal boundary. Layer 1 was defined by the extent of the highly permeable, glaciofluvial deposits of sands and gravels. The clay lense divides the unconsolidated glaciofluvial deposits into confined and unconfined components. The clay layer was modeled as a lens in Layer 1. Layer 2 was defined by the extent of Tertiary basalt and the Latah Formation. The impermeable basal boundary was defined by the upper surface of the crystalline basement rock.
Spatial distribution and hydraulic conductivity values estimated by CH2M HILL (1998a, b) and Buchanan (2000) for the SVRP aquifer were used to develop initial values for the model. Conductivity values were adjusted during model calibration. A horizontal to vertical anisotropy ratio of hydraulic conductivity of 3:1 was initially used. These values also were refined during calibration. Hydraulic properties of the Grande Ronde basalt were used to represent the Tertiary unit in Layer 2. Boese and Buchanan (1996) provided estimates of hydraulic conductivity for the Tertiary basalts. A hydraulic conductivity of 0.5 ft/d initially was used to represent the layer, but was adjusted to 15.9 ft/d after model calibration. Whiteman and others (1994) and Boese and Buchanan (1996) provided vertical hydraulic conductivity estimates of the basalts ranging from 0.0005 ft/d to 3.5 ft/d. A vertical hydraulic conductivity of 1.59 ft/d initially was used in the model.
Only the primary lakes and rivers were modeled in the watersheds. The MIKE 11 HD model required data on surface-water channels including altitude, length, blocking structures, cross sections, roughness, and leakage coefficients. The cross sections provided altitude and channel shape for the model. An average range of Mannings roughness values was used depending on the population density, slope, and land cover in the basin. Structures within the Spokane River were modeled as weirs. Estimates of leakage between the Spokane and Little Spokane Rivers and the aquifer developed in previous modeling efforts were used for the initial leakage values in this model where possible, these values were adjusted during calibration.
Ground water withdrawal included purveyor, industrial, commercial, irrigation, and residential uses. Abstraction rates for purveyor, industrial, and commercial wells were obtained from Spokane County. A total of 191 wells were modeled with a total annual withdrawal of 52,663 Mgal/yr. Abstraction rates for exempt wells were modeled at estimated consumptive use rates to account for lawn or agricultural irrigation. Abstraction or consumptive use associated with the exempt wells was modeled in the same cell as where the irrigation occurred.
Precipitation and temperature data were input into the model as daily, 4 km gridded data using a method outlined in Bauer and Vaccaro (1987). This method involves the comparison of point measurements and PRISM (Parameter-elevation Regressions on Independent Slopes Model) data (http://www.ocs.oregonstate.edu/). A degree-day snow melt factor of 2.0 mm snow/day/ºC and a snowmelt threshold of 1.0ºC were used to simulate snowmelt within the model.
Evapotranspiration was calculated using spatial and time varying components of potential evapotranspiration (PET) along with relative evapotranspiration potential of vegetative cover, unsaturated zone characteristics, water surfaces, snow coverage, and saturated zone characteristics. Monthly PET estimates were calculated using the Blaney Criddle (Doorenboos and Pruit, 1977) method. Monthly temperature estimates were obtained from PRISM (http://www.ocs.oregonstate.edu/).
Irrigation estimates were based on previous work from Golder Associates, Inc. (2003). Irrigation was modeled to take place during the months of April through October. Lawn irrigation was estimated at 3.69 ft/yr and agricultural irrigation in the Middle Spokane River watershed was estimated at 1.58 ft/yr. Coverages of irrigated lawns and agricultural lands were provided by Spokane County. The use of wastewater discharge for irrigation purposes (such as for golf courses) also was included as irrigated areas. The source of irrigation water was assumed to be municipal if the area was in current water district boundaries and from a shallow well in the same cell if the area was outside of the water district boundaries.
The model was calibrated using ground water levels available throughout the Spokane Valley and discharge data available from four points along the Spokane River. Data were more limited for the Little Spokane River watershed. Parameters that were adjusted during model calibration included: vertical and horizontal hydraulic conductivity of the aquifer layers; unsaturated zone saturated hydraulic conductivity; unsaturated zone moisture retention and hydraulic conductivity curves; river bed lining leakage coefficients, snow melt degree-day coefficient and melting point temperature; overland run-off parameters; and drainage time constants. Spokane watershed model calibration primarily was done through aquifer property modifications.
For more information about USGS activities in Washington, visit the USGS Washington District home page .
U.S. Department of the Interior |
U.S. Geological Survey
URL: https://pubsdata.usgs.gov/pubs/sir/2005/5227/section7.html
Page Contact Information:
Page Last Modified: Thursday, 01-Dec-2016 18:45:52 EST
