Scientific Investigations Report 2006–5137
U.S. GEOLOGICAL SURVEY
Scientific Investigations Report 2006–5137
The water resources of the mid-Snake region have been studied in investigations of regional and local hydrologic systems of Idaho, including investigations of the Snake River Plain in southern Idaho that were part of the USGS RASA and National Water-Quality Assessment (NAWQA) programs. The RASA and NAWQA programs in southern Idaho have produced an assemblage of geologic, hydrologic, and water-quality information about the region’s ground-water systems, as well as computer simulation of those systems. For the mid-Snake region study, we used this information, augmented with data collected in more recent and ongoing water-quality assessments of the area, to develop the Snake River Nitrate Scenario Simulator (SRNSS). Because similar information resources are available for much of the Snake River Plain, the investigation and tool-building approach developed in this study might be applicable to other priority water-quality management areas of the plain.
We based the configuration of the subregional model used in this study on a regional, three-dimensional ground-water flow model of the eastern Snake River Plain (ESRP model) that was developed and calibrated during the USGS RASA program (Garabedian, 1992). We derived flow boundary conditions and initial estimates of aquifer properties for the new subregional model from the calibrated ESRP model. Because the ESRP model and its use influenced our modeling approach, it is important to note here some of the ESRP modeling techniques described by Garabedian (1992) that we used to develop the subregional model for this study.
The goal of the eastern Snake River Plain modeling was to reproduce known aquifer conditions (for example, head and spring discharge) within reasonable ranges of values, and to examine the hydrologic effects of changes in model input data through sensitivity testing. Initial values of aquifer hydraulic properties, recharge, discharge, and pumpage were estimated for model input based on geologic and hydrologic information. Model output was compared with known aquifer conditions to determine whether the initial estimates were reasonable. The model was tested to determine its sensitivity to changes in transmissivity, storage coefficient, aquifer leakance, recharge, riverbed or spring-outlet conductance, ground-water pumpage, and boundary flux. Input data were adjusted within reasonable ranges to achieve a better fit to known conditions. Adjustments were made to least-known and to most-sensitive model parameters. Simulation results indicated how the aquifer might have responded to past stresses, such as increased pumping and reduced recharge, and how the aquifer might respond to hypothetical stresses in the future.
To model the regional aquifer system, assumptions were made about aquifer properties, hydraulic fluxes, and initial conditions for transient analysis. The three-dimensional finite-difference ground-water-flow model MODFLOW-88 (McDonald and Harbaugh, 1988) was used most extensively. Vertical variations in head within each model layer were assumed to be negligible, and head losses between layers were assumed to be controlled by confining beds near the base of each layer. Therefore, model-simulated heads were an approximate average of heads within that aquifer layer.
According to Garabedian (1992), the regional aquifer system was subdivided vertically into four model layers. The layers were assigned equal thicknesses because differentiating the predominantly basalt aquifer system into distinct geohydrologic units was not possible. Layer 1 represents the upper 200 ft of the aquifer system; layer 2 is the next 300 ft below. Layers 1 and 2 contain Quaternary basalt of the Snake River Group and Tertiary basalt. Layers 3 and 4, are of lesser areal extent, and they are present only where basalt of the Snake River Group and interlayered sedimentary rocks are greater than 500-ft thick. Layer 3 is ≤500 ft in thickness, and is present only across the central part of the plain. Layer 4 ranges in thickness from 0 to about 3,000 ft in the central part of the plain.
Local aquifers perched above the regional aquifer system were not simulated, although recharge to perched aquifers was assumed to reach the regional aquifer. Vertical hydraulic conductivity was assumed to be anisotropic, owing to low-hydraulic-conductivity basalt between permeable flow tops and fine-grained layers within sand and gravel zones. Horizontal hydraulic conductivity was assumed to be isotropic because two-dimensional simulations indicated small differences in modeling results between isotropic and anisotropic conditions.
A steady-state model was constructed using 1980 water-year fluxes (Garabedian, 1992, table 16). The 1980 fluxes were assumed to approximate the average annual flux for 1950–80. During that period, hydrologic conditions were stable relative to conditions from 1880 to 1950. Ground-water recharge was assumed to equal discharge (steady-state conditions) from 1950 to 1980 because irrigation diversions and ground-water levels were assumed to be stable relative to conditions from 1880 to 1950 (Garabedian, 1992, fig. 8 and pl. 4). Declines in ground-water-level due to pumping, climatic variations, and decreased surface-water diversions during that period were generally small, and changes in storage were accordingly small (about 1 percent). These small changes in storage were included as a part of the recharge term by including them in the computation of approximate steady-state fluxes. Recharge from irrigation was assumed to take place directly below surface-water-irrigated areas.
Initial conditions were approximated by the steady-state solution obtained after removing recharge due to surface-water irrigation and discharge by ground-water pumping from the calibrated 1950–80 steady-state model. Surface-water altitudes, used in the model for river-leakage simulations, were corrected for pre-reservoir conditions. The estimated pre-irrigation steady-state condition was a stable initial condition for transient simulation. Therefore, simulated changes in head and discharge were assumed to result from changes in model input fluxes rather than from non-equilibrium initial conditions.
To develop the subregional model used in this report, we converted the calibrated ESRP flow model from Modflow-88 to MODFLOW 2000 (Harbaugh and others, 2000). Initial model runs were made on the calibrated regional model to calculate the quantity of ground-water flow across a newly-defined upgradient subregional model boundary. This boundary was set near the western boundary of Minidoka County at column 25 in the ESRP model (Garabedian, 1992).
We then “clipped” the ESRP RASA model grid to this new upgradient boundary. The original ESRP model grid was divided into square cells approximately 4 mi per side. After the model area was reduced, the model area was rediscretized with each grid cell divided into approximately 1-mi square cells. Therefore, the number of grid cells increased by 16 times. A schematic drawing of this operation is illustrated in figure 2.
The subregional model includes only the western portion of the ESRP model. Flow values across column 25 of the ESRP model were input as constant flow into the subregional model cells along the upgradient boundary. The flow was input into these cells as line wells. The model was run with this new configuration and resultant water levels and budget values were checked against the original model. Results were nearly identical to the original model.
No changes were made to the model configuration that affected ground-water elevations, discharge from springs, recharge sources, and other critical values or conditions. Also, we did not change the hydraulic parameters from the original ESRP model. The only change to the subregional flow model resulted from the use of a much smaller grid cell size than the one used in the ESRP model. The ESRP model incorporated the use of wells as recharge sources. Because of the smaller grid size used in the subregional model, we needed to increase the number of cells into which recharge was input in order to more accurately incorporate the recharge volume; otherwise, the model would not match the original results. Therefore, the total volume of recharge remained the same as in the ESRP model, but this volume was distributed over 16 times the number of recharge wells. Wherever possible, we placed these additional recharge wells in the rediscretized adjacent grid cell blocks. The volume of each original recharge well was evenly distributed between the 16 smaller cells. A few minor adjustments in the locations of these new recharge wells were needed so the new simulations matched the original simulations, especially along the model boundary.
We simulated the transport of dissolved nitrate in the mid-Snake region’s ground-water system by using the USGS three-dimensional solute-transport model MOC3D (Konikow and others, 1996, Kipp and others, 1998), which was developed for use with MODFLOW 2000. MOC3D simulates three-dimensional solute transport in flowing ground water, computing changes in concentrations of a dissolved constituent over time that can be caused by factors such as advective transport, hydrodynamic dispersion, and mixing or dilution. MOC3D uses the method of characteristics to solve the transport equation for a given time step on the basis of hydraulic gradients computed with MODFLOW (Konikow and others, 1996, Kipp and others, 1998). Particle tracking is used within the model to represent advective transport. The conceptual model of nitrogen input developed by Skinner and Donato (2003) was used in the transport model as a GIS, nitrogen-input data layer. This nitrogen source layer was used to define nitrogen flux at the surface in the flow and transport model. The flow and transport models were simulated using steady-state conditions in three dimensions.
Data collected or compiled during recent and ongoing studies that define ground-water nitrate concentrations were used to refine the new model and to evaluate model calibration. All concentrations of nitrate are presented consistently as milligrams per liter of nitrite plus nitrate as nitrogen with negligible nitrite concentrations in the area.
Data used for nitrogen input into the solute transport model were calculated as a nitrogen-input layer during a separate study, and described in Skinner and Donato (2003). The input was updated to account for differences in study area extent between the two studies. This study’s area includes additional areas of Bingham, Blaine, Butte, Camas, Elmore, and Power Counties; and less area of Cassia, Gooding, Lincoln, and Twin Falls Counties. The areas of the original nitrogen-input layer not in this study’s area were simply “clipped off.” The expanded portions of this study area are primarily upgradient to the solute transport model, and represent background conditions for nitrogen input. We assigned nitrogen-input values to the expanded study area portions that corresponded to similar approaches in the adjacent counties (fig. 3).
Input from the nitrogen source layer was incorporated in the MOC3D transport model by using the area wells function in all active model cells in layer 1. The area wells function allows for the nitrogen source layer load to be introduced into the flow model by assigning a volumetric stress to each cell. The stress applied is a small volume of water that contains a specified concentration of nitrogen derived from the source layer. The initial pre-agriculture background concentration applied throughout the model was 0.5 mg/L in the upper model layer, 0.4 mg/L in the second and third underlying layers, and 0.3 mg/L in the fourth underlying layer. A very small stress term was initially given, and this term was adjusted based on available nitrogen in ground-water data in the USGS National Water Information System (NWIS) database. These data were used to refine the new model and to evaluate its calibration. The calibration points included about 300 wells located throughout the model area, and 11 springs that discharge into the Snake River. The well data usually included one analysis per well. The data for most of the bedrock springs discharging from the aquifer included many analyses over numerous years (Clark and Ott, 1996). The model continued to be simulated as steady state because the additional stress term was small in relation to overall stresses.
Calibration of the transport results was accomplished by comparing areal nitrate concentrations results from the simulations with nitrate data from the wells and springs. No rigorous calibration techniques, such as parameter estimation, were used in the calibration process because such techniques likely would have required changes to the underlying hydraulic parameters in the already-calibrated ESRP model, which was not part of the process. Therefore, calibration was done mostly using trial and error techniques.
Calibration began by taking the nitrogen source layer as representing current source conditions, and then adjusting it to start in 1910. County agricultural statistics for the six counties that overlap with Skinner and Donato (2003) were used to determine the total numbers of agricultural acres and dairy cattle. Those statistics were used to approximate the percentage of current nitrogen load that could be used to represent the 90-year period. The amount of available nitrogen loading was increased during sequential model time steps by a factor related to those statistics (fig. 4). The result of the simulation, therefore, represents the end of a 90-year historical period simulation based on the increase in the area’s crop farming and the introduction and increase of dairy farms.
Comparing maps of measured ground-water nitrate concentration data from the Skinner and Donato (2003) study area and simulated ground-water nitrate concentration data from this study and graphs of measured and simulated spring nitrate concentration data showed that the model was reasonably able to simulate the general pattern of measured nitrate concentrations (fig. 5).
The springs flow from a bedrock wall, in a range of elevation of hundreds of feet and are therefore represented in model layers 1 and 2. We input the vertical locations of these springs into the model, with most springs located in model layer 1. The results of simulations indicated that the model was capable of approximating the nitrate concentrations for springs in both model layers (fig. 6). Simulated nitrate concentration was about 0.5 mg/L higher than data collected for a 3-year period from Box Canyon Spring and showed similar intra-annual variability. Simulated results represent nitrogen loading conditions from the same period as the measured data. The simulated results were acceptable for the objective of this project and model to simulate regional nitrate variability in response to changes in proportions of nitrogen sources.
Simulation results indicated that the nitrogen source layer was not necessarily the driving factor in the resulting nitrate concentration map; rather, the ground-water velocity appeared to be the prime controlling factor (fig. 7). Some areas with large nitrogen loads and large ground-water velocities resulted in relatively small simulated nitrate concentrations, and other areas with relatively small nitrogen loads and small ground-water velocities had large simulated ground-water nitrate concentrations. The simulated results reasonably represent the measured data.
Data gathered from previous studies in the mid-Snake region, including USGS RASA and NAWQA programs, allow for the development and assessment of detailed conceptual and computer simulation models. These data and regional-scale RASA ground-water flow models are available for much of the Snake River Plain. Therefore, the approach used in this study is one potential assessment option for nitrate contamination issues in other areas of southern Idaho. However, the accuracy of such models will be limited because of the simplifying assumptions used in the models and uncertainty in model parameters. These limitations must be understood and considered when applying the model to resource-management issues.
U.S.
Department of the Interior | U.S. Geological
Survey
Persistent URL: https://pubs.water.usgs.gov/sir20065137
Page Contact Information: Publications Team
Page Last Modified: Thursday, 01-Dec-2016 19:16:13 EST
