- Email: [email protected]
An unconfined groundwater model of the Death Valley Regional Flow System and a comparison to its confined predecessor
Journal of Hydrology 373 (2009) 316–328
Contents lists available at ScienceDirect
Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol
An unconfined groundwater model of the Death Valley Regional Flow System and a comparison to its confined predecessor Rosemary W.H. Carroll *, Greg M. Pohll, Ronald L. Hershey Division of Hydrologic Sciences, Desert Research Institute, 2215 Raggio Parkway, Reno, NV 89512, United States
a r t i c l e
i n f o
Article history: Received 12 November 2008 Received in revised form 21 April 2009 Accepted 4 May 2009
This manuscript was handled by P. Baveye, Editor-in-Chief, with the assistance of Hongbin Zhan, Associate Editor Keywords: Death Valley Nevada Test Site MODFLOW SURFACT
s u m m a r y The MODFLOW version of the United States Geological Survey (USGS) Death Valley Regional Flow System (DVRFS) in California and Nevada is conceptually inaccurate in that it models an unconfined aquifer as a confined system and does not accurately simulate unconfined drawdown in transient pumping simulations. The transfer of geologic and hydrologic information from the confined MODFLOW DVRFS model to an unconfined MODFLOW–SURFACT (SURFACT) version was accomplished by maintaining cell structure between models and computing effective cell properties to translate the HUF2 package used in MODFLOW to the BCF4 package used by SURFACT. The confined version of the DVRFS was compared to the unconfined SURFACT version by examining head contour maps and the ability of the SURFACT model to match the 4900 observations of hydraulic head/drawdown, 49 observations of groundwater discharge, and 15 estimates of groundwater fluxes into/out of the model domain. Resultant weighted root mean squared error (xRMSE) for the unconfined SURFACT model was lower than the USGS confined model. Despite a lower xRMSE, unconfined conditions simulated with SURFACT did produce greater heads in mountainous regions compared to the confined MODFLOW with differences most pronounced in regions where cell thickness is large, horizontal conductivity small and recharge large. Difference in computed heads reflects computation schemes employed by both models to estimate interblock conductance. Specifically, interblock conductance for the unconfined SURFACT model is dependent on the relative saturation of a modeled cell while MODFLOW’s confined system is not. Despite head differences, SURFACT simulates comparable flux estimates to MODFLOW (e.g. observed ET, groundwater spring flow, and groundwater flux across model boundaries), while significantly improving transient well drawdown estimates. SURFACT is also capable of producing more realistic estimates of water availability from proposed groundwater development and resultant potential impacts to the region. Ó 2009 Elsevier B.V. All rights reserved.
Introduction There have been numerous water rights applications submitted throughout southern Nevada to offset the needs of growth in Las Vegas. Several applications include groundwater withdrawals adjacent to the Nevada Test Site (NTS) where underground nuclear testing was conducted. If large quantities of groundwater are pumped adjacent to the NTS, the groundwater system could change dramatically. Potential impacts from groundwater pumping include decreasing water levels, reduction in groundwater resources on the NTS, reduction in spring flows adjacent to proposed pumping centers, and the alteration of groundwater flowpaths.
* Corresponding author. Address: Division of Hydrologic Sciences, Desert Research Institute, 2215 Raggio Parkway, Reno, NV 89512 (P.O. Box 4196, Crested Butte, CO 81224), United States. Tel.: +1 970 349 0356; fax: +1 775 673 7363. E-mail addresses: [email protected] (R.W.H. Carroll), [email protected] (G.M. Pohll), [email protected] (R.L. Hershey). 0022-1694/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2009.05.006
In 1998, the US Department of Energy (DOE) Nevada Site Office funded the United States Geological Survey (USGS) to improve upon two previous groundwater flow models of the Death Valley region with the initial intention of understanding groundwater flowpaths and travel times associated with potential movement of radioactive material from the NTS as well as to characterize the groundwater system in the vicinity of Yucca Mountain, and address effects on users down-gradient from the NTS and Yucca Mountain (Belcher, 2004). The first of these earlier models was developed by DOE for the National Nuclear Security Administration/Nevada Site Office (NNSA/NSO) Underground Test Area (UGTA) project (IT Corporation, 1996). The second was developed by the Office of Civilian Radioactive Waste Management’s (OCRWM) Yucca Mountain Project (YMP) and the NNSA/NSO Hydrologic Resources Management Program (HRMP). The resultant USGS Death Valley Regional Groundwater Flow System Model (DVRFSM) has improved upon the two previous models by using newly acquired data and modeling tools (Belcher, 2004). Studies compiled and used in the DVRFSM construction
317
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
include reassessing groundwater discharge via evapotranspiration (ET) (Laczniak et al., 1999, 2001; Reiner et al., 2002; DeMeo et al., 2003) and spring flow (refer to Table C-2, Belcher, 2004), cataloging historical groundwater pumping from 1913 through 1998 (Moreo et al., 2003), reinterpreting groundwater recharge as net infiltration (Hevesi et al., 2002, 2003), assessing model boundary inflows and outflows from regional hydraulic gradients, developing a water budget (Belcher, 2004), and finally incorporating hydraulic conductivity relationships as a function of depth (Belcher et al., 2001, 2002). Carroll et al. (2006) used the DVRFSM to test the impacts of proposed groundwater development south of the NTS and found that, despite increased level of geologic detail and improved water budget accounting, the DVRFSM produced drawdowns on the order of thousands of meters after only a few decades of pumping. Such large drawdowns are not realistic and serve to highlight the most significant limitation of the DVRFSM: the model was built using the confined layer assumption to improve numeric stability. The use of a confined, rather than a ‘‘convertible” (i.e., unconfined), layer type within MODFLOW assumes that the saturated thickness remains constant throughout the entire simulation. Therefore, cells are not allowed to dry or become inactive as the water level decreases below the bottom of a cell. The confined approach was adopted by the USGS because the model was computationally unstable when cells were allowed to convert between confined
and unconfined. Unfortunately the confined approach, while numerically stable, allows unrealistic estimates of drawdown to occur in transient simulations. The primary objective of this study was to produce a numerically stable, unconfined version of the Death Valley Regional Flow System (DVRFS) for future use as a tool to more accurately address potential withdrawals in the NTS region and to determine the potential differences in the solution of the groundwater flow equation given confined and unconfined assumptions.
Numerical models Death Valley Regional Flow System model: MODFLOW The DVRFS and the associated USGS model domain (DVRFSM), with the NTS superimposed, are shown in Fig. 1. The DVRFS is approximately 100,000 km2 in Nevada and California and is bounded by latitudes 35°000 N and 38°150 N and by longitudes 115°000 W and 118°000 W. This system encompasses flow between recharge areas in the mountains of central and southern Nevada and discharge areas of wet playas and springs south and west of the NTS and in Death Valley, California. The flow is strongly influenced by a complex geologic framework and the USGS DVRFSM incorporates the distribution of the flow system’s principal aqui-
Death Valley Region Flow System and Model Boundary
Nevada California
A
A’
0
45 km
Fig. 1. Location of the Death Valley groundwater flow system with the Nevada Test Site (NTS), Death Valley, Spring Mountains, and the Nevada–California state line marked. Inset shows the site with respect to the western United States. Select mountain ranges labeled with white text and select basins labeled in black text. Cross section A–A0 marked.
318
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
NTS
Hypothetical Well
N
0
45 km
Fig. 2. Cell structure in the DVRFSM with the location of a hypothetical well (used for pumping tests) marked.
fers and confining units as well as the principal geologic structures that might affect subsurface flow (Belcher, 2004). The DVRFSM was built with MODFLOW-2000 (McDonald and Harbaugh, 1988; Harbaugh et al., 2000), a public domain, finite-difference code developed by the USGS that is capable of simulating groundwater flow in transient, three-dimensional, anisotropic and heterogeneous systems. MODFLOW’s governing three-dimensional flow equation for a confined aquifer combines Darcy’s Law and the principle of conservation of mass via
@ @h @ @h @ @h S @h K xx þ K yy þ K zz G¼ @x @x @y @y @z @Z b @t
ð1Þ
where x, y, and z are Cartesian coordinates (L), Kxx, Kyy and Kzz are the principal components of saturated hydraulic conductivity along the x, y, and z axes, respectively (L/T), h is hydraulic head (L), G is a volumetric source/sink term (1/T), S is stortivity (dimensionless), b is saturated thickness (L) which is equivalent to cell thickness, and t is time (T). Conversely, Eq. (2) defines groundwater flow given unconfined conditions, where h (L) is the phreatic surface elevation and Sy is specific yield (dimensionless).
@ @h @ @h @ @h @h K xx h þ K yy h þ K zz h G ¼ Sy @x @x @y @y @z @z @t
ð2Þ
Solving for h in Eq. (2) is non-linear problem because h depends on hydraulic conductivity and saturated thickness, with saturated thickness dependent on the water table’s elevation. To address
non-linearity inherent within unconfined aquifer systems, MODFLOW’s updated HUF2 package maintains a confined system, with the top elevation equal to the approximate elevation of the water table, and substitutes Sy for S into Eq. (1). The HUF2 approach to simplifying the non-linearity of an unconfined solution is only applicable when there is little change in saturated thickness. Belcher (2004) describes the DVRFSM in detail, however a consolidated review of the DVRFSM is provided here for clarity. The DVRFSM uses a constant-grid spacing in the horizontal plane of 1500 m by 1500 m. The finite-difference grid is oriented north– south, with 160 columns, 194 rows, and 16 layers for a total of 314,784 active cells (Fig. 2). Cell thicknesses range between 1 m and more than 3000 m, with layer 16 (the bottom layer) reaching 4000 m below sea level. The top elevation of layer one represents the water table (and modeled as a potentiometric surface) and not the actual land surface elevation with cell thickness in layer one increasing substantially in mountain ranges (Fig. 3). Fig. 4 shows the model layer configuration along cross section A–A0 (refer to Fig. 1 for cross section location) in comparison to land surface elevation and simulated hydrogeologic units. Twenty-seven hydrogeologic units were modeled in the DVRFSM via the MODFLOW2000 HUF2 package. A brief description of these units is provided in Table 1. In general, the model’s upper layers portray relatively shallow basin-fill sediments and volcanics, while deeper layers define the regional carbonate aquifer and confining units. Model cells containing more than one hydrogeologic unit are simulated with vertically averaged hydraulic properties.
319
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
Cell Thickness (m) 150 – 500 100 – 150 N
50 – 100
0
1 - 50
45 km
Fig. 3. Cell thickness (m) in layer one.
Model calibration (Belcher, 2004) for steady state conditions prior to year 1912 and for transient conditions for years 1913– 1998 was accomplished by best matching 4899 observations of hydraulic head and drawdown as well as 49 estimates of groundwater discharge and 15 estimates of regional flux into/out of the model domain using MODFLOW-2000 parameter estimation processes (Hill et al., 2000). One hundred parameters were estimated in the calibration process including hydraulic conductivity, vertical anisotropy, drain conductance, fault-associated hydraulic conductivity, recharge zone multipliers, depth decay coefficients for hydraulic conductivity, and storage parameters (S and Sy). Fig. 5 provides an example of where boundary conditions, applied stresses (e.g. wells, drains) and faults (i.e., hydrologic flow barriers) are designated in the DVRFSM’s top layer during the final modeled timestep (year 1998). Fig. 6 shows recharge rates (m/d) applied throughout the modeled domain. Recharge is highest in mountainous regions, in particular the Spring Mountains contribute the greatest volume. Recharge was not modeled in most lowland areas, but those with designated recharge represent irrigation recharge or the re-infiltration of groundwater springs (e.g. Death Valley). Lateral flows into and out of the model boundary were specified as constant head cells that were interpolated from the regional potentiometric surface. Only the region bounding the Spring Mountains was defined as a no flow boundary. Drains were used to define both groundwater ET and spring flows. Discharge was defined by observation and drain conductances were determined by parameter estimation while drain elevations were
set equal to 10 m below the lowest land-surface altitudes for each group of drain cells modeled. These elevations were assumed to adequately represent extinction depths for ET as well as account for springs being located in land-surface depressions (Belcher, 2004). Faults acting as barriers to flow were defined in MODFLOW-2000 with the Horizontal-Flow Barrier (HFB) package. Faults were assumed to extend the entire thickness of the model and the width of each fault was assumed 1 m, while the hydraulic conductivity of each fault was determined by parameter estimation. Nine faults (i.e., horizontal flow barriers) were found to exert some influence on hydraulic heads and are shown in Fig. 5. The DVRFSM simulated hydraulic heads in layer one for the final time step (year 1998) are shown in Fig. 7. SURFACT SURFACT (HGL, 2007) is a fully integrated flow and transport numeric code based on the USGS finite-difference groundwater model MODFLOW. SURFACT potentially improves upon MODFLOW with additional modules. Specifically, SURFACT’s modified BCF4 package performs a complex saturated-unsaturated subsurface flow analysis using a pseudo-soil function for a rigorous treatment of unconfined flow. This is done using robust numerical and matrix solution techniques for greater stability (PCG4) and an adaptive time-stepping scheme (ATO) to improve model efficiency. The three-dimensional movement of water in variably saturated media is given by (Huyakorn et al., 1986),
320
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
(a)
4,000 A Top of Layer 1 = simulated potentiometric surface
3,000 Panamint Mountains Land Surface
Spring Mountains
A'
2,000
Elevation (m)
Greenwater-Black Mountains
1,000
Death Valley
Amargosa Desert
0 -1,000 -2,000 -3,000
Bottom of Layer 16
-4,000
(b)
4,000 0 A 3,000
Elevation (m)
2,000
20
40
60
XCU (K1) LCA (K2) LCA_T1 (K2) VSU_UPPER (K4) MISC Basin fill (K4)
80
100
120
LCCU (K1) LCCU_T1 (K1) VSU_LOWER (K4) LFU (K3)
A'
1,000 0 -1,000 -2,000 -3,000 -4,000 0
30
60
90
120
Distance (km) Fig. 4. Cross section A–A0 (a) configuration of 16 model layers and (b) hydrogeologic units along cross section with units and major rock type (K1, K2, K3 and K4) defined in Table 1.
@ @h @ @h @ @h K xx krw þ K yy krw þ K zz krw G @x @x @y @y @z @z ¼U
@S @h þ Sw SS @t @t
ð3Þ
where additional terms relative to Eq. (1) include krw equal to relative permeability (a function of water saturation, dimensionless), u equal to a drainable porosity (i.e., specific yield Sy), and Sw representing the degree of water saturation (a function of pressure head, dimensionless). When fully saturated conditions exist, then Sw = 1.0, krw = 1.0, and Eq. (3) reduces to Eq. (1), the confined groundwater equation used by MODFLOW. Consequently, the conventional flow equation is used by SURFACT below the water table and in confined systems. While MODFLOW’s other flow packages (BCF and LPF) support several averaging schemes (e.g. harmonic mean, arithmetic mean, arithmetic mean thickness-logarithmic mean hydraulic conductivity) for computing effective interblock transmissivity, MODFLOW’s HUF2 package is limited to only the harmonic mean (i.e., harmonic mean of transmissivity between cell i and cell j is defined as Ti,j = 2TiTj/(Ti + Tj)). The effective transmissivity is equal to the saturated hydraulic conductivity multiplied by the saturated block thickness. For confined systems, the saturated thickness is equal to cell thickness (b). The harmonic mean, by definition, removes the impacts of large outliers, aggregates smaller values and is best used to average
fluxes. However, its estimation automatically underestimates the equivalent interblock conductivity (Romeu and Noetinger, 1995) by biasing toward the lower block value, and in the extreme case never allows a desaturated block (i.e saturated thickness of zero) to resaturated. In contrast, SURFACT’s BCF4 package computes interblock conductance as a product of the weighted harmonic mean of block saturated hydraulic conductivities (Kxx, Kyy, Kzz), relative permeability (krw), and mean flow area. The SURFACT approach to interblock conductance is more stable compared to MODFLOW’s unconfined solution, does not bias toward the lower block value and does not inactivate dry cells during desaturation (HGL, 2007). SURFACT provides a fracture well (FWL) technique for more realistic simulation of groundwater pumping withdrawals compared to the original MODFLOW well (WEL) package. The original WEL package required the user to prescribe the contribution of water from each node (i.e., layer) in a multi-node well and did not automatically adjust pumping rates when water levels dropped below screened intervals. Instead, pumping produced withdrawals into the nonphysical realm. On the other hand, FWL automatically computes (and updates) relative nodal contributions of total well flux rates. For an unconfined system, the FWL package adjusts water withdrawal to feasible volumes from a potentially over-pumped well such that heads do not drop below a well’s screened interval.
321
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328 Table 1 The 27 hydrogeologic units modeled with HUF2 package in the DVRFM and corresponding major rock types. Table modified from Belcher (2004) page 272 Table F-3. Major rock type
Abbreviation
Name
Confining units (K1)
LCCU LCCU_T1 UCCU XCU ICU SCU
Lower clastic-rock confining unit Lower clastic-rock thrust Upper clastic-rock confining unit Crystalline-rock confining unit Intrusive-rock confining unit Sedimentary-rock confining unit
Carbonate rock aquifer (K2)
LCA LCA_T1 UCA
Lower carbonate-rock aquifer Lower carbonate-rock thrust Upper carbonate-rock aquifer
Volcanic rock units (K3)
LFU OVU PVA TMVA CHVU CFTA CFBCU CFPPA YVU WVU BRU
Lava-flow unit Older volcanic-rock unit Paintbrush volcanic-rock aquifer Thirsty canyon-timber mountain volcanic-rock aquifer Calico Hills volcanic-rock unit Crater flat-tram aquifer Crater flat-bullfrog confining unit Crater flat-prow pass aquifer Younger volcanic-rock unit Wahmonie volcanic-rock unit Belted range unit
Basin fill sediments (K4)
OACU YACU YAA LA OAA VSU_LOWER VSU_UPPER
Older alluvial confining unit Younger alluvial confining unit Younger alluvial aquifer Limestone aquifer Older alluvial aquifer Lower volcanic and sedimentary-rock unit Upper volcanic and sedimentary-rock unit
Fig. 5. An example of designated boundary conditions and cell stresses within the DVRFSM (year 1998, top layer).
322
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
Recharge (m/d) 1e-5 – 1e-3 1e-7 – 1e-3 1e-8 – 1e-5
N
0
0 – 1e-8
45 km
Fig. 6. Recharge (m/d) applied to the model’s top layer.
Methods Translation of the DVRFS from MODFLOW to SURFACT Translation of MODFLOW’s geological model packages HUF2 and KDEP (Anderman and Hill, 2000; 2003), which accounts for hydraulic conductivity as a function of depth, into SURFACT’s modified BCF4 package was accomplished by extracting effective cell properties from the original DVRFSM (e.g. horizontal hydraulic conductivity, vertical hydraulic conductivity, and specific storage) using a Fortran code developed for this purpose and compiled as a subroutine of MODFLOW. Specific yield in the DVRFSM was originally calibrated between 0.19 and 0.20 depending on the hydrogeologic unit considered. A value of 0.19 was assumed reasonable and assigned across the entire SURFACT model domain. Fig. 8 shows hydraulic conductivity (K) in the model’s top layer as an example of effective values used by SURFACT. Comparison between MODFLOW and SURFACT The SURFACT DVRFS model was run using USGS-calibrated parameters for the MODFLOW DVRFSM. Comparison between the unconfined and confined modeling approaches included: (1) areal head distributions at the end of the transient model (year 1998), (2) ability to predict observed values of heads, drawdown, groundwater fluxes to springs (i.e., modeled drains) and across modeled boundaries (i.e., constant heads), and (3) estimated drawdowns related to 50 years additional pumping and 50 years recovery begin-
ning in year 1998 from a hypothetical well south of the NTS, a region often considered for groundwater development to support Las Vegas growing need for water. Results and discussion Predicted heads Computed heads at the end of the transient simulation (year 1998) found both MODFLOW and SURFACT consistent in their estimates of flow magnitude and direction. However, differences in predicted heads did occur. In general, the unconfined modeled heads were higher than those predicted by the confined model in mountainous regions with greatest differences occurring in mountain ranges defined with thick model cells (b), low hydraulic conductivity (K) and high recharge flux (w). Specifically, SURFACT predicted heads in excess of 100 m compared to MODFLOW in the Spring Mountains, the Sheep Range, Pintwater Range, Timpahute Range, and Last Chance Range. SURFACT predicted heads on the order of 10–100 m above MODFLOW predicted heads occurred within the Amargosa Range, Black Mountains, Greenwater Range, Belted Range, Black Mountain and the Yucca Mountain region. These regions are defined with thinner cell thicknesses, higher hydraulic conductivities and lower recharge values than other mountain regions. SURFACT and MODFLOW predicted heads were relatively similar for most lowland regions in the model, with SURFACT predicted heads less than MODFLOW in Cactus Flat, Penoyer Valley, Tikaboo Valley, the southern end of the Amargosa Desert
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
323
Fig. 7. DVRFSM-simulated heads (m) for the model’s top layer given 1998 pumping conditions. Contour interval is 200 m.
(i.e., Amargosa Flat) and the northern end of Death Valley. Differences were on the order of 1.0–10.0 m. To understand discrepancies between modeling approaches it is necessary to revisit the groundwater flow equation. The analytical solution for a confined aquifer defined by Eq. (1) and assuming 1dimensional, steady state flow given homogenous, isotropic conditions is,
h ¼ ho
wx2 h1 h0 2L x þ þ Kb 2Kb L
ð4Þ
where h0 and h1 equal heads (L) at x0 = 0 and x1 = L, respectively and w is recharge (L/T). Assuming h0 = h1 = 0 and that the groundwater divide ðq0x ¼ Kbðdh=dxÞ ¼ 0Þ occurs at x = L/2 will produce a maximum head (hmax) equal to,
hmax;c ¼
wL2 8Kb
ð5Þ
Given an unconfined aquifer with identical conditions and using Dupuit assumptions, the 1-dimensional, analytical solution becomes (Fetter, 1994),
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ðh h1 Þx w 2 h ¼ h0 þ 0 þ ðL xÞx L K with hmax at a groundwater divide x = L/2 equal to,
ð6Þ
hmax;uc
sffiffiffiffiffiffiffiffiffi wL2 ¼ 4K
ð7Þ
Fig. 9a compares the analytical difference between hmax,uc and hmax,c for various ratios of w/K and b/L. At low values of w/K and for all values of b/L, hmax,uc > hmax,c differences in solutions are small. Differences become slightly larger with increasing values of b/L. With increased values of w/K and low ratios of b/L, the relationship between unconfined and confined changes such that hmax,uc < hmax,c with differences in heads becoming increasingly significant at larger and larger w/K. In contrast, unconfined solutions for large b/L surpass confined heads given all examined quantities of w/K, with substantial differences occurring at large w/K. Fig. 9b and c compare 1- and 3-dimensional numeric results, respectively, for unconfined heads computed with SURFACT (hmax,SF) and confined heads computed with MODFLOW (hmax,MF). Cell length in all simulations was held constant (L = 1500 m) while cell thickness (b) was allowed to vary. Recharge (w) and vertical hydraulic conductivity (Kzz) were held constant while Kxx, Txx and MODFLOW’s Vcond (Kzz/b) were adjusted accordingly. Layer one bottom elevations were defined at the same elevation as h0 and h1 (i.e., 0 m). This simplified the analysis since hmax,SF = 0 implies a krw = 0, and an hmax,SF P b implies krw = 1.0. Differences between analytic and numeric results reflect computation schemes employed by both models to estimate interblock
324
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
Fig. 8. Effective hydraulic conductivity (K) in the model’s top layer.
conductance and the relationship of predicted hmax,SF to the top elevation of the modeled cell (in this example top elevation of layer 1 = b). MODFLOW’s confined system uses a constant saturated thickness and a constant saturated hydraulic conductivity. SURFACT’s unconfined approach, however, allows for variable saturation of the modeled cell. For the situation where hmax,SF < b, then krw < 1.0 and the effective hydraulic conductivity used by SURFACT is less than the pre-prescribed saturated hydraulic conductivity used by MODFLOW (both in x- and/or z- directions). Thereby, SURFACT produces greater heads than MODFLOW, with greatest differences observed at higher w/K. The difference between hmax,SF and hmax,MF is smaller than its counterpart in the analytical solution (hmax,uchmax,c) because of MODFLOW’s use of transmissivities and subsequent under estimation of interblock conductance. Fig. 9b, b/L = 1.0 illustrates this point. For the situation where hmax,SF P b, then SURFACT uses krw = 1.0 and the flow area between cells is limited to the cell thickness, even if predicted heads lie above the top elevation of the modeled cell. Therefore, the unconfined and confined systems estimate equivalent cell transmissivities. Given these circumstances, unconfined heads predicted by SURFACT begin to converge at increasing values of w/K toward solutions in which MODFLOW heads are slightly larger than SURFACT heads. Despite having numerically equal transmissivities, heads do not converge toward an equal solution because SURFACT computes interblock conductance with
hydraulic conductivity, not transmissivity, and conductance is less biased toward low conductance compared to MODFLOW. Convergence toward nearly equal solutions at large w/K was not witnessed in the analytical solution because the area of flow conveyance for the unconfined system was not limited by b, but allowed to increase to the predicted value hmax,uc. The resultant larger conductance forced the unconfined system to predict significantly smaller heads compared to the confined solution. For the 1-dimensional numeric model, conductance was limited to the horizontal direction (Kxx), while the 3-dimensional numeric model included vertical conductance (Kzz). Resultant leakance from the top layer to deeper layers caused heads in the top layer of the 3-dimensional model to be lower than those in the 1-dimensional model for any given value of w/K. The impact was that SURFACT heads remained less than cell thickness for larger and larger values of w/K. The principles discussed above were applied to the cross section A–A0 (refer to Fig. 1), which transects several different mountain ranges and low lying basins. Recharge, hydraulic conductivity, land surface elevation, cell top and bottom elevations as well as modeled heads were extracted from layer one on a cell-by-cell basis. Quantities were plotted in Fig. 10. Included is a plot of relative cell saturation for unconfined conditions, (hSF cell bottom elevation)/b, (refer to Fig. 10d). A constant head boundary condition forced MODFLOW and SURFACT to predict the same heads at location A,
325
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328 600
400
Elevation (m)
500
b/L = 1
300
h max,uc - h max,c
200 100 0
b/L = 0.1
-100
4,000 A 3,500 3,000 2,500 Panamint Mountains Greenwater-Black 2,000 Ranges Land Surface 1,500 Amargosa Desert 1,000 Death 500 Valley 0 MODFLOW SURFACT -500 020406080100120140 0.20
Distance (km)
A
-200
Spring Mountains
A'
A'
0.15
-300
b/L
b/L = 0.01
-400 -500 -600 0.0001
0.10 0.05
0.001
0.01
0.1
0.00
1
w/K
020406080100120140
1.E+00
w/K
400
b/L = 1
300
A'
Distance (km)
A'
1.E-04 1.E-05 1.E-06
200 100
1.E-03
0.0 1.E-07
b/L = 0.1
020406080100120140
b/L = 0.01
-100 -200 -300 -400 -500
0.001
0.01
0.1
6.00 4.00 2.00
-2.00
h SF-h MF (m)
400
75
b/L = 1
50
0.0
0.00
500 020406080100120140
1
w/K 100
A
Krw = 0.0
Relative Saturation
0
Krw = 1.0
h max,SF - h max,MF
Distance (km)
1.E-02
500
-600 0.0001
A
1.E-01
600
A
A'
Distance (km)
300 200 100 0
h max,SF - h max,MF
-100
25
0
b/L = 0.1
20
40
60
80
100
120
140
Distance (km)
0
b/L = 0.01
Fig. 10. Cross section A–A0 layer one cell values at the end of the transient simulation (year 1998) for: (a) land surface elevation, MODFLOW predicted heads (hMF) and SURFACT (hSF) predicted heads. The shaded area depicts distance between top and bottom elevation of the cell, (b) b/L, (c) w/L, (d) relative saturation and (e) hSFhMF.
-25 -50 -75 -100 0.0001
0.001
0.01
0.1
1
w/K Fig. 9. Difference in steady state, maximum head versus w/K for various values of b/ L: (a) 1-dimensional analytical solution for unconfined (hmax,uc) and confined (hmax,c) systems, (b) 1-dimensional numeric solution for unconfined SURFACT (hmax,SF) and confined MODFLOW (hmax,MF), (c) 3-dimensional numeric solutions for unconfined SURFACT and confined MODFLOW.
while a no-flow boundary allowed heads to vary at location A0 . Low values of w/K kept MODFLOW and SURFACT heads nearly equivalent throughout the Penamint range. SURFACT heads exceeded MODFLOW heads only in regions where SURFACT’s relative saturation was less than 1.0, thereby reducing krw and the equivalent SURFACT conductivity in relationship to MODFLOW. More significant deviation occurred along the western edge of the Greenwater and Black Ranges. In this region, SURFACT produced a krw < 1.0 and
Table 2 A comparison of MODFLOW and SURFACT weighted root mean squared errors (xRMSE) for each type of observation.
xRMSE Observation
Number
MODFLOW
SURFACT
Heads Drawdown Drains (i.e., ET) Constant heads Total
1930 2969 49 15 4963
100.19 18.95 3.93 5.40 64.17
100.82 14.06 10.36 4.94 63.89
values of w/K were large enough to produce SURFACT heads 10– 50 m above MODFLOW heads. The greatest difference occurred at a distance of 30.0 km, when relative saturation dropped below zero (heads below top layer), pegging krw to a value of 0.0. In the Spring Mountains, the ratio of w/K was much higher than in the Green-
326
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
water-Black Range forcing SURFACT heads to exceed the top elevation of layer one despite substantial values of b/L. None the less, SURFACT heads were in excess of 400 m compared to MODFLOW heads. Given krw = 1.0, it is presumed that SURFACT heads would have been larger if the top elevation of the cell was raised to land surface elevations and boundary conditions maintained. For the unconfined system, placing layer one’s top elevation at land surface would be preferable, however the USGS calibrated depth decay coefficients in the HUF2 package were based on layer one’s current top elevation. Redefining the top elevation to land surface would negate this calibration. SURFACT and MODFLOW predicted heads were nearly equivalent in low lying areas along A–A0 since w/K was generally low and krw was approximately 1.0.
where Xo is the observed value, Xp is the predicted value, n is the number of observations, and x is the inverse variance. Predicted values derived from multiple layers were computed using fractional inputs from each contributing layer as defined by Belcher (2004) and discussed by Hill (1998) and Hill et al. (2000). A comparison of xRMSE values is provided in Table 2 and also shown in Fig. 11a for heads, Fig. 11b for drawdown, Fig. 11c for drain fluxes and Fig. 11d for constant head fluxes. Both models perform similarly with the SURFACT model having a slightly smaller total xRMSE. In addition the xRMSE for drawdown was significantly smaller than the MODFLOW model. The larger xRMSE for drawdown in the MODFLOW model is the result of extremely large drawdown values that were simulated by MODFLOW using the confined governing equation.
Observations Well withdrawal Model results were compared to observed data using the weighted root mean squared error (xRMSE),
Pn i
xRMSE ¼
(a)
xi ðX 0i X pi Þ2
!0:5 ð8Þ
n 5,500
(b)
MODFLOW
300
SURFACT
250
200
Modeled Drawdown (m)
4,500
Modeled Heads (m)
Last, a hypothetical well was placed in the model’s top layer to the south of the NTS (row 124, column 99, refer to Fig. 1). This well’s location was chosen for its proximity to the NTS and because predicted heads by the MODFLOW and SURFACT models were nearly
3,500
2,500
1,500
150 100 50
0
500 MODFLOW
-50
SURFACT
-500 -500
1,500
500
2,500
3,500
-100 -100
4,500
Observed Heads (m)
(c)
-50
0
50
100
150
200
250
300
Observed Drawdown (m)
(d)
0
90,000
MODFLOW
MODFLOW
-5,000
SURFACT
-10,000
3
-15,000 -20,000 -25,000 -30,000 -35,000
40,000
-10,000
-60,000
90,000
40,000
-10,000
0
-5,000
-10,000
-15,000
-20,000
-25,000
-30,000
-35,000
-40,000
Observed Drains(m 3/d)
-60,000
-110,000
-40,000
-110,000
Modeled Drains (m /d)
3
Modeled Constant Head (m /d)
SURFACT
Observed Constant Head (m 3/d)
Fig. 11. Comparison of MODFLOW- and SURFACT-predicted values to observed values used in the calibration of the MODFLOW model: (a) heads, (b) drawdowns, (c) drain fluxes, and (d) constant head fluxes.
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
not realistically represent the actual unconfined system, which cannot support drawdown of heads below the well unless an adjacent and deeper well is pumped simultaneously.
MODFLOW
785
SURFACT
Head in Well (m)
760 735
Conclusions
Bottom Elevation of Well
710 685 660 Pumping Ends
635
(a)
610 585 0
10
20
30
40
50
60
70
80
90
100
70
80
90
100
Years of Pumping 6,000
5,000
Well Flux (m 3/d)
327
4,000
Prescribed
3,000
2,000 Actual
1,000
(b)
0 0
10
20
30
40
50
60
Years of Pumping Fig. 12. (a) Comparison of head drawdown in a hypothetical well south of the NTS between the MODFLOW (confined) and SURFACT (unconfined) models given 50 years pumping at 5000 m3/d (1480 acre ft/yr) followed by 50 years of recovery. (b) Prescribed pumping rate compared with the SURFACT model actual withdrawal.
identical at the end of transient simulation (year 1998), which simplifies the analysis of subsequent drawdown with respect to the bottom of the well. The bottom of the well was established at the bottom elevation of the cell containing the well (729 m above mean sea level). A pumping rate of 5000 m3/d (1480 acre ft/yr) was initiated at the end of the transient model and then allowed to recover for an additional 50 years with no pumping. Results (Fig. 12) show that the SURFACT model can only meet prescribed pumping rates during the first 10 years of pumping, after which drawdown approached the base of the well and pumping rates were reduced to feasible quantities that the aquifer could supply to maintain the water table at the base of the well. Maximum modeled drawdown in the well reached only 60.6 m, or the bottom elevation of the well. After pumping ceased, slow rebound in the water table resumed such that, after 50 years, water in the well rebounded approximately 5.5 m to an elevation of 734.5 m. Unlike the SURFACT model, MODFLOWs confined conditions allowed full withdrawal demand to be met for the entire simulation. Divergence in SURFACT and MODFLOW simulated drawdown began after 10 years when the SURFACT model reduced the amount of stress on the well to keep heads just above the bottom elevation of the well while the MODFLOW model maintained the prescribed stress. At approximately 15 years of pumping, the MODFLOW model’s predicted heads dropped below the bottom of the well, and after 50 years of pumping, predicted heads fell to 586.4 m, or more than 142 m below the bottom of the well. While large drawdowns may be numerically possible in a confined aquifer, they do
Representation of the DVRFS using SURFACT is an improvement over its MODFLOW predecessor for transient simulations. The original USGS model assumed that the entire model domain was in a confined state to ensure numerical stability. The uppermost layer used a specific yield storage parameter to represent unconfined conditions, but the cell thickness does not change in response to water level changes. The result was projected drawdown estimates from well withdrawals were unrealistic. On the other hand, SURFACT’s saturated-unsaturated subsurface flow analysis using a pseudo-soil function produced a numerically stable model of the DVRFS using unconfined conditions. In addition, SURFACT’s multi-node well package FWL4 produced more realistic drawdowns as evidenced during the 1913–1998 transient simulation (refer to Fig. 11b) as well as during well withdrawals from a hypothetical well (refer to Fig. 12a). Specifically, SURFACT allowed pumping rates to meet full withdrawal demand only in the context of aquifer capabilities. As soon as the modeled water table reached the bottom elevation of a simulated well, then SURFACT automatically adjusted pumping stresses such that heads were fixed to well-bottom conditions. Consequently, the unconfined SURFACT version of the DVRFS did not produce unrealistic drawdowns as evidenced by the confined MODFLOW DVRFSM. While SURFACT provides a more robust and viable model of the DVRFS in terms of drawdown, it does produce higher heads than the MODFLOW model, with the most notable differences occurring in the mountainous regions. Mountain ranges in the DVRFSM are generally defined with large cell thicknesses, relatively large recharge rates and low hydraulic conductivity, which, in combination, can produce significantly larger unconfined heads compared to those modeled with a confined system. Differences are a function of the relative saturation of a modeled cell with reduced interblock conductance and elevated heads in SURFACT driven by the term krw and the ratio of recharge to hydraulic conductivity. The confined DVRFSM is currently used by several agencies (e.g. DOE, USGS) and is considered the most advanced tool for assessing groundwater resources in the NTS area despite admitted conceptual limitations and consequential errors in drawdown estimates (Carroll et al., 2006; Faunt, USGS written communication, 2006). The unconfined SURFACT model represents an evolution of the USGS model toward greater conceptual accuracy and improved stability. The unconfined model uses calibrated values put forth by the USGS and does and equal, if not better, job in model assessment than its confined predecessor. In essence, the unconfined approach helps validate the estimates of these calibrated parameters. The USGS is currently working on a new version of MODFLOW that incorporates a Newton–Raphson non-linear solver and a more robust unconfined solution algorithm. These modifications are similar to the SURFACT algorithm and should increase the capability of MODFLOW to simulate a wide variety of conditions. Acknowledgments Funding for this project was provided by the United States Department of Energy, National Nuclear Security Administration, Nevada Site Office, Hydrologic Resources Management Program. A special thanks to Rich Niswonger and Dave Prudic (retired) at the USGS, Carson City, NV office for their referral to interblock conductance averaging techniques to explain head differences in modeling approaches.
328
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
References Anderman, E.R., Hill, M.C., 2000. MODFLOW-2000, The US Geological Survey Modular Ground-water Model – Documentation of the Hydrogeologic-Unit Flow (HUF) Package: US Geological Survey Open-File Report 00-342, pp. 89. Anderman, E.R., Hill, M.C., 2003. MODFLOW-2000, The US Geological Survey Modular Ground-water Model – Three Additions to the Hydrogeologic-Unit Flow (HUF) Package: Alternative Storage for the Uppermost Active Cells, Flows in Hydrogeologic Units, and the Hydraulic-conductivity Depth-Dependence (KDEP) Capability: US Geological Survey Open-File Report 03-347, pp. 36. Belcher, W.R. (Ed.), 2004. Death Valley Regional Ground-water Flow System, Nevada and California—Hydrogeologic Framework and Transient Ground-water Flow Model: US Geological Survey Scientific Investigations Report 2004-5205, pp. 408.
Region, Nevada and California: US Geological Survey Water-Resources Investigation Report 03-4090, pp. 161. Hill, M.C., 1998. Methods and Guidelines for Effective Model Calibration: US Geological Survey Water-Resources Investigations Report 98-4005, pp. 90. (with application to UCODE, a computer code for universal inverse modeling, and MODFLOWP, a computer code for inverse modeling with MODFLOW). Hill, M.C., Banta, E.R., Harbaugh, A.W., Anderman, E.R., 2000. MODFLOW-2000, The US Geological Survey Modular Ground-water Model – User Guide to the Observation, Sensitivity, and Parameter-estimation Processes and Three Postprocessing Programs: US Geological Survey Open-File Report 00-184, pp. 210. Huyakorn, P.S., Springer, E.P., Guvanasen, V., Wadsorth, T.D., 1986. A three dimensional finite element model for simulating water flow in variably saturated porous media. Water Resources Research 22 (12), 1790–1808. Hydrogeologic, Inc. (HGL), 2007. MODFLOW–SURFACT Software (Version 3.0) Overview: Installation, Registration and Running Procedures, Herndon, VA.. IT Corporation, 1996. Underground Test Area Subproject, Phase I, Data Analysis Task, Vol. VI. Groundwater Flow Model Data Documentation Package: Las Vegas, Nevada. Report ITLV/10972-81 prepared for the US Department of Energy. Laczniak, R.J. Smith, J.L., Elliott, P.E., DeMeo, G.A., Chatigny, M.A., Roemer, G.J., 2001. Ground-water Discharge Determined from Estimates of Evapotranspiration, Death Valley Regional Flow System, Nevada and California: US Geological Survey Water-Resources Investigation Report 01-4195, pp. 51. Laczniak, R.J., DeMeo, G.A., Reiner, S.R., Smith, J.L., Nylund, W.E., 1999. Estimates of Ground-water Discharge as Determined from Measurements of Evapotranspiration, Ash Meadows Area, Nye County, Nevada: US Geological Survey Water-Resources Investigations Report 99-4079, pp. 70. McDonald, M.G., Harbaugh, A.W., 1988. A Modular Three-dimensional Finite Difference Ground-water Flow Model: US Geological Survey Techniques of Water Resources Investigations Book 6. Moreo, M.T., Halford, K.J., La Camera, R.J., Laczniak, R.J., 2003. Estimated groundwater withdrawals from the Death Valley regional flow system, Nevada and California, 1913-98: U.S. Geological Survey Water-Resources Investigations Report 03-4245, 28 p. Reiner, S.R., Laczniak, R.J., DeMeo, G.A., Smith, J.L. Elliott, P.E., Nylund, W.E., Fridrich, C.J., 2002. Ground-water Discharge Determined from Measurements of Evapotranspiration, Other Available Hydrologic Components and Shallow Water-Level Changes, Oasis Valley, Nye County, Nevada: US Geological Survey Water-Resources Investigations Report 99-4239, pp. 65. Romeu, R.H., Noetinger, R., 1995. Calculation of inter nodal transmissivities in finite difference models of flow in heterogeneous porous media. Water Resources Research 31 (4), 943–959.
Contents lists available at ScienceDirect
Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol
An unconfined groundwater model of the Death Valley Regional Flow System and a comparison to its confined predecessor Rosemary W.H. Carroll *, Greg M. Pohll, Ronald L. Hershey Division of Hydrologic Sciences, Desert Research Institute, 2215 Raggio Parkway, Reno, NV 89512, United States
a r t i c l e
i n f o
Article history: Received 12 November 2008 Received in revised form 21 April 2009 Accepted 4 May 2009
This manuscript was handled by P. Baveye, Editor-in-Chief, with the assistance of Hongbin Zhan, Associate Editor Keywords: Death Valley Nevada Test Site MODFLOW SURFACT
s u m m a r y The MODFLOW version of the United States Geological Survey (USGS) Death Valley Regional Flow System (DVRFS) in California and Nevada is conceptually inaccurate in that it models an unconfined aquifer as a confined system and does not accurately simulate unconfined drawdown in transient pumping simulations. The transfer of geologic and hydrologic information from the confined MODFLOW DVRFS model to an unconfined MODFLOW–SURFACT (SURFACT) version was accomplished by maintaining cell structure between models and computing effective cell properties to translate the HUF2 package used in MODFLOW to the BCF4 package used by SURFACT. The confined version of the DVRFS was compared to the unconfined SURFACT version by examining head contour maps and the ability of the SURFACT model to match the 4900 observations of hydraulic head/drawdown, 49 observations of groundwater discharge, and 15 estimates of groundwater fluxes into/out of the model domain. Resultant weighted root mean squared error (xRMSE) for the unconfined SURFACT model was lower than the USGS confined model. Despite a lower xRMSE, unconfined conditions simulated with SURFACT did produce greater heads in mountainous regions compared to the confined MODFLOW with differences most pronounced in regions where cell thickness is large, horizontal conductivity small and recharge large. Difference in computed heads reflects computation schemes employed by both models to estimate interblock conductance. Specifically, interblock conductance for the unconfined SURFACT model is dependent on the relative saturation of a modeled cell while MODFLOW’s confined system is not. Despite head differences, SURFACT simulates comparable flux estimates to MODFLOW (e.g. observed ET, groundwater spring flow, and groundwater flux across model boundaries), while significantly improving transient well drawdown estimates. SURFACT is also capable of producing more realistic estimates of water availability from proposed groundwater development and resultant potential impacts to the region. Ó 2009 Elsevier B.V. All rights reserved.
Introduction There have been numerous water rights applications submitted throughout southern Nevada to offset the needs of growth in Las Vegas. Several applications include groundwater withdrawals adjacent to the Nevada Test Site (NTS) where underground nuclear testing was conducted. If large quantities of groundwater are pumped adjacent to the NTS, the groundwater system could change dramatically. Potential impacts from groundwater pumping include decreasing water levels, reduction in groundwater resources on the NTS, reduction in spring flows adjacent to proposed pumping centers, and the alteration of groundwater flowpaths.
* Corresponding author. Address: Division of Hydrologic Sciences, Desert Research Institute, 2215 Raggio Parkway, Reno, NV 89512 (P.O. Box 4196, Crested Butte, CO 81224), United States. Tel.: +1 970 349 0356; fax: +1 775 673 7363. E-mail addresses: [email protected] (R.W.H. Carroll), [email protected] (G.M. Pohll), [email protected] (R.L. Hershey). 0022-1694/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2009.05.006
In 1998, the US Department of Energy (DOE) Nevada Site Office funded the United States Geological Survey (USGS) to improve upon two previous groundwater flow models of the Death Valley region with the initial intention of understanding groundwater flowpaths and travel times associated with potential movement of radioactive material from the NTS as well as to characterize the groundwater system in the vicinity of Yucca Mountain, and address effects on users down-gradient from the NTS and Yucca Mountain (Belcher, 2004). The first of these earlier models was developed by DOE for the National Nuclear Security Administration/Nevada Site Office (NNSA/NSO) Underground Test Area (UGTA) project (IT Corporation, 1996). The second was developed by the Office of Civilian Radioactive Waste Management’s (OCRWM) Yucca Mountain Project (YMP) and the NNSA/NSO Hydrologic Resources Management Program (HRMP). The resultant USGS Death Valley Regional Groundwater Flow System Model (DVRFSM) has improved upon the two previous models by using newly acquired data and modeling tools (Belcher, 2004). Studies compiled and used in the DVRFSM construction
317
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
include reassessing groundwater discharge via evapotranspiration (ET) (Laczniak et al., 1999, 2001; Reiner et al., 2002; DeMeo et al., 2003) and spring flow (refer to Table C-2, Belcher, 2004), cataloging historical groundwater pumping from 1913 through 1998 (Moreo et al., 2003), reinterpreting groundwater recharge as net infiltration (Hevesi et al., 2002, 2003), assessing model boundary inflows and outflows from regional hydraulic gradients, developing a water budget (Belcher, 2004), and finally incorporating hydraulic conductivity relationships as a function of depth (Belcher et al., 2001, 2002). Carroll et al. (2006) used the DVRFSM to test the impacts of proposed groundwater development south of the NTS and found that, despite increased level of geologic detail and improved water budget accounting, the DVRFSM produced drawdowns on the order of thousands of meters after only a few decades of pumping. Such large drawdowns are not realistic and serve to highlight the most significant limitation of the DVRFSM: the model was built using the confined layer assumption to improve numeric stability. The use of a confined, rather than a ‘‘convertible” (i.e., unconfined), layer type within MODFLOW assumes that the saturated thickness remains constant throughout the entire simulation. Therefore, cells are not allowed to dry or become inactive as the water level decreases below the bottom of a cell. The confined approach was adopted by the USGS because the model was computationally unstable when cells were allowed to convert between confined
and unconfined. Unfortunately the confined approach, while numerically stable, allows unrealistic estimates of drawdown to occur in transient simulations. The primary objective of this study was to produce a numerically stable, unconfined version of the Death Valley Regional Flow System (DVRFS) for future use as a tool to more accurately address potential withdrawals in the NTS region and to determine the potential differences in the solution of the groundwater flow equation given confined and unconfined assumptions.
Numerical models Death Valley Regional Flow System model: MODFLOW The DVRFS and the associated USGS model domain (DVRFSM), with the NTS superimposed, are shown in Fig. 1. The DVRFS is approximately 100,000 km2 in Nevada and California and is bounded by latitudes 35°000 N and 38°150 N and by longitudes 115°000 W and 118°000 W. This system encompasses flow between recharge areas in the mountains of central and southern Nevada and discharge areas of wet playas and springs south and west of the NTS and in Death Valley, California. The flow is strongly influenced by a complex geologic framework and the USGS DVRFSM incorporates the distribution of the flow system’s principal aqui-
Death Valley Region Flow System and Model Boundary
Nevada California
A
A’
0
45 km
Fig. 1. Location of the Death Valley groundwater flow system with the Nevada Test Site (NTS), Death Valley, Spring Mountains, and the Nevada–California state line marked. Inset shows the site with respect to the western United States. Select mountain ranges labeled with white text and select basins labeled in black text. Cross section A–A0 marked.
318
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
NTS
Hypothetical Well
N
0
45 km
Fig. 2. Cell structure in the DVRFSM with the location of a hypothetical well (used for pumping tests) marked.
fers and confining units as well as the principal geologic structures that might affect subsurface flow (Belcher, 2004). The DVRFSM was built with MODFLOW-2000 (McDonald and Harbaugh, 1988; Harbaugh et al., 2000), a public domain, finite-difference code developed by the USGS that is capable of simulating groundwater flow in transient, three-dimensional, anisotropic and heterogeneous systems. MODFLOW’s governing three-dimensional flow equation for a confined aquifer combines Darcy’s Law and the principle of conservation of mass via
@ @h @ @h @ @h S @h K xx þ K yy þ K zz G¼ @x @x @y @y @z @Z b @t
ð1Þ
where x, y, and z are Cartesian coordinates (L), Kxx, Kyy and Kzz are the principal components of saturated hydraulic conductivity along the x, y, and z axes, respectively (L/T), h is hydraulic head (L), G is a volumetric source/sink term (1/T), S is stortivity (dimensionless), b is saturated thickness (L) which is equivalent to cell thickness, and t is time (T). Conversely, Eq. (2) defines groundwater flow given unconfined conditions, where h (L) is the phreatic surface elevation and Sy is specific yield (dimensionless).
@ @h @ @h @ @h @h K xx h þ K yy h þ K zz h G ¼ Sy @x @x @y @y @z @z @t
ð2Þ
Solving for h in Eq. (2) is non-linear problem because h depends on hydraulic conductivity and saturated thickness, with saturated thickness dependent on the water table’s elevation. To address
non-linearity inherent within unconfined aquifer systems, MODFLOW’s updated HUF2 package maintains a confined system, with the top elevation equal to the approximate elevation of the water table, and substitutes Sy for S into Eq. (1). The HUF2 approach to simplifying the non-linearity of an unconfined solution is only applicable when there is little change in saturated thickness. Belcher (2004) describes the DVRFSM in detail, however a consolidated review of the DVRFSM is provided here for clarity. The DVRFSM uses a constant-grid spacing in the horizontal plane of 1500 m by 1500 m. The finite-difference grid is oriented north– south, with 160 columns, 194 rows, and 16 layers for a total of 314,784 active cells (Fig. 2). Cell thicknesses range between 1 m and more than 3000 m, with layer 16 (the bottom layer) reaching 4000 m below sea level. The top elevation of layer one represents the water table (and modeled as a potentiometric surface) and not the actual land surface elevation with cell thickness in layer one increasing substantially in mountain ranges (Fig. 3). Fig. 4 shows the model layer configuration along cross section A–A0 (refer to Fig. 1 for cross section location) in comparison to land surface elevation and simulated hydrogeologic units. Twenty-seven hydrogeologic units were modeled in the DVRFSM via the MODFLOW2000 HUF2 package. A brief description of these units is provided in Table 1. In general, the model’s upper layers portray relatively shallow basin-fill sediments and volcanics, while deeper layers define the regional carbonate aquifer and confining units. Model cells containing more than one hydrogeologic unit are simulated with vertically averaged hydraulic properties.
319
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
Cell Thickness (m) 150 – 500 100 – 150 N
50 – 100
0
1 - 50
45 km
Fig. 3. Cell thickness (m) in layer one.
Model calibration (Belcher, 2004) for steady state conditions prior to year 1912 and for transient conditions for years 1913– 1998 was accomplished by best matching 4899 observations of hydraulic head and drawdown as well as 49 estimates of groundwater discharge and 15 estimates of regional flux into/out of the model domain using MODFLOW-2000 parameter estimation processes (Hill et al., 2000). One hundred parameters were estimated in the calibration process including hydraulic conductivity, vertical anisotropy, drain conductance, fault-associated hydraulic conductivity, recharge zone multipliers, depth decay coefficients for hydraulic conductivity, and storage parameters (S and Sy). Fig. 5 provides an example of where boundary conditions, applied stresses (e.g. wells, drains) and faults (i.e., hydrologic flow barriers) are designated in the DVRFSM’s top layer during the final modeled timestep (year 1998). Fig. 6 shows recharge rates (m/d) applied throughout the modeled domain. Recharge is highest in mountainous regions, in particular the Spring Mountains contribute the greatest volume. Recharge was not modeled in most lowland areas, but those with designated recharge represent irrigation recharge or the re-infiltration of groundwater springs (e.g. Death Valley). Lateral flows into and out of the model boundary were specified as constant head cells that were interpolated from the regional potentiometric surface. Only the region bounding the Spring Mountains was defined as a no flow boundary. Drains were used to define both groundwater ET and spring flows. Discharge was defined by observation and drain conductances were determined by parameter estimation while drain elevations were
set equal to 10 m below the lowest land-surface altitudes for each group of drain cells modeled. These elevations were assumed to adequately represent extinction depths for ET as well as account for springs being located in land-surface depressions (Belcher, 2004). Faults acting as barriers to flow were defined in MODFLOW-2000 with the Horizontal-Flow Barrier (HFB) package. Faults were assumed to extend the entire thickness of the model and the width of each fault was assumed 1 m, while the hydraulic conductivity of each fault was determined by parameter estimation. Nine faults (i.e., horizontal flow barriers) were found to exert some influence on hydraulic heads and are shown in Fig. 5. The DVRFSM simulated hydraulic heads in layer one for the final time step (year 1998) are shown in Fig. 7. SURFACT SURFACT (HGL, 2007) is a fully integrated flow and transport numeric code based on the USGS finite-difference groundwater model MODFLOW. SURFACT potentially improves upon MODFLOW with additional modules. Specifically, SURFACT’s modified BCF4 package performs a complex saturated-unsaturated subsurface flow analysis using a pseudo-soil function for a rigorous treatment of unconfined flow. This is done using robust numerical and matrix solution techniques for greater stability (PCG4) and an adaptive time-stepping scheme (ATO) to improve model efficiency. The three-dimensional movement of water in variably saturated media is given by (Huyakorn et al., 1986),
320
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
(a)
4,000 A Top of Layer 1 = simulated potentiometric surface
3,000 Panamint Mountains Land Surface
Spring Mountains
A'
2,000
Elevation (m)
Greenwater-Black Mountains
1,000
Death Valley
Amargosa Desert
0 -1,000 -2,000 -3,000
Bottom of Layer 16
-4,000
(b)
4,000 0 A 3,000
Elevation (m)
2,000
20
40
60
XCU (K1) LCA (K2) LCA_T1 (K2) VSU_UPPER (K4) MISC Basin fill (K4)
80
100
120
LCCU (K1) LCCU_T1 (K1) VSU_LOWER (K4) LFU (K3)
A'
1,000 0 -1,000 -2,000 -3,000 -4,000 0
30
60
90
120
Distance (km) Fig. 4. Cross section A–A0 (a) configuration of 16 model layers and (b) hydrogeologic units along cross section with units and major rock type (K1, K2, K3 and K4) defined in Table 1.
@ @h @ @h @ @h K xx krw þ K yy krw þ K zz krw G @x @x @y @y @z @z ¼U
@S @h þ Sw SS @t @t
ð3Þ
where additional terms relative to Eq. (1) include krw equal to relative permeability (a function of water saturation, dimensionless), u equal to a drainable porosity (i.e., specific yield Sy), and Sw representing the degree of water saturation (a function of pressure head, dimensionless). When fully saturated conditions exist, then Sw = 1.0, krw = 1.0, and Eq. (3) reduces to Eq. (1), the confined groundwater equation used by MODFLOW. Consequently, the conventional flow equation is used by SURFACT below the water table and in confined systems. While MODFLOW’s other flow packages (BCF and LPF) support several averaging schemes (e.g. harmonic mean, arithmetic mean, arithmetic mean thickness-logarithmic mean hydraulic conductivity) for computing effective interblock transmissivity, MODFLOW’s HUF2 package is limited to only the harmonic mean (i.e., harmonic mean of transmissivity between cell i and cell j is defined as Ti,j = 2TiTj/(Ti + Tj)). The effective transmissivity is equal to the saturated hydraulic conductivity multiplied by the saturated block thickness. For confined systems, the saturated thickness is equal to cell thickness (b). The harmonic mean, by definition, removes the impacts of large outliers, aggregates smaller values and is best used to average
fluxes. However, its estimation automatically underestimates the equivalent interblock conductivity (Romeu and Noetinger, 1995) by biasing toward the lower block value, and in the extreme case never allows a desaturated block (i.e saturated thickness of zero) to resaturated. In contrast, SURFACT’s BCF4 package computes interblock conductance as a product of the weighted harmonic mean of block saturated hydraulic conductivities (Kxx, Kyy, Kzz), relative permeability (krw), and mean flow area. The SURFACT approach to interblock conductance is more stable compared to MODFLOW’s unconfined solution, does not bias toward the lower block value and does not inactivate dry cells during desaturation (HGL, 2007). SURFACT provides a fracture well (FWL) technique for more realistic simulation of groundwater pumping withdrawals compared to the original MODFLOW well (WEL) package. The original WEL package required the user to prescribe the contribution of water from each node (i.e., layer) in a multi-node well and did not automatically adjust pumping rates when water levels dropped below screened intervals. Instead, pumping produced withdrawals into the nonphysical realm. On the other hand, FWL automatically computes (and updates) relative nodal contributions of total well flux rates. For an unconfined system, the FWL package adjusts water withdrawal to feasible volumes from a potentially over-pumped well such that heads do not drop below a well’s screened interval.
321
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328 Table 1 The 27 hydrogeologic units modeled with HUF2 package in the DVRFM and corresponding major rock types. Table modified from Belcher (2004) page 272 Table F-3. Major rock type
Abbreviation
Name
Confining units (K1)
LCCU LCCU_T1 UCCU XCU ICU SCU
Lower clastic-rock confining unit Lower clastic-rock thrust Upper clastic-rock confining unit Crystalline-rock confining unit Intrusive-rock confining unit Sedimentary-rock confining unit
Carbonate rock aquifer (K2)
LCA LCA_T1 UCA
Lower carbonate-rock aquifer Lower carbonate-rock thrust Upper carbonate-rock aquifer
Volcanic rock units (K3)
LFU OVU PVA TMVA CHVU CFTA CFBCU CFPPA YVU WVU BRU
Lava-flow unit Older volcanic-rock unit Paintbrush volcanic-rock aquifer Thirsty canyon-timber mountain volcanic-rock aquifer Calico Hills volcanic-rock unit Crater flat-tram aquifer Crater flat-bullfrog confining unit Crater flat-prow pass aquifer Younger volcanic-rock unit Wahmonie volcanic-rock unit Belted range unit
Basin fill sediments (K4)
OACU YACU YAA LA OAA VSU_LOWER VSU_UPPER
Older alluvial confining unit Younger alluvial confining unit Younger alluvial aquifer Limestone aquifer Older alluvial aquifer Lower volcanic and sedimentary-rock unit Upper volcanic and sedimentary-rock unit
Fig. 5. An example of designated boundary conditions and cell stresses within the DVRFSM (year 1998, top layer).
322
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
Recharge (m/d) 1e-5 – 1e-3 1e-7 – 1e-3 1e-8 – 1e-5
N
0
0 – 1e-8
45 km
Fig. 6. Recharge (m/d) applied to the model’s top layer.
Methods Translation of the DVRFS from MODFLOW to SURFACT Translation of MODFLOW’s geological model packages HUF2 and KDEP (Anderman and Hill, 2000; 2003), which accounts for hydraulic conductivity as a function of depth, into SURFACT’s modified BCF4 package was accomplished by extracting effective cell properties from the original DVRFSM (e.g. horizontal hydraulic conductivity, vertical hydraulic conductivity, and specific storage) using a Fortran code developed for this purpose and compiled as a subroutine of MODFLOW. Specific yield in the DVRFSM was originally calibrated between 0.19 and 0.20 depending on the hydrogeologic unit considered. A value of 0.19 was assumed reasonable and assigned across the entire SURFACT model domain. Fig. 8 shows hydraulic conductivity (K) in the model’s top layer as an example of effective values used by SURFACT. Comparison between MODFLOW and SURFACT The SURFACT DVRFS model was run using USGS-calibrated parameters for the MODFLOW DVRFSM. Comparison between the unconfined and confined modeling approaches included: (1) areal head distributions at the end of the transient model (year 1998), (2) ability to predict observed values of heads, drawdown, groundwater fluxes to springs (i.e., modeled drains) and across modeled boundaries (i.e., constant heads), and (3) estimated drawdowns related to 50 years additional pumping and 50 years recovery begin-
ning in year 1998 from a hypothetical well south of the NTS, a region often considered for groundwater development to support Las Vegas growing need for water. Results and discussion Predicted heads Computed heads at the end of the transient simulation (year 1998) found both MODFLOW and SURFACT consistent in their estimates of flow magnitude and direction. However, differences in predicted heads did occur. In general, the unconfined modeled heads were higher than those predicted by the confined model in mountainous regions with greatest differences occurring in mountain ranges defined with thick model cells (b), low hydraulic conductivity (K) and high recharge flux (w). Specifically, SURFACT predicted heads in excess of 100 m compared to MODFLOW in the Spring Mountains, the Sheep Range, Pintwater Range, Timpahute Range, and Last Chance Range. SURFACT predicted heads on the order of 10–100 m above MODFLOW predicted heads occurred within the Amargosa Range, Black Mountains, Greenwater Range, Belted Range, Black Mountain and the Yucca Mountain region. These regions are defined with thinner cell thicknesses, higher hydraulic conductivities and lower recharge values than other mountain regions. SURFACT and MODFLOW predicted heads were relatively similar for most lowland regions in the model, with SURFACT predicted heads less than MODFLOW in Cactus Flat, Penoyer Valley, Tikaboo Valley, the southern end of the Amargosa Desert
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
323
Fig. 7. DVRFSM-simulated heads (m) for the model’s top layer given 1998 pumping conditions. Contour interval is 200 m.
(i.e., Amargosa Flat) and the northern end of Death Valley. Differences were on the order of 1.0–10.0 m. To understand discrepancies between modeling approaches it is necessary to revisit the groundwater flow equation. The analytical solution for a confined aquifer defined by Eq. (1) and assuming 1dimensional, steady state flow given homogenous, isotropic conditions is,
h ¼ ho
wx2 h1 h0 2L x þ þ Kb 2Kb L
ð4Þ
where h0 and h1 equal heads (L) at x0 = 0 and x1 = L, respectively and w is recharge (L/T). Assuming h0 = h1 = 0 and that the groundwater divide ðq0x ¼ Kbðdh=dxÞ ¼ 0Þ occurs at x = L/2 will produce a maximum head (hmax) equal to,
hmax;c ¼
wL2 8Kb
ð5Þ
Given an unconfined aquifer with identical conditions and using Dupuit assumptions, the 1-dimensional, analytical solution becomes (Fetter, 1994),
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ðh h1 Þx w 2 h ¼ h0 þ 0 þ ðL xÞx L K with hmax at a groundwater divide x = L/2 equal to,
ð6Þ
hmax;uc
sffiffiffiffiffiffiffiffiffi wL2 ¼ 4K
ð7Þ
Fig. 9a compares the analytical difference between hmax,uc and hmax,c for various ratios of w/K and b/L. At low values of w/K and for all values of b/L, hmax,uc > hmax,c differences in solutions are small. Differences become slightly larger with increasing values of b/L. With increased values of w/K and low ratios of b/L, the relationship between unconfined and confined changes such that hmax,uc < hmax,c with differences in heads becoming increasingly significant at larger and larger w/K. In contrast, unconfined solutions for large b/L surpass confined heads given all examined quantities of w/K, with substantial differences occurring at large w/K. Fig. 9b and c compare 1- and 3-dimensional numeric results, respectively, for unconfined heads computed with SURFACT (hmax,SF) and confined heads computed with MODFLOW (hmax,MF). Cell length in all simulations was held constant (L = 1500 m) while cell thickness (b) was allowed to vary. Recharge (w) and vertical hydraulic conductivity (Kzz) were held constant while Kxx, Txx and MODFLOW’s Vcond (Kzz/b) were adjusted accordingly. Layer one bottom elevations were defined at the same elevation as h0 and h1 (i.e., 0 m). This simplified the analysis since hmax,SF = 0 implies a krw = 0, and an hmax,SF P b implies krw = 1.0. Differences between analytic and numeric results reflect computation schemes employed by both models to estimate interblock
324
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
Fig. 8. Effective hydraulic conductivity (K) in the model’s top layer.
conductance and the relationship of predicted hmax,SF to the top elevation of the modeled cell (in this example top elevation of layer 1 = b). MODFLOW’s confined system uses a constant saturated thickness and a constant saturated hydraulic conductivity. SURFACT’s unconfined approach, however, allows for variable saturation of the modeled cell. For the situation where hmax,SF < b, then krw < 1.0 and the effective hydraulic conductivity used by SURFACT is less than the pre-prescribed saturated hydraulic conductivity used by MODFLOW (both in x- and/or z- directions). Thereby, SURFACT produces greater heads than MODFLOW, with greatest differences observed at higher w/K. The difference between hmax,SF and hmax,MF is smaller than its counterpart in the analytical solution (hmax,uchmax,c) because of MODFLOW’s use of transmissivities and subsequent under estimation of interblock conductance. Fig. 9b, b/L = 1.0 illustrates this point. For the situation where hmax,SF P b, then SURFACT uses krw = 1.0 and the flow area between cells is limited to the cell thickness, even if predicted heads lie above the top elevation of the modeled cell. Therefore, the unconfined and confined systems estimate equivalent cell transmissivities. Given these circumstances, unconfined heads predicted by SURFACT begin to converge at increasing values of w/K toward solutions in which MODFLOW heads are slightly larger than SURFACT heads. Despite having numerically equal transmissivities, heads do not converge toward an equal solution because SURFACT computes interblock conductance with
hydraulic conductivity, not transmissivity, and conductance is less biased toward low conductance compared to MODFLOW. Convergence toward nearly equal solutions at large w/K was not witnessed in the analytical solution because the area of flow conveyance for the unconfined system was not limited by b, but allowed to increase to the predicted value hmax,uc. The resultant larger conductance forced the unconfined system to predict significantly smaller heads compared to the confined solution. For the 1-dimensional numeric model, conductance was limited to the horizontal direction (Kxx), while the 3-dimensional numeric model included vertical conductance (Kzz). Resultant leakance from the top layer to deeper layers caused heads in the top layer of the 3-dimensional model to be lower than those in the 1-dimensional model for any given value of w/K. The impact was that SURFACT heads remained less than cell thickness for larger and larger values of w/K. The principles discussed above were applied to the cross section A–A0 (refer to Fig. 1), which transects several different mountain ranges and low lying basins. Recharge, hydraulic conductivity, land surface elevation, cell top and bottom elevations as well as modeled heads were extracted from layer one on a cell-by-cell basis. Quantities were plotted in Fig. 10. Included is a plot of relative cell saturation for unconfined conditions, (hSF cell bottom elevation)/b, (refer to Fig. 10d). A constant head boundary condition forced MODFLOW and SURFACT to predict the same heads at location A,
325
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328 600
400
Elevation (m)
500
b/L = 1
300
h max,uc - h max,c
200 100 0
b/L = 0.1
-100
4,000 A 3,500 3,000 2,500 Panamint Mountains Greenwater-Black 2,000 Ranges Land Surface 1,500 Amargosa Desert 1,000 Death 500 Valley 0 MODFLOW SURFACT -500 020406080100120140 0.20
Distance (km)
A
-200
Spring Mountains
A'
A'
0.15
-300
b/L
b/L = 0.01
-400 -500 -600 0.0001
0.10 0.05
0.001
0.01
0.1
0.00
1
w/K
020406080100120140
1.E+00
w/K
400
b/L = 1
300
A'
Distance (km)
A'
1.E-04 1.E-05 1.E-06
200 100
1.E-03
0.0 1.E-07
b/L = 0.1
020406080100120140
b/L = 0.01
-100 -200 -300 -400 -500
0.001
0.01
0.1
6.00 4.00 2.00
-2.00
h SF-h MF (m)
400
75
b/L = 1
50
0.0
0.00
500 020406080100120140
1
w/K 100
A
Krw = 0.0
Relative Saturation
0
Krw = 1.0
h max,SF - h max,MF
Distance (km)
1.E-02
500
-600 0.0001
A
1.E-01
600
A
A'
Distance (km)
300 200 100 0
h max,SF - h max,MF
-100
25
0
b/L = 0.1
20
40
60
80
100
120
140
Distance (km)
0
b/L = 0.01
Fig. 10. Cross section A–A0 layer one cell values at the end of the transient simulation (year 1998) for: (a) land surface elevation, MODFLOW predicted heads (hMF) and SURFACT (hSF) predicted heads. The shaded area depicts distance between top and bottom elevation of the cell, (b) b/L, (c) w/L, (d) relative saturation and (e) hSFhMF.
-25 -50 -75 -100 0.0001
0.001
0.01
0.1
1
w/K Fig. 9. Difference in steady state, maximum head versus w/K for various values of b/ L: (a) 1-dimensional analytical solution for unconfined (hmax,uc) and confined (hmax,c) systems, (b) 1-dimensional numeric solution for unconfined SURFACT (hmax,SF) and confined MODFLOW (hmax,MF), (c) 3-dimensional numeric solutions for unconfined SURFACT and confined MODFLOW.
while a no-flow boundary allowed heads to vary at location A0 . Low values of w/K kept MODFLOW and SURFACT heads nearly equivalent throughout the Penamint range. SURFACT heads exceeded MODFLOW heads only in regions where SURFACT’s relative saturation was less than 1.0, thereby reducing krw and the equivalent SURFACT conductivity in relationship to MODFLOW. More significant deviation occurred along the western edge of the Greenwater and Black Ranges. In this region, SURFACT produced a krw < 1.0 and
Table 2 A comparison of MODFLOW and SURFACT weighted root mean squared errors (xRMSE) for each type of observation.
xRMSE Observation
Number
MODFLOW
SURFACT
Heads Drawdown Drains (i.e., ET) Constant heads Total
1930 2969 49 15 4963
100.19 18.95 3.93 5.40 64.17
100.82 14.06 10.36 4.94 63.89
values of w/K were large enough to produce SURFACT heads 10– 50 m above MODFLOW heads. The greatest difference occurred at a distance of 30.0 km, when relative saturation dropped below zero (heads below top layer), pegging krw to a value of 0.0. In the Spring Mountains, the ratio of w/K was much higher than in the Green-
326
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
water-Black Range forcing SURFACT heads to exceed the top elevation of layer one despite substantial values of b/L. None the less, SURFACT heads were in excess of 400 m compared to MODFLOW heads. Given krw = 1.0, it is presumed that SURFACT heads would have been larger if the top elevation of the cell was raised to land surface elevations and boundary conditions maintained. For the unconfined system, placing layer one’s top elevation at land surface would be preferable, however the USGS calibrated depth decay coefficients in the HUF2 package were based on layer one’s current top elevation. Redefining the top elevation to land surface would negate this calibration. SURFACT and MODFLOW predicted heads were nearly equivalent in low lying areas along A–A0 since w/K was generally low and krw was approximately 1.0.
where Xo is the observed value, Xp is the predicted value, n is the number of observations, and x is the inverse variance. Predicted values derived from multiple layers were computed using fractional inputs from each contributing layer as defined by Belcher (2004) and discussed by Hill (1998) and Hill et al. (2000). A comparison of xRMSE values is provided in Table 2 and also shown in Fig. 11a for heads, Fig. 11b for drawdown, Fig. 11c for drain fluxes and Fig. 11d for constant head fluxes. Both models perform similarly with the SURFACT model having a slightly smaller total xRMSE. In addition the xRMSE for drawdown was significantly smaller than the MODFLOW model. The larger xRMSE for drawdown in the MODFLOW model is the result of extremely large drawdown values that were simulated by MODFLOW using the confined governing equation.
Observations Well withdrawal Model results were compared to observed data using the weighted root mean squared error (xRMSE),
Pn i
xRMSE ¼
(a)
xi ðX 0i X pi Þ2
!0:5 ð8Þ
n 5,500
(b)
MODFLOW
300
SURFACT
250
200
Modeled Drawdown (m)
4,500
Modeled Heads (m)
Last, a hypothetical well was placed in the model’s top layer to the south of the NTS (row 124, column 99, refer to Fig. 1). This well’s location was chosen for its proximity to the NTS and because predicted heads by the MODFLOW and SURFACT models were nearly
3,500
2,500
1,500
150 100 50
0
500 MODFLOW
-50
SURFACT
-500 -500
1,500
500
2,500
3,500
-100 -100
4,500
Observed Heads (m)
(c)
-50
0
50
100
150
200
250
300
Observed Drawdown (m)
(d)
0
90,000
MODFLOW
MODFLOW
-5,000
SURFACT
-10,000
3
-15,000 -20,000 -25,000 -30,000 -35,000
40,000
-10,000
-60,000
90,000
40,000
-10,000
0
-5,000
-10,000
-15,000
-20,000
-25,000
-30,000
-35,000
-40,000
Observed Drains(m 3/d)
-60,000
-110,000
-40,000
-110,000
Modeled Drains (m /d)
3
Modeled Constant Head (m /d)
SURFACT
Observed Constant Head (m 3/d)
Fig. 11. Comparison of MODFLOW- and SURFACT-predicted values to observed values used in the calibration of the MODFLOW model: (a) heads, (b) drawdowns, (c) drain fluxes, and (d) constant head fluxes.
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
not realistically represent the actual unconfined system, which cannot support drawdown of heads below the well unless an adjacent and deeper well is pumped simultaneously.
MODFLOW
785
SURFACT
Head in Well (m)
760 735
Conclusions
Bottom Elevation of Well
710 685 660 Pumping Ends
635
(a)
610 585 0
10
20
30
40
50
60
70
80
90
100
70
80
90
100
Years of Pumping 6,000
5,000
Well Flux (m 3/d)
327
4,000
Prescribed
3,000
2,000 Actual
1,000
(b)
0 0
10
20
30
40
50
60
Years of Pumping Fig. 12. (a) Comparison of head drawdown in a hypothetical well south of the NTS between the MODFLOW (confined) and SURFACT (unconfined) models given 50 years pumping at 5000 m3/d (1480 acre ft/yr) followed by 50 years of recovery. (b) Prescribed pumping rate compared with the SURFACT model actual withdrawal.
identical at the end of transient simulation (year 1998), which simplifies the analysis of subsequent drawdown with respect to the bottom of the well. The bottom of the well was established at the bottom elevation of the cell containing the well (729 m above mean sea level). A pumping rate of 5000 m3/d (1480 acre ft/yr) was initiated at the end of the transient model and then allowed to recover for an additional 50 years with no pumping. Results (Fig. 12) show that the SURFACT model can only meet prescribed pumping rates during the first 10 years of pumping, after which drawdown approached the base of the well and pumping rates were reduced to feasible quantities that the aquifer could supply to maintain the water table at the base of the well. Maximum modeled drawdown in the well reached only 60.6 m, or the bottom elevation of the well. After pumping ceased, slow rebound in the water table resumed such that, after 50 years, water in the well rebounded approximately 5.5 m to an elevation of 734.5 m. Unlike the SURFACT model, MODFLOWs confined conditions allowed full withdrawal demand to be met for the entire simulation. Divergence in SURFACT and MODFLOW simulated drawdown began after 10 years when the SURFACT model reduced the amount of stress on the well to keep heads just above the bottom elevation of the well while the MODFLOW model maintained the prescribed stress. At approximately 15 years of pumping, the MODFLOW model’s predicted heads dropped below the bottom of the well, and after 50 years of pumping, predicted heads fell to 586.4 m, or more than 142 m below the bottom of the well. While large drawdowns may be numerically possible in a confined aquifer, they do
Representation of the DVRFS using SURFACT is an improvement over its MODFLOW predecessor for transient simulations. The original USGS model assumed that the entire model domain was in a confined state to ensure numerical stability. The uppermost layer used a specific yield storage parameter to represent unconfined conditions, but the cell thickness does not change in response to water level changes. The result was projected drawdown estimates from well withdrawals were unrealistic. On the other hand, SURFACT’s saturated-unsaturated subsurface flow analysis using a pseudo-soil function produced a numerically stable model of the DVRFS using unconfined conditions. In addition, SURFACT’s multi-node well package FWL4 produced more realistic drawdowns as evidenced during the 1913–1998 transient simulation (refer to Fig. 11b) as well as during well withdrawals from a hypothetical well (refer to Fig. 12a). Specifically, SURFACT allowed pumping rates to meet full withdrawal demand only in the context of aquifer capabilities. As soon as the modeled water table reached the bottom elevation of a simulated well, then SURFACT automatically adjusted pumping stresses such that heads were fixed to well-bottom conditions. Consequently, the unconfined SURFACT version of the DVRFS did not produce unrealistic drawdowns as evidenced by the confined MODFLOW DVRFSM. While SURFACT provides a more robust and viable model of the DVRFS in terms of drawdown, it does produce higher heads than the MODFLOW model, with the most notable differences occurring in the mountainous regions. Mountain ranges in the DVRFSM are generally defined with large cell thicknesses, relatively large recharge rates and low hydraulic conductivity, which, in combination, can produce significantly larger unconfined heads compared to those modeled with a confined system. Differences are a function of the relative saturation of a modeled cell with reduced interblock conductance and elevated heads in SURFACT driven by the term krw and the ratio of recharge to hydraulic conductivity. The confined DVRFSM is currently used by several agencies (e.g. DOE, USGS) and is considered the most advanced tool for assessing groundwater resources in the NTS area despite admitted conceptual limitations and consequential errors in drawdown estimates (Carroll et al., 2006; Faunt, USGS written communication, 2006). The unconfined SURFACT model represents an evolution of the USGS model toward greater conceptual accuracy and improved stability. The unconfined model uses calibrated values put forth by the USGS and does and equal, if not better, job in model assessment than its confined predecessor. In essence, the unconfined approach helps validate the estimates of these calibrated parameters. The USGS is currently working on a new version of MODFLOW that incorporates a Newton–Raphson non-linear solver and a more robust unconfined solution algorithm. These modifications are similar to the SURFACT algorithm and should increase the capability of MODFLOW to simulate a wide variety of conditions. Acknowledgments Funding for this project was provided by the United States Department of Energy, National Nuclear Security Administration, Nevada Site Office, Hydrologic Resources Management Program. A special thanks to Rich Niswonger and Dave Prudic (retired) at the USGS, Carson City, NV office for their referral to interblock conductance averaging techniques to explain head differences in modeling approaches.
328
R.W.H. Carroll et al. / Journal of Hydrology 373 (2009) 316–328
References Anderman, E.R., Hill, M.C., 2000. MODFLOW-2000, The US Geological Survey Modular Ground-water Model – Documentation of the Hydrogeologic-Unit Flow (HUF) Package: US Geological Survey Open-File Report 00-342, pp. 89. Anderman, E.R., Hill, M.C., 2003. MODFLOW-2000, The US Geological Survey Modular Ground-water Model – Three Additions to the Hydrogeologic-Unit Flow (HUF) Package: Alternative Storage for the Uppermost Active Cells, Flows in Hydrogeologic Units, and the Hydraulic-conductivity Depth-Dependence (KDEP) Capability: US Geological Survey Open-File Report 03-347, pp. 36. Belcher, W.R. (Ed.), 2004. Death Valley Regional Ground-water Flow System, Nevada and California—Hydrogeologic Framework and Transient Ground-water Flow Model: US Geological Survey Scientific Investigations Report 2004-5205, pp. 408.
Region, Nevada and California: US Geological Survey Water-Resources Investigation Report 03-4090, pp. 161. Hill, M.C., 1998. Methods and Guidelines for Effective Model Calibration: US Geological Survey Water-Resources Investigations Report 98-4005, pp. 90. (with application to UCODE, a computer code for universal inverse modeling, and MODFLOWP, a computer code for inverse modeling with MODFLOW). Hill, M.C., Banta, E.R., Harbaugh, A.W., Anderman, E.R., 2000. MODFLOW-2000, The US Geological Survey Modular Ground-water Model – User Guide to the Observation, Sensitivity, and Parameter-estimation Processes and Three Postprocessing Programs: US Geological Survey Open-File Report 00-184, pp. 210. Huyakorn, P.S., Springer, E.P., Guvanasen, V., Wadsorth, T.D., 1986. A three dimensional finite element model for simulating water flow in variably saturated porous media. Water Resources Research 22 (12), 1790–1808. Hydrogeologic, Inc. (HGL), 2007. MODFLOW–SURFACT Software (Version 3.0) Overview: Installation, Registration and Running Procedures, Herndon, VA.