Groundwater-Flow Simulation


The primary goal of the groundwater-flow simulation analysis was to evaluate the relative contributions of pumping and commingling to the persistent post-1970 declines in groundwater levels in the OWRD administrative area. Following development of a numerical model of groundwater flow using MODFLOW-2000 (Harbaugh and others, 2000) the aquifer system analysis was conducted in three sequential phases using the model-independent parameter estimation software, PEST (Doherty, 2005; Doherty, 2010). At the completion of this analysis, a set of “best estimated” aquifer parameters were selected for presentation in this report, and limitations of the resulting model were identified. The three phases and the principal conclusions of each are:


1. Rough calibration of a pre-development condition steady‑state model – Initial values of hydraulic parameters were developed from literature values, followed by calibration to the earliest available groundwater level measurements in each aquifer and region of the groundwater model area. Even though there are only two pre-1950 (pre-development) wells with data, it was reasonable to use other early data because documented groundwater declines are relatively small when compared to the range of water-level elevations across the model area (75–1,750 ft). Use of this data allowed a rough calibration of the model, providing reasonable estimates of parameters for use in subsequent analyses. The calibrated model fit the data reasonably well, implying that the conceptual model of flow and the representation of hydrogeology in the groundwater-flow simulation model are reasonable. 

2. Transient modeling of the groundwater-flow system – A time-varying (transient) version of the MODFLOW model was developed to simulate system hydraulic response during the period of development of water resources. This version of the model performed poorly due to limitations of MODFLOW-2000 when representing the large effect of commingling wells and due to the complex storage characteristics of the aquifer system. To conduct a fully transient analysis, a different groundwater-flow simulation code may need to be used.

3. Modified transient analysis – To separate the complex storage problem from the analysis of the relative effects of pumping and commingling, the transient problem was divided into four steady-state simulations to represent the four steady-state conditions shown in figure 2. This formulation of the problem removed the uncertainty associated with estimating the storage terms. Model results indicated that greater than 80 percent of the observed aquifer declines in the OWRD administrative area are attributable to commingling, with the remainder being the result of pumping. A subsequent analysis of water storage mechanisms shows that the long-term linear declines may be the result of draining higher hydraulic head aquifers to supply the confined aquifer system in the administrative area.


The following sections summarize the groundwater-flow simulation model construction (methods used to simulate hydrologic boundaries and assign hydraulic properties to hydrogeologic units) , the modified transient analysis and results, and subsequent groundwater-flow simulation analyses designed to aid in selection of management actions intended to restore groundwater levels. Additional details of the simulation analyses are provided in appendix E. 


Model Discretization and Boundaries


The groundwater-flow simulation model was created by dividing the model area into 500-ft on a side square model cells of variable thickness and representing springs, streams, rivers, and wells (fig. 22). The thickness of each flow model cell was defined using the hydrogeologic framework, with each layer corresponding to a single hydrogeologic unit in most cases. Model boundary conditions are shown for each layer in figures A9 through A22. The relation between groundwater-flow model layers and the geologic and hydrogeologic layers is shown in figure 3.


Lateral Boundaries


The extent of the model was based on natural hydraulic boundaries. The model area is bounded to the north by the Columbia River, to the south-east and west by ridgelines with adjacent faults, and to the south by a segment of the Chenoweth thrust fault (fig. 2). Water is allowed to flow to the Columbia River, but all other boundaries are considered to be impediments to flow. 


Where aquifers are in connection with the Columbia River, a general head boundary allows water to flow to or from the river at a rate proportional to the difference in hydraulic head between the aquifer and the river (fig. 22). Where aquifers are exposed above the river along its shore, a drain boundary allows water to drain out of aquifers with the controlling drainage elevation being the aquifer bottom elevation at the modeled outcrop. This assumption allows simulation of springs along the shore of the river.


The axis of the Columbia Hills anticline forms the ridgetop to the southeast. Immediately to the southeast of this ridgetop is the Chenoweth thrust fault, which is draped down the slope into the next valley. This fold-fault combination is likely to form an effective flow barrier and is simulated using a no-flow boundary condition. Similarly, the western boundary is a ridgetop immediately to the east of a wrench fault of significant offset, also forming a flow barrier and being simulated using a no-flow boundary condition. 


The southern boundary corresponds to a mapped extension of the Chenoweth thrust fault system. Two factors support using this feature as a flow model boundary for the groundwater-flow model:


1. The geologic modeling results indicate that many of the upper aquifer units pinch out to the north of the Chenoweth thrust fault (figs. 2, 4, and 5). As a result, most aquifers of interest are not continuous to this point. 

2. The Rocky Prairie thrust fault is known to be a major impediment to groundwater flow, so by analogy, the Chenoweth thrust fault is hypothesized to be a flow impediment. Analysis of streamflow data and watershed modeling indicate that most of the groundwater is forced into the stream system above this point, providing evidence that the thrust fault is an effective flow barrier (see Discharge to Surface Water). 


Initially, the Chenoweth thrust fault was modeled as a no-flow boundary, with the expectation that later modeling might use a general head boundary to test the effects of flux across this boundary. Results from the modified transient analysis precluded the need for adding additional water at this boundary, because additional water would require that commingling fluxes be even higher. Because commingling was already shown to be the principal cause of declines, the no-flow condition was the conservative case for estimating commingling effects. 


Faults


The importance of faults as a control on groundwater flow was evaluated using the groundwater-flow model. In the groundwater-flow model area and below the overburden, the faults with significant offset (potentially juxtaposing aquifers and confining units) which were identified during geologic modeling were simulated as possible impediments to lateral flow through the basalt aquifers (figs. A9 through A22) using the Horizontal Flow Barrier (HFB) package (Hsieh and Freckleton, 1993). Values of conductance of the faults were allowed to vary during calibration and uncertainty analysis. Faults were divided into segments laterally and vertically to test the possibility that water may move more easily through parts of the fault based on fault geometry. For example, the Maupin wrench fault (fig. 2) exhibits high offset to the southeast, and little or no offset near the Rocky Prairie thrust fault, indicating that the hydraulic conductance will vary along the fault. Similarly, the Rocky Prairie thrust fault can have different properties with depth, especially because the overthrust part of the fault likely only contains Priest Rapids Basalt and younger strata. Fault segment hydraulic conductance values were regularized by indicating that adjacent (vertically or horizontally) segments within the same fault (for example, the Rocky Prairie thrust fault) likely have similar values of hydraulic characteristic (the MODFLOW parameter defining fault conductance).


Streams


In the study area, the aquifers contribute groundwater to streams more efficiently than they gain water from streams. Study area streams lose water to the groundwater system only when the relatively thin aquifers are in direct connection with the stream and stream stage is above the aquifer hydraulic head. If aquifer head is above the stream stage at these locations, then water will flow the other direction, but water may also drain out of the aquifer system into streams through springs and seeps that are above the stream stage. A combination of the MODFLOW drain and stream packages (McDonald and Harbaugh, 1984; Prudic, 1989, respectively) was used to simulate study area creeks. Groundwater-flow simulation model cells containing perennial streams were simulated using both drains and streams, and only drains were used in cells containing ephemeral streams (fig. 22). Details are provided in appendix E.2.


Only two perennial streams originate outside the groundwater model area. These are Rock and Mosier Creeks, originating to the south of the model boundary. Because the reach of Mosier Creek immediately to the north of the Chenoweth thrust fault loses water seasonally, estimated streamflow at the groundwater model boundary was required. This flow was estimated as the total average annual recharge to the drainage area contributing to the stream above the point where it enters the groundwater model area. Inflow at the groundwater-flow simulation model boundary was estimated to be 16.5 ft³/s for Mosier Creek and 0.911 ft³/s for Rock Creek. The Mosier Creek estimate is approximately 80 percent of the average annual base flow estimated at the Mosier Creek gaging station (fig. 1), agreeing well with seepage run observations (see section Discharge to Surface Water).


In addition to using drains in model cells intersecting streams, drains were used to represent springs and seeps that would form where aquifers are exposed at land surface if groundwater levels were sufficiently high. Flow from these cells is assumed to flow into streams further downslope. Total groundwater discharge to streams along specified reaches were computed by summing all drain and stream cell fluxes in the contributing drainage area using ZONEBUDGET (Harbaugh, 1990). 


Recharge


Recharge estimates were derived from the PRMS watershed modeling results (see Recharge in the Groundwater‑Flow System section, appendix B, figure B4). Monthly recharge was estimated for use as input during transient modeling, but the average annual recharge rate was used for the modified transient analysis and all subsequent management scenario simulations. PRMS recharge estimates were made for the entire groundwater model area except for a small relatively low recharge area along the Columbia River (fig. B4). Recharge in this area was estimated as a simple average of recharge for adjacent PRMS model units. 


Wells


Several model-area wells are known to commingle the basalt aquifers. Of principal concern to this study is the extent to which vertical intraborehole flow contributes to groundwater-level declines, so representation of the about 500 wells in the model was a critical piece of model formulation. Of these wells, 25 wells account for about 80 percent of the total pumping (see Pumping of Groundwater section). These wells also are typically some of the deepest wells in the area, so accurate representation of the wells is important.


Historical well construction practices in the area have resulted in open boreholes between multiple aquifer zones. Even wells with a casing or liner installed commonly have no functional seal except near the ground surface, allowing hydraulic connections between aquifers through the annular space between the borehole and the casing or liner. The geometry of an important subset of wells was examined in detail, and aquifers pumped and commingled by these wells were assigned to groundwater-flow model units using the geologic model (diagram of this relation shown in figure 3), well logs, and best professional judgment. These wells include all of the high capacity wells and wells installed before and a few years after the declines of the 1970s began. Understanding the geometry of these wells is critical for understanding system response. Groundwater-flow model layers for the remainder of the wells were assigned based on the method described in the Commingling Wells section, providing a reasonable distribution of pumping and commingling for the flow model. 


The MODFLOW Multi-Node Well package (Halford and Hanson, 2002) was initially used to simulate commingling wells, but numerical instability of the model resulted. To correct this problem, intraborehole flow between aquifers was simulated by using high vertical conductivity of cells containing commingling wells, and distributing pumping between aquifers based on the transmissivity of the aquifers (see appendix E.2 for details). If the aquifers were not commingled to the extent assumed in the parameterization described previously, then during the calibration process, vertical conductance of well cells was decreased until the net effect of commingling was represented. If commingling was occurring in many wells, the net effect was an increased vertical conductance. In this way, the calibration process revealed the role of commingling wells in controlling study area flow patterns. 


Flow and Storage Properties of Hydrogeologic Units


Each hydrogeologic unit was assigned a hydraulic conductivity value, which controls the ease with which water flows through the unit, a specific storage value, which represents how much water is stored under confined-flow conditions due to the compression of water and rock, and a specific yield value, which represents the water that would be released if it were drained under unconfined-flow conditions. These values were assumed to be constant across the study area for each hydrogeologic unit. Additionally, the hydraulic properties are assumed to be similar between units of similar geologic character. For example, CRBG flow interiors are assumed to transmit and store water similarly to each other. Details are provided in appendixes E.3 and E.4. 


Groundwater-Flow Model Analyses


The groundwater-flow simulation model was used to evaluate the importance of hydrogeologic controls on groundwater flow and storage. Adjustable parameters were selected for evaluation during the modified transient analysis, and a set of values for these parameters were selected for use in subsequent analyses. The adjustable parameters, selected based on the conceptual model of the system, were horizontal and vertical hydraulic conductivity of all flow model units, hydraulic conductance of the Rocky Prairie thrust fault and the Maupin wrench fault, conductance of streambeds and drains, and precipitation-derived recharge rate.


Parameterization is the process of dividing, grouping, or fixing adjustable parameters based on likely similarities or differences in hydraulic behavior. Tikhonov regularization (Doherty, 2005) was used to represent likely relations between independently adjustable parameters. Parameterization and regularization for major parameter groups include:


1. Horizontal and vertical hydraulic conductivity of groundwater-flow model units – Every flow model unit was initially assumed to have isotropic hydraulic conductivity, because CRBG aquifer system horizontal to vertical anisotropy is assumed to result from the contrast between the aquifers and confining units, and these units were separated and represented explicitly in the Mosier area groundwater-flow model. However, the Upper Undifferentiated Overburden confining unit was represented as anisotropic for the final management scenarios analyses to prevent groundwater levels from being simulated above land surface in this unit in the western part of the study area. This was done to prevent bias in estimates of changes in groundwater storage. Conductivity of each model unit was independently adjustable, but regularization was used to create groups with similar values. These groups are the basalt aquifers, basalt confining units, and interbeds. Vertical conductance of commingling well cells was an independent parameter.

2. Hydraulic conductance of faults – The Rocky Prairie thrust fault and the Maupin wrench fault were divided into 8 segments in map view and 6 segments vertically (one for each simulated aquifer-confining unit combination, fig. 3), resulting in 48 independently adjustable parameters. Regularization was used to indicate that adjacent fault segments likely have similar properties, while allowing fault conductance to vary with depth or along the fault trace.

3. Streambed and drain conductance – Adjustable stream conductance and drain conductance parameters were defined for each groundwater-flow model unit. Regularization was used to create groups with likely similar behavior (for example, drains and streams intersecting basalt aquifers are likely more similar to each other than drains and streams intersecting basalt confining units). During calibration, the conductance of drains located in several low permeability units was tied to the conductance of drains located in the adjacent aquifers, preventing these insensitive parameters from becoming arbitrarily large. This is consistent with the assumption that more water will drain through the aquifers into streams than will drain directly from low-permeability confining units into the stream (though water can drain from confining units into aquifers). Because altering streambed and drain conductance have similar effects on model results, regularization was used to minimize streambed conductance. This regularization condition allows streambed conductance to be sufficiently large to account for stream loss to the aquifer system where data support the need for stream leakage, while allowing drains to account preferentially for most of the base flow contributions.

4. Precipitation-driven recharge – Recharge was divided into six hydrogeologic zones (based on surficial geology), each with a parameter to adjust the fraction of PRMS-derived recharge to use as groundwater-flow model input. 


Analysis of Persistent Decline of Groundwater Levels


Following poor performance of a fully-transient simulation model, three steady state models were used to simulate the final steady conditions that would result under continued pre-development, early-time (circa 1970), and late-time (circa 2006) stress conditions (fig. 20). This analysis allows examination of the magnitude of declines observed in Group 1 wells, independently of the aquifer-system storage parameters, and is referred to hereafter as the modified transient analysis. Model parameters were identical for all three models, with the only difference being configuration of commingling wells and pumping rates (table 2). The three configurations were: 


1. Pre-development (prior to 1950)–No wells.

2. Early time (the few years prior to onset of persistent declines)–Pumping rates and commingling well geometries present during 1970.

3. Late time (in the future)–Pumping rates and commingling well geometries during 2006.


Calibration


Calibration is the process of finding reasonable values of model parameters that result in the best model fit to the measurements. Model fit is a measure of how well the simulated values match the measured values (observations). For a complex model with a relatively large amount of data, fit will seldom be perfect, but it is often good enough to allow the model to be useful to gain understanding of system behavior. 


The three steady state models were calibrated to groundwater level and streamflow data representative of each period. Because the groundwater system was apparently approaching dynamic steady-state prior to 1970, groundwater‑level calibration targets for the pre-development and early-time periods were estimated as the median value of winter groundwater-level measurements during the corresponding period, resulting in 2 pre-development and 12 early-time groundwater-level calibration targets. Winter levels have been collected from many wells historically by the OWRD and are assumed to represent the long-term effects of pumping rather than the seasonal drawdowns associated with summer irrigation. Late-time steady-state groundwater-level data do not exist, although groundwater levels must be at or less than their current levels. Current winter groundwater levels were selected to provide 46 calibration targets, with the understanding that modeled values should be equal to or lower than these estimated values. Calibration target weights were lowered to account for the larger number of late-time targets and the lower confidence that the target values represent true steady-state conditions. Groundwater levels were assumed to represent the hydraulic head at the bottom of each well, which corresponds to the deepest aquifer intersected. For a large rate of flow between commingled aquifers, effective vertical conductance of well cells is large, and the hydraulic head difference between the commingled aquifers is correspondingly small, which makes the assumption that groundwater levels represent the deepest aquifer a reasonable assumption. The assumption that head differential between aquifers at commingling wells is small was examined at each stage of parameter estimation and prediction, and proved to be reasonable for all sets of parameters that result in a calibrated groundwater-flow model. 


Spatially distributed average annual PRMS-derived recharge (see section Recharge for details) was used for all steady-state simulations; 100 percent of PRMS-derived recharge was initially used, resulting in adequate model performance, so these parameters were not adjusted for most of the analyses. Following determination that commingling was the dominant cause of declines, recharge was reduced to 90 percent to evaluate the uncertainty associated with possible errors in estimation of recharge. The commingling effect was smaller but still dominant, indicating study results are valid unless PRMS derived recharge estimates are much too high. 


The average annual PART-derived base flow estimate upstream of the Mosier Creek gaging station for the period 1964–81 (see section Discharge to Surface Water for additional detail) was used as a streamflow target for all three steady-state simulations. Because groundwater levels declined during this period, the simulated base flow will be higher than the PART-derived estimate for the predevelopment time period. Similarly, the late-time period simulated base flow will be lower than the PART-derived estimate. PART-derived base flow values were used as calibration targets instead of seepage data because seepage data are available only during the dry season and do not represent the average annual values being simulated (fig. B2). However, simulation results were consistent with spatial patterns of measured seepage.


The modified transient model was calibrated (fig. 23) using PEST (Doherty, 2005), providing reasonable best estimates for model parameters (“Best Estimates for Modified Transient Analysis” in fig. 24). Observations were divided into four groups for use with PEST: pre-development groundwater levels, early-time groundwater levels, late-time groundwater levels, and all stream base flows. Initially, observation weights for groundwater levels for each simulation period were assigned inversely proportional to the number of observations in the group so that periods (groups) with fewer observations had higher weights per observation, providing similar importance to calibration data from each period. This strategy ensured that each groundwater-level observation group (each representing a time-period) provided a non-negligible contribution to the calibration objective function. The stream base-flow group weights were reduced relative to groundwater-level observations to account for measurement unit differences, while ensuring the contribution of this group to the calibration objective function was also non-negligible. This strategy ensured that each observation group provided constraints on model calibration and predictive uncertainty analysis. Pre-development, early-time, and late-time base flow were placed in a single observation group because the PRMS estimate of average annual base flow (1955–2007) was used for all three periods with the expectation that the pre-development and early-time groundwater flow simulation model estimates will be larger than the PRMS value and the late-time estimate will be smaller than the PRMS value.


Adjustment of weights assigned to individual observations within groups was accomplished on a case-by-case basis, with lower weights assigned to observations that were not representative of the processes being examined. Of the 63 measurements and estimates used for calibration, model fit for 2 early-time and about 12 late-time groundwater levels were persistently too low, otherwise fit was good (fig. 23). The reasons for poor fit at these points were examined, and the bias was deemed acceptable for the modeling purposes of analyzing Group 1 groundwater-level declines. The poorly fit groundwater levels were de-emphasized in the automated calibration method by using low observation weights (fig. 23B). Three of the poorly fit points (1 early time and 2 late time) represent water levels in the undifferentiated overburden. Because this unit is thick and highly heterogeneous, the measured groundwater levels are assumed to represent perched groundwater that is recharging the deeper aquifer system. The remaining 11 poorly fit points represent the mid-slope area between the Rocky Prairie thrust fault and the ridgetop (fig. 5A–A′). Calibration with equal observation weights, improved model fit in this area by reducing hydraulic conductivities of the aquifers, but caused excessively large horizontal hydraulic gradients in the OWRD administrative area. This result suggests that there is a conductivity contrast between the upslope and downslope portions of the system. This contrast may be the result of a flow barrier created by a fault or fold, or there may be gradational lateral changes in aquifer hydraulic conductivity created by varying depositional characteristics of the intra-canyon basalt flows. The flow margin of all three upper aquifers is mid-slope (fig. 5A–A′).


Because the location and nature of the mid-slope conductivity contrast is not known, representation in the flow model would be uncertain. Rather than developing an uncertain model to represent the conductivity change, the observation weights were reduced for the poorly fit data. De-emphasizing these data is reasonable because (1) the modeling objective is to represent the system behavior in the OWRD administrative area, and (2) the amount of water flowing into the OWRD administrative area is controlled by the prescribed recharge. Streamflow measurements collected for seepage analysis (fig. 11) showed no obvious barriers to flow forcing large amounts of water into Mosier Creek in the area of the likely conductivity contrast (mid-slope), providing evidence that groundwater flow into the OWRD administrative area does not strongly depend on the conductivity contrast. 


The resulting calibrated model is appropriate for testing the effects of changes in recharge, pumping, and commingling in the OWRD administrative area. Although the model generally represents aquifer-system behavior, groundwater levels predicted outside of the administrative area are less reliable. The computed groundwater budget is reasonable (table 2) with the model simulating that early time commingling of the upper aquifers resulted in increased streamflow to local streams and reduced flow to the Columbia River. Conversely, for late time, commingling is predicted to increase flow to the Columbia River at the expense of local stream base flow contributions. These results are reasonable given system geometry, but because a predictive uncertainty analysis was not performed, these results are somewhat uncertain, and it is reasonable to assume that combinations of parameters might exist that show reductions in Columbia River base flow for early and late times with corresponding base flow increases to local streams. 


Evaluation of Model Parameters


Calibrated parameters match the conceptual model and previously collected data. In addition to the sets of parameters selected as “best” for the previous analysis and for management scenarios (fig. 24), the range of possible parameter combinations was explored during a series of calibrations where starting parameter values and calibration strategy were varied to assess the sensitivity of the estimated parameter values. The analysis was not exhaustive, but did provide confidence that the general values of parameters and relations between parameter values were reasonable. The most sensitive parameters varied between automated calibration runs, with some parameters being highly sensitive at the end of one run but being relatively insensitive at the end of another run.


The CRBG aquifers are the most permeable units, with the younger units generally more permeable than the older units. Lite and Grondin (1988) conducted two one‑day pumping tests of the upper aquifers in the OWRD administrative area, estimating transmissivity that is equivalent to hydraulic conductivity values on the order of a few thousand feet per day. These data were used as prior information in the model calibration with low weight, and the calibrated values of hydraulic conductivity for the upper CRBG aquifers were somewhat lower. The lower values are consistent with the observation that upslope aquifer hydraulic conductivity may be lower. Calibrated values of hydraulic conductivity for CRBG flow interiors were 5–6 orders of magnitude less than their associated flow top aquifer values (fig. 24A). 


Modeling results indicated that flow barriers associated with faults are important for controlling groundwater flow (fig. 24A). The Rocky Prairie thrust fault consistently had the lowest conductance, with permeability possibly increasing to the east. The Maupin wrench fault was also a barrier with a general trend of higher conductance to the north. Calibrated values of fault conductance indicated that some faults might be less permeable at depth, although trends were not strong.


All of the aquifer and confining unit hydraulic conductivities and fault conductances were sensitive because of the regularization constraints defining the likely relations between units of similar or contrasting properties (fig. 24B). However, the effective conductivity of commingling wells was not always sensitive. When the well conductivity became too large, it became insensitive (compare fig. 24A and B) because, at high values, aquifer hydraulic conductivity limited flow to the commingling wells. The drain and stream conductances for aquifers and overburden had a range of sensitivities (fig. 25), most falling within the range shown in figure 24B. Stream and drain conductances for the low permeability confining units were frequently insensitive, so these were tied to adjacent aquifer stream and drain conductances. Because both stream and drain cells can control the rate at which groundwater leaves the aquifer system, they may act as surrogates for each other during the calibration process. To prevent surrogacy, regularization was used to find the minimum values of stream conductance that result in a calibrated model (appendix E). This allows drain conductance to control the rate of leakage from aquifers to streams, except in reaches where data suggest that aquifers are gaining water from streams. 


Separation of Pumping and Commingling Effects


The effects of pumping and commingling were separated by simulating the cessation of all pumping in the model area, but leaving the commingling wells in place. This was accomplished by adding a fourth steady-state simulation to the modified transient analysis model and estimating the groundwater level recovery in all of the Group 1 wells (conceptually shown in fig. 21). Using the calibrated best estimates from the modified transient analysis, the flow model predicted very poor recovery of groundwater levels in the 24 Group 1 observation wells (fig. 26), indicating that commingling may be a large contributor to groundwater-level declines.


PEST was then used in predictive mode (Doherty, 2005) to find the set of hydraulic parameters that still match the data (to a 95-percent confidence, details in appendix E), but that provide the maximum recovery of Group 1 groundwater levels. During this analysis, PEST allows model fit to degrade in order to find the set of parameters that predict the maximum recovery of water levels. The best recovery predicted by the model was typically less than 30 ft for Group 1 wells (fig. 26), with the maximum predicted recovery still accounting for less than one-half the current declines. Some of the parameters attained unlikely values during the analysis, indicating that the recovery predicted by the best-calibrated model (“best estimated parameters” on fig. 26) may be a better estimator. In conclusion, commingling is likely the dominant cause of declines, accounting for greater than 85 percent of the observed declines in Group 1 wells.


The predictive uncertainty analysis above did not consider changes in the zero-flux condition at the Chenoweth thrust fault boundary. However, addition of more water into the aquifer system across this boundary would require that the Maupin wrench fault be less permeable or that the commingling fluxes be even higher to remove additional water. For this reason, the zero-flux condition was a conservative assumption, preventing the need to simulate possible flow at this boundary.


Limitations of the Groundwater-Flow Simulation Model


In addition to the limitations and assumptions described above, the resulting groundwater-flow simulation model and calibration data were examined to identify limitations that may affect future application of the model. This was accomplished by exploring a range of starting parameter values and model calibration strategies. These limitations may be framed in the context of the data required to improve the model for future uses. Six groups of data were identified:


1. There is insufficient data to constrain groundwater-flow conditions to the west of the Maupin wrench fault (fig. 2). In particular, groundwater levels in the thick sedimentary overburden to the southwest, far from the area of current interest, are simulated as being above land surface. This is the result of having only a few groundwater‑level measurements for the low-permeability upper undifferentiated confining unit in the eastern part of the study area where data constrain the vertical hydraulic conductivity, which controls the rate at which water is transmitted through the overburden to the aquifers. However, this unit is heterogeneous and stratified, and the effective horizontal hydraulic conductivity controls the rate at which water is transmitted laterally to streams in the western part of the study area (for example, West Fork Mosier Creek (fig. 4). Increasing the horizontal hydraulic conductivity of this unit would allow groundwater levels in this unit to be at or below land surface at all locations.

2. There are insufficient groundwater-level measurements to constrain groundwater-flow conditions through the Rocky Prairie thrust fault and to the Columbia River for the Frenchman Springs and Grande Ronde aquifers. Although the Rocky Prairie thrust fault is a flow barrier to the stratigraphically higher units, the thrust detachment is likely above the Frenchman Springs and Grande Ronde units. However, there may still be deformation of the deep aquifers due to compression that results in restriction of lateral flow. Similarly, there is no data to constrain the connection between these deep units and the Columbia River. As a result, elevated groundwater levels in these deeper aquifers to the southeast of the thrust fault may be explained by poor hydraulic connection across the thrust fault or with the Columbia River.

3. There are insufficient groundwater-level measurements to constrain groundwater-flow conditions in the deep Grande Ronde aquifer system, with only a few groundwater level measurements made in this unit near the crest of the Columbia Hills anticline (fig. 2). Although there is documented and anecdotal evidence that hydraulic‑head gradients indicated upward flow in the OWRD administrative area, this evidence is for the upper Frenchman Springs aquifers and above (younger strata), and does not extend to the older, deeper Grande Ronde aquifers. The absence of groundwater-level measurements for the Grande Ronde aquifer in the administrative area, coupled with the absence of Grand Ronde unit exposure in aquifer recharge areas upslope (fig. 2), indicate that not only are groundwater-flow conditions in this unit poorly understood, but that the Grand Ronde may not be a reliable long-term source of groundwater in the study area. Additionally, it is possible that the hydraulic gradients and flow are downward into the Grande Ronde in the OWRD administrative area.

4. There are insufficient groundwater-level measurements to establish the pre-development vertical hydraulic head differences between the aquifers in the OWRD administrative area. As a result, there are nonunique sets of model parameters that can be used to calibrate the groundwater-flow simulation model where the parallel aquifers act as surrogates for each other, transmitting water from the upper to the lower parts of the study area. This generally is accomplished with relatively modest changes in hydraulic conductivity of the confining units. This is because most of the calibration data does not represent the extent to which the aquifers were hydraulically separated before installation of commingling wells. Most groundwater-level measurements represent late-time conditions where the restriction of vertical flow between adjacent aquifers due to the confining units has been greatly diminished by commingling wells.

5. There are insufficient detailed data to determine which wells simulated as commingling are commingling in reality. For this reason, the groundwater-flow simulation model is not an effective tool for evaluating the optimal sequence of well repair.

6. There were too few stream gaging stations with continuous records to allow the use of streamflow data to test hypotheses about changes in spatial distribution of base flow resulting from pumping and commingling. The model predicts a late-time reduction in average annual base flow at the Mosier Creek gaging station (fig. 1) that falls within the noise of the natural variability of the data. Additional streamflow observations at key locations along Mosier Creek likely would have allowed an analysis of changes in base flow upstream and downstream of the current gaging station between the Rocky Prairie and Chenoweth thrust faults. Additionally, continuous streamflow measurements above and below the thrust faults likely would have improved the understanding of the role of these faults as barriers to groundwater flow.


Not all of the previous data limitations are equal for any given future modeling purpose. Of these identified weaknesses, the principal improvement in knowledge necessary to aid in restoration of groundwater levels is the distribution of wells that are commingling in reality.


Evaluation of Potential Management Options


Because commingling likely is the dominant cause of groundwater declines, the groundwater-flow simulation model was used to identify areas most vulnerable to commingling wells and to evaluate combinations of commingling well repair and artificial aquifer recharge that might be used to restore groundwater levels. In all cases, the flow model was used to calculate total change in CRBG aquifer system storage, and the percentage change in CRBG aquifer storage was used as the metric of comparison. Because this performance measure was selected, assumptions were made to address some of the limitations of the groundwater-flow simulation model (see section Limitations of the Groundwater‑Flow Simulation Model). These assumptions resulted in two additional constraints and one additional adjustable parameter, and the model was recalibrated to estimate a new set of “best” model parameters for use in evaluating management scenarios. The constraints reflect anecdotal information, best professional judgment, and information gained during the modified transient model analysis. Unlike for the modified transient analysis of commingling and pumping, a predictive uncertainty analysis was not performed. As a result, estimates of change in storage should be used only for comparison between scenarios and not to make absolute estimates of expected groundwater-level recovery at any single point.


The first additional constraint is that pre-development water levels in the overburden must be at or below land surface. This constraint was added because groundwater levels above land surface were simulated near the western boundary of the study area frequently during the calibration process, and these simulated levels are likely erroneous. High simulated groundwater levels in this area probably result from the modeling assumption that the upper undifferentiated overburden is isotropic (despite its high heterogeneity) and from the lack of calibration data to the west to constrain simulated groundwater levels in the overburden to below land surface. The additional constraint was imposed to prevent erroneous groundwater levels in this area from unduly biasing computation of changes in total basalt aquifer storage. 


The previous calibrated values of hydraulic conductivity for the undifferentiated overburden unit were controlled by the vertical permeability of this unit near the OWRD administrative area where the unit is less extensive but few calibration data exist. In this area, the primary function of the overburden is to transmit recharge into the deeper aquifer system. To the west, in the trough created by the Mosier syncline (fig. 2), the Mt. Hood volcaniclastic deposits are likely stratified with alternating sequences of more and less permeable materials. In this area, the permeable layers are likely important for controlling lateral flow to several creeks (fig. 2). To allow the groundwater-flow model calibration process to concurrently simulate lower groundwater levels in the undifferentiated overburden to the west and match calibration data near the OWRD management area, the upper undifferentiated overburden-confining unit was simulated as anisotropic, adding one additional adjustable parameter to the calibration process.


The second additional model calibration constraint assumes an almost uniform upward gradient in the basalt aquifers in the OWRD administrative area. Because commingling is significant, the Lite and Grondin (1988) upward hydraulic head gradient of about 70 ft across the Selah Interbed is assumed to be a lower bound for the pre‑development hydraulic head differential. Anecdotal evidence suggests a significant upward gradient was encountered each time wells were drilled into deeper aquifers in the OWRD administrative area, so a pre-development head difference across each confining unit was assumed to be 100 ft. This estimate is uncertain because virtually all measurements made below the Pomona basalt aquifers have been affected by commingling of the aquifers. Previous calibrations that disregarded the anecdotal evidence resulted in some calibrated models where hydraulic heads between aquifers were smaller than anecdotal evidence supports. This effect likely is an artifact of the sparse pre-development and early time data, and is not a reflection of the model’s ability to replicate this behavior. Inclusion of this second constraint (artificially) rectifies this problem. The uniform upward gradient constraint was enforced by adding a vertical string of artificial observations, one in each basalt aquifer above and including the Frenchman Springs aquifer, immediately to the southeast of the Rocky Prairie thrust fault in the OWRD administrative area. The 100 ft head difference across each confining unit was weighted equally large for the upper three units and 4 orders of magnitude smaller for the confining unit between the Rosalia and Frenchman Springs aquifers. No head difference condition was applied between the Frenchman Springs and Grande Ronde aquifers.


Even though these additional constraints and parameters may improve the groundwater-flow simulation model for many purposes, the previous model (without these additional constraints) is more conservative for evaluating the roles of commingling and pumping. This is true because the extra parameter (a possible new degree of freedom for the predictive analysis) and head observation constraining the upper undifferentiated overburden are largely independent of the declines in the OWRD administrative area (indicating neither the degree of freedom nor the observations constrain the predictive analysis), but the uniform upward gradient requirement (extra observations that constrain the predictive analysis) narrows the range of acceptable models that match the system response (see section Establishing Confidence Intervals for Predictions in appendix E.4). Because the previous model was conservative, the analysis separating commingling from pumping effects need not be repeated for the model with the new constraints.


When the additional constraints and parameter were added to the groundwater-flow simulation model, the model calibrated easily with comparable error to the previous best calibration (compare fig. 27 to fig. 23A). The resulting parameter values (fig. 24) were used for all management scenario analyses. 


Commingling Well Vulnerability Maps


To understand where the aquifer system is vulnerable to commingling, the groundwater simulation model was used to compute change in basalt aquifer system storage resulting from leakage through commingling wells. Vulnerability maps were created by assuming that a larger loss in storage resulting from commingling equates to a more vulnerable area. Because final groundwater-levels in the aquifers are controlled by the elevation of the point where the water leaks out of the system and both aquifers slope upward toward the uplands, placement of a commingling well in a topographically low part of the system would likely result in lower final heads in many aquifers than placement of a similar well farther upslope (fig. 15B). As a result, the low land surface elevation areas of the aquifer system may be considered more vulnerable.


Vulnerability maps were constructed by computing the total change in storage in the basalt aquifers from pre‑development conditions (no wells) to final steady-state conditions resulting from placement of a nonpumping commingling well in a single row and column of the model where more than one aquifer exists. The Grande Ronde flow top aquifer was excluded from this computation because few wells pump from it, and its hydraulic head is poorly defined over most of the study area. The rank, from largest to smallest, of the value of the change in storage was plotted at the center of each model cell (fig. 28). This is a single commingling well vulnerability map, and synergistic effects of commingling at multiple locations are not examined. The pattern of vulnerability shows that the area most vulnerable to commingling is generally coincident with the OWRD management area, which contains most of the Group 1 wells. The vulnerability map also indicates that the area to the west of the Maupin wrench fault may be vulnerable to commingling wells, although this area is currently largely undeveloped.


Simulation of Well Repair Options and Artificial Recharge/Aquifer Storage and Recovery


To consider the potential value of repairing commingling wells, priority well repair zones were identified, and all wells in a particular zone were simulated as repaired. The modeling assumption that, unless a detailed interpretation of well construction information provides evidence to the contrary, each well passing through multiple aquifers commingles all aquifers intersected by the well, works well when predicting aquifer response resulting from the net effect of commingling in many wells. However, there is no guarantee that any single well is actually commingling, because actual commingling depends on the local geology and the well construction method. Without knowledge about which wells are actually commingling, repair zones were defined (fig. 29) using the single well vulnerability map (fig. 28) and all wells in a defined region were simulated as repaired. This approach assumes that the exact location of commingling wells can be unknown, but the net effect of repairing any commingling wells in the zone still can be evaluated.


Two possible methods of recharging aquifers were considered: artificial recharge (AR) and aquifer storage and recovery (ASR) (State of Oregon, 1996). Artificial recharge is accomplished by applying water at the land surface and allowing it to infiltrate into aquifers. Aquifer storage and recovery is the process of injecting water through an injection well, and subsequently removing the water for use. Aquifer recharge from AR and ASR was simulated using the MODFLOW well package (McDonald and Harbaugh, 1984). Because the simulation method is the same for AR and ASR, the term AR/ASR is used in subsequent discussions of simulation results. 


For evaluation of combinations of commingling well repair and aquifer recharge scenarios, the change in groundwater storage in the basalt aquifers was computed (assuming a simple confined storage coefficient), and expressed in terms of percent recovery (table 3). For example, if all groundwater levels recover to pre-development conditions, then it was assumed that the system was restored to 100 percent. The late time simulated conditions representing current pumping and commingling were assumed to be “restored” to 0 percent (not restored). The computed maximum amount of recovery due to repair of all commingling wells was 85 percent, where the remaining 15 percent is attributed to continued pumping. An uncertainty analysis was not done on the computed recovery, so the actual recovery due to repairs may vary from values summarized in table 3. However, the percent recovery may be used to evaluate the relative value of well repair and AR/ASR. The “Repairs only” column in table 3 summarizes the percent recovery of storage resulting from repairing the wells in high priority zones (fig. 29). The number of possibly commingling wells that is simulated as repaired for each scenario is listed in “Number of wells repaired” column (table 3).


The predicted recovery from AR/ASR was estimated by simulating recharge to each model cell representing an aquifer and computing the resulting change in aquifer system storage. Because the faults are strong flow barriers and the goal is to restore groundwater levels in the OWRD groundwater management area, AR/ASR was simulated into each CRBG aquifer model cell south of the Rocky Prairie thrust fault and east of the Maupin wrench fault (fig. 4). For all AR/ASR simulations, simulated injection rate was 250 gal/min for the 6 winter months, which corresponds to approximately 202 acre-ft of water each year or 21.5 percent of total annual pumping for the entire study area. Storage change resulting from AR/ASR in each model cell may be plotted for each basalt aquifer (for example, see fig. 30). Assuming AR/ASR would be at the optimal location, the maximum computed change in storage from all cells in all layers is reported in table 3 for comparison. The “Value of AR/ASR” column (table 3) shows the additional benefit of conducting AR/ASR in conjunction with some amount of well repair.


As commingling wells are repaired, the value of AR/ASR generally improves because less water is lost through intra-borehole flow. Additionally, the location of greatest change in storage varies depending on which wells have been repaired (fig. 30). Without well repair, injection in the OWRD administrative area yields little improvement in aquifer storage, but after repairing wells in the area, it becomes the best location to enhance aquifer-system storage using AR/ASR. 


In all cases, simulations showed that application of 
AR/ASR to the Frenchman Springs aquifer would provide the most benefit to system storage overall owing to its large extent and the fact that commingling waters leaking from this unit will flow upward through shallower units. The maximum benefit location for the “No Repair” scenario is predicted to be near the crest of the Columbia Hills anticline (fig. 30A). Note, however, that the model boundary in this location prevents water leakage out of the model to the south. In reality, the overlying limb of the Chenoweth thrust fault to the south may allow water to leak out, indicating that the estimated 5.9 percent maximum recovery is likely an overestimation of the potential benefit of AR/ASR for the “No Repair” scenario.


Evaluation of the Value of Repairs Targeting a Single Confining Unit


The final scenario considered was to repair the integrity of a single confining unit. In practice, a well open to three or more aquifers separated by confining units would not be repaired at a single confining unit, but the simulation results illustrate aquifer response that may be used when designing repair strategies. For example, in an area with multiple confining units, if one confining unit is penetrated by fewer wells than the other confining units, repair of this subset of wells effectively would hydraulically isolate the aquifers above and below the repaired confining unit, potentially providing substantial benefit. For the simulations, all commingling wells in repair zones 1, 2, and 3 were partially repaired. In each case, the repair separates two adjacent CRBG aquifers, but all other commingling still is assumed to exist.


Of the 82 wells simulated as commingling in these zones, at least one-half this number penetrates each confining unit (table 4). The percent recovery attributed to AR/ASR is fairly constant, but the percent recovery attributed to repairs starts out high for repair of the upper two layers, decreases for the third, but then increases again for the fourth. This irregular pattern can be explained with a simple conceptual model of storage.


The single commingling well problem used to generate the vulnerability maps is conceptually illustrated by figure 31. Because water is added to each aquifer by recharge, if the commingling well did not exist, groundwater levels would rise until they reach the point where they naturally spill over (fig. 31A). The physical spill point for aquifers in the study area is where an aquifer intersects a spring or creek. Repair of any single confining unit allows water to fill any cross‑connected aquifers until the lowest spill point is reached (fig. 31C through F). The amount of recovery depends on the geometry of the aquifer being filled.


In the Mosier aquifer system, there are many cross‑connecting points (commingling wells), so complete repair of all wells in a zone (for example, fig. 29) will only provide partial recovery of the system (fig. 32). A more complex distribution of commingling wells (fig. 33) provides a simple explanation for the simulation results (table 4) that relies on the locations of commingling wells instead of the size and geometry of the aquifers. 


In conclusion, repair of wells that raise the effective outfall for one or more aquifers will raise water levels in only those aquifers. Because the aquifer system is sloped, this indicates that sealing the system from the low end to the high end will provide a systematic improvement in water levels in all aquifers. Because the system currently is not in equilibrium, this strategy may also slow the rate of declines, because the gradient driving declines will be diminished.


Limitations of Management Option Analysis


A single set of model parameters was used for all management scenarios. This allowed a rapid comparative analysis of different management options, but the uncertainty in the magnitude of benefit has not been evaluated. Prior to adopting a particular management strategy, a predictive uncertainty analysis of the range of expected benefits could be performed.


First posted March 1, 2012