Assessing Potential Effects of Changes in Water Use in the Middle Carson River Basin with a Numerical Groundwater-Flow Model, Eagle, Dayton, and Churchill Valleys, West-Central Nevada
- Document: Report (18 MB pdf) , HTML , XML
- Data Releases:
- Data for the report Geologic Framework and Hydrogeology of the middle Carson River basin, Eagle, Dayton, and Churchill Valleys, West-Central Nevada
- Data for the report assessing potential effects of changes in water use in the middle Carson River Basin with a numerical groundwater-flow model, Eagle, Dayton, and Churchill Valleys, west-central Nevada
- MODFLOW-NWT model used to simulate potential effects of changes in water use in the middle Carson River Basin, Eagle, Dayton, and Churchill Valleys, west-central, Nevada
- Download citation as: RIS | Dublin Core
The author would like to thank the Bureau of Reclamation for their long-running support of not only the modeling effort of the middle Carson River Basin, but also for their continued partnership in addressing water-related issues throughout the entire Carson River Basin. We are especially grateful to Kenneth Parr, Terri Edwards, Tom Scott, Mandy Brinnand, and Caryn Huntt DeCarlo of the Bureau of Reclamation for their patience during the model development and calibration process. We would also like to thank Dave Wathen with the Federal Water Master’s office in Reno, Nevada, for his prompt responses to our inquiries, sharing of available data, and patience as we tried to understand the administration of the Alpine Decree. The following individuals are thanked for providing data, background information, or lending general support: Ed James with the Carson Water Subconservancy District; Matt Dillon with the Nevada State Engineer’s office, Mike Workman with Lyon County Utilities, Rit Palmer and Randall Gray with Carson City Public Works, and Chris Mahannah with Mahannah & Associates.
Doug Maurer (retired); Toby Welborn, Rose Medina, and Joe Gardner (retired), all with the U.S. Geological Survey, are acknowledged for their help and support keeping the report moving along. Reviews by Phil Gardner, Dan Bright, Jill Frankforter, Jason Bellino, and Andy Leaf are gratefully acknowledged.
During the economic boom of the mid part of the first decade of the 2000s in northwestern Nevada, municipal and housing growth increased use of the water resources of this semi-arid region. In 2008, when the economy slowed, new housing development stopped, and immediate pressure on groundwater resources abated. The U.S. Geological Survey, in cooperation with the Bureau of Reclamation, began a hydrogeologic study of the middle Carson River Basin. The first half of the study reviewed and synthesized previous geologic studies and contributed new datasets that served as a foundation for a three-dimensional, transient numerical model of groundwater and surface-water flow for the middle Carson River Basin extending from Eagle Valley to Churchill Valley. The model can be used to evaluate the effects of proposed alternative management strategies on groundwater sustainability, flows in the Carson River, and routine operation of Lahontan Reservoir and can also provide a basis for basin-wide investigations seeking to quantitatively evaluate the effects of climate change or yet-to-be-determined alternative management strategies.
The middle Carson model was constructed using the U.S. Geological Survey groundwater modeling software MODFLOW-NWT. MODFLOW is widely used groundwater modeling software and is well-suited for evaluating groundwater and surface-water interactions. The model uses 550-feet square grid cells that align with the previously published model for Carson Valley (adjacent upstream valley). Six grid layers with more finely resolved vertical resolution near the perimeter of the active model domain and near surface-water features, compared to other areas of the active model domain, hone the simulated groundwater and surface-water exchanges. In addition to simulating groundwater and surface-water interaction, crop and phreatophyte evapotranspiration, lake evaporation, mountain-front recharge, recharge from irrigation return flows, and groundwater pumping are also simulated. Surface-water flow entering the model domain, including the Carson River, tributary inflow from perennial streams in Eagle Valley, and trans-basin imports through the Truckee Canal (surface water diverted from the Truckee River) are specified according to U.S. Geological Survey streamgage records. Groundwater pumpage and surface-water diversions to 10 agricultural ditches and the managed release from Lahontan Reservoir, at the end of the middle Carson River Basin, are specified according to water-manager records.
The model simulation period extended from 2000 through 2010 (January 1, 2000, to December 31, 2010) using 574 weekly stress periods, with a single steady-state stress period at the beginning of the simulation that establishes initial conditions by approximating average conditions during the transient simulation period. All available observations for this period were used during the model calibration process, performed using automated parameter-estimation software. Calibration targets included observations of groundwater elevations in wells, streamflow, differences in observed streamflow between successive streamgages and actual evapotranspiration from irrigated lands. Among all 5,296 simulated and observed groundwater level pairs, the mean error was 1.42 feet; the mean absolute error, 7.71 feet; and the percent bias was −0.1 percent.
Three alternative management scenarios, run using the entire period of analysis (2000–10), were simulated to improve understanding of the potential effects of (1) loss of irrigated agricultural lands following conversion of water-rights to municipal groundwater rights; (2) reclaiming treated wastewater with induction wells; and (3) exercising permitted but under-utilized groundwater rights. Scenarios 2 and 3 were further explored using two and four subscenarios, respectively. Simulated scenario results ranged from having little effect on the groundwater system relative to a baseline simulation to having spatially extensive and large groundwater-level declines (10 to 20 feet) compared to the baseline simulation. None of the simulated scenarios increased delivery of river flows to Lahontan Reservoir. On the contrary, one of the subscenarios under alternative management scenario 3 led to surface-water delivery shortfalls of more than 10,000 acre-feet per year.
Future model improvements may include an extension of the model simulation period backward and forward in time and directly linking it to the upstream Carson Valley groundwater model. Furthermore, converting this MODFLOW model to a GSFLOW model, which fully integrates groundwater and surface-water flows including precipitation runoff and infiltration, may provide an improved tool for comprehensive management of water-resources in the middle Carson River Basin.
The Carson River as well as the alluvial aquifers over which it flows are vital water resources for the communities in Eagle, Dayton, and Churchill Valleys (figs. 1, 2). Population growth in all three valleys (fig. 3) increased pressure on water resources. Moreover, population growth is projected to continue; for example, in Carson City, the Nevada State Demographer’s Office predicted that the 2012 population of about 55,400 would increase to about 64,000 by 2032, representing a 15-percent increase in total population (Hardcastle, 2013). Neighboring counties that also rely on Carson River water, at least in part, are projected to grow as well. Douglas, Lyon, and Churchill Counties are forecasted to grow by 5, 16, and 26 percent, respectively (Hardcastle, 2013). As competition for water continues to mount in the middle Carson River (MCR) basin and surrounding areas, groundwater pumping is likely to continue to increase (Maurer and others, 2008, figs. 9 and 10). In addition to spurring additional pumping, competition for water could lead to transfers of surface-water rights traditionally used for agriculture, thereby reducing an important source of recharge (sometimes referred to as irrigation return flows) to the alluvial aquifers.
Table 1.Federal Water Master streamgages on diversions from the middle Carson River.
[USGS, U.S. Geological Survey; NV, Nevada; USFWM, U.S. Federal Water Master]
The “Exploration and Graphics for River Trends” (EGRET) tool (Hirsch and De Cicco, 2015) was applied to evaluate long-term streamflow at the U.S. Geological Survey (USGS) Carson River near Carson City, NV streamgage (10311000 [U.S. Geological Survey, 2018], fig. 2A and table 1). Results of the locally weighted scatterplot smoothing method (LOWESS) in the EGRET analysis (fig. 4) showed a clear downward trend in the annual mean flow, indicative of a diminished river-water supply entering the middle Carson study area, which is consistent with the interpretation of Maurer and others (2008). From 1941 to 1970, the mean annual Carson River flow decreased by 14 percent before experiencing a 4.5-percent increase from 1970 to 1975. Between 1975 and 2011, the mean annual Carson River flow decreased by an additional 20 percent. From 1941 to 2011, the decline in the mean daily inflow rate to the MCR was roughly 0.4 percent per year, equating to an overall decrease of 28 percent (fig. 4). Low-flow conditions became more common during May and June, which previously had typically been the peak runoff period. For example, the average annual decline in the 7-day minimum flow was about 1 percent at the Carson River near Carson City, Nevada streamgage between 1941 and 2011 (fig. 5), supporting the conclusion of long-term diminished flow entering the MCR. Decreases in the 7-day minimum flow rate at the Carson River near Carson City, Nevada streamgage between 1941 and 1970 and between 1975 and 2011 were 38 and 34 percent, respectively. Between 1970 and 1975, the 7-day minimum flow rate was generally level.
The water supply in the Carson River basin is dominated by melting snow and will likely experience earlier runoff in a warming climate even without any changes in precipitation intensity (Barnett and others, 2005). Moreover, local temperature records indicate a warming trend during the past century that could lead to increased water demand for irrigation (Pederson and others, 2010). For example, the first day of the year on which temperatures exceed 90 degrees Fahrenheit (°F) occurs roughly 10 days earlier after 2000 than it did before 1910 (fig. 6). In addition, the last day on which temperatures exceed 90 °F has shown a trend toward occurring later in the year. Moreover, the first and last frost days are later and earlier in the year, respectively, compared to a century ago (fig. 7). Together, these trends indicate longer growing seasons and an associated increase in crop-water demands.
As surface-water diversions are changed, concomitant changes to groundwater-flow dynamics related to altered crop consumptive use, irrigation returns (that is, recharge from flood irrigation), and, subsequently, groundwater discharge to the Carson River are expected, such that historic in-stream flow rates also could be altered. Depending on how water needs are addressed, increased groundwater pumping from existing wells or pumping from new wells could potentially decrease surface-water flows in the MCR. This, in turn, could affect downstream water users by diminishing deliveries to Lahontan Reservoir and threaten the sustainability of groundwater supplies in the MCR. Development of a calibrated groundwater-flow model, using a computer code like MODFLOW (Harbaugh, 2005) to represent the MCR system, may assist water-resource managers by providing an appraisal of current conditions as well as identification of potential effects of alternative management strategies before implementation.
The USGS, in cooperation with the Bureau of Reclamation (Reclamation), began this study in 2008 to develop a numerical groundwater model to simulate groundwater and surface-water interactions in the Carson River basin upstream of Lahontan Dam and downstream of Carson Valley. The first part of the study collected and compiled hydrogeologic information that was published in Maurer (2011). Additionally, a companion study (Jeton and Maurer, 2011) documents the development of precipitation-runoff models of selected watersheds in the middle Carson River basin to estimate runoff and groundwater recharge. Both documents serve as a basis for the numerical groundwater flow model documented herein.
Purpose and Scope
This report presents the development of a numerical groundwater-flow model calibrated to historical observations related to the ground- and surface-water systems along the MCR. The model is designed to assess possible changes to available water supplies in the MCR, largely driven by residential and light industrial development and associated water use. Changes to water availability are assessed by simulating the effects of alternative water-management scenarios that vary the timing and distribution of groundwater pumping and surface-water deliveries.
The numerical model was developed using MODFLOW-NWT (Niswonger and others, 2011) and parameter-estimation (PEST) software utilities (Doherty, 2010a, b, c; Doherty and Hunt, 2010). The model was designed, through appropriate selection of packages available with MODFLOW-NWT (Niswonger and others, 2011), to help water managers assess the future effects of alternative management decisions made today in response to a growing need for water. The model synthesizes a broad body of information about the MCR that was collected as part of a cooperative study between the USGS and Reclamation that began in 2008, the results of which are reported in Jeton and Maurer (2011) and Maurer (2011). Data from the cooperative studies were compiled and used to develop (1) the hydrogeologic framework; (2) distributions of precipitation, runoff, and infiltration in 12 perennial and ephemeral tributaries along the MCR; (3) key characteristics of the surface-water network, including identification of gaining and losing stretches of the middle Carson River study area; and (4) surface-water flow and groundwater hydraulic head observations (Jeton and Maurer, 2011; Maurer, 2011). Thus, the MCR modeling study codifies results from multiple studies between the USGS and Reclamation to help guide resource managers with the long-term sustainable management of water resources in the region.
Situated on the western edge of the Great Basin, the headwaters of the Carson River are in Alpine County, California, at elevations above 10,000 feet (ft). The East and West Forks of the Carson River flow out of the mountain block near the USGS streamgages of East Fork Carson River near Gardnerville, Nev. (10309000), and West Fork Carson River at Woodfords, Calif. (10310000; fig. 1). The confluence of the East and West Forks of the Carson River is just east of Genoa, Nev., on the west side of Carson Valley. From this point, the mainstem Carson River flows north out of Carson Valley and into the MCR study area at the USGS Carson River near Carson City, Nev. streamgage (10311000; fig. 1). From this point, the Carson River flows north along the westernmost edge of the Dayton Valley hydrographic area (HA 103; figs. 1, 2A). The Carson River never enters the Eagle Valley hydrographic area. It traverses the eastern edge of Eagle Valley but remains in the Riverview subbasin of the Dayton Valley hydrographic area (fig. 1).
Maurer (2011) describes four subbasins in the Dayton Valley hydrographic area: (1) Riverview, (2) Moundhouse, (3) Carson Plains, and (4) Bull Canyon. Maurer (2011, p. 4) is careful to point out there is no single valley named “Dayton Valley,” the term is used locally to refer to the Carson Plains subbasin because the town of Dayton—the largest town in the topographic valley—is there. This report refers to the entire HA 103 area as “Dayton Valley,” unless a specific subbasin is mentioned (fig. 1).
Upon entering the study area at the Carson River near the Carson City streamgage (10311000, fig. 1), the Carson River flows north along the eastern side of Eagle Valley in the Riverview subbasin (of the Dayton Valley hydrographic area). The river then flows east into the Moundhouse subbasin of the Dayton Valley hydrographic area. Mexican Ditch is the only diversion from the Carson River to Eagle Valley. The Carson River at Deer Run Road near Carson City streamgage (10311400) is at the divide between the Riverview and Moundhouse subbasins (fig. 2A).
From the Carson River at Deer Run Road near Carson City streamgage, the Carson River flows through a relatively narrow canyon in the Moundhouse subbasin upstream from the town of Dayton. The river exits the canyon near the Moundhouse–Carson Plains subbasin divide. After exiting the canyon, six diversion dams divert flow to irrigation delivery ditches on both sides of the river in the Carson Plains subbasin (table 2; fig. 2A). After traversing Carson Plains, the Carson River enters a second canyon at the Carson Plains–Bull Canyon subbasin divide (fig. 2A), where two more ditches, the Houghman-Howard and Buckland ditches, divert water from the river for irrigation (table 2; fig. 2B). The river continues flowing eastward past Table Mountain and Churchill Butte, entering the Churchill Valley hydrographic area near the USGS Carson River near Fort Churchill streamgage (10312000; fig. 2B). From this point, the river flows approximately 8 miles (mi) east to Lahontan Reservoir (fig. 2B).
Table 2.Active and retired diversions in the middle Carson River model area, including estimates of the irrigated, fallow, retired, and urbanized land areas (in acres). Ditches are listed in order from upstream to downstream.
[Ditches are listed in order from upstream to downstream. Abbreviations: no., number; iseg, streamflow routing package segment number]
The Newlands Project, the first reclamation project built in the United States because of the Reclamation Act of 1902, was completed in 1915 and led to the construction of Lahontan Reservoir (Wilds, 2010). The need for the reservoir arose from the considerable variability in annual streamflow in the Carson River, making it an unreliable resource during dry years. Lahontan Reservoir receives streamflow from the Carson River and the Truckee River through the Truckee Canal (figs. 1, 2B). Releases from the reservoir provide water to approximately 56,000 irrigated acres through a complex network of canals (Maurer and others, 2008, p. 4). The reservoir reaches capacity at 295,591 acre-feet (acre-ft) at a corresponding stage of 4,162 ft, of which 91 acre-ft is dead storage below the reservoir outlet elevation of 4,070 ft. At full capacity, the surface area of Lahontan reservoir is 13,470 acres (Stockton and others, 2003). Mean annual flows into Lahontan reservoir during the period of analysis as measured by the Carson River near Fort Churchill streamgage (10312000), the last streamgage before flow enters the reservoir, ranged between 108,700 acre-feet per year (acre-ft/yr), or 150 cubic feet per second (ft3/s), in 2007 (15th smallest year on record out of 81 years) and 550,600 acre-ft/yr (760 ft3/s) in 2006 (11th largest year on record) during the study period. In 1977, the driest year on record at the Fort Churchill streamgage since streamflow monitoring began in 1912, the total gaged volume was 26,300 acre-ft. Only 6 years later, in 1983, more than 800,000 acre-ft was measured at the Fort Churchill streamgage because of one of the largest snowpacks on record.
Elevation exerts a strong effect on the spatial distribution of precipitation in the Carson River Basin (Houghton and others, 1975). In addition to elevation, Maurer and Halford (2004) note the rain shadow in the Carson River Basin created by the Sierra Nevada, indicating longitude also is associated with the spatial distribution of precipitation in the Carson River Basin. Annual precipitation at the top of the Carson Range is approximately 30 inches per year (in/yr) and drops to roughly 15 in/yr near the top of the Pine Nut Mountains, the next mountain range to the east (National Oceanic and Atmospheric Administration, 2002). The average annual precipitation in Eagle Valley, which separates the Northern Carson Range from the Northern Pine Nut Mountains, is about 10 in/yr (National Oceanic and Atmospheric Administration, 2002). In Churchill Valley, near the eastern edge of the study area, precipitation is closer to 5 in/yr (National Oceanic and Atmospheric Administration, 2002).
Among the tributaries to the Carson River in the middle Carson River study area, only three, Clear Creek, Ash Canyon, and Kings Canyon, sustain perennial streamflow (fig. 2A). The headwaters of all three are in the higher elevations of the Carson Range, where winter snowpack persists. Tributary streamflow from 10 ephemeral tributaries that contribute flow during spring runoff in excessively wet years or during extreme precipitation are east of the Riverview subbasin and include (from upstream to downstream) Brunswick, Eureka, Hackett, Daney, Eldorado, Gold, Sixmile, Bull-Mineral, Churchill, and Ramsey Canyons (figs. 2A–B).
Owing to the aridity of the region and considerable variability in the annual surface-water supplies, pumped groundwater is often used to meet domestic, municipal, agricultural, industrial, and commercial needs. The total number of domestic wells in Eagle Valley decreased during the simulation period (2000–10), whereas in Dayton and Churchill Valleys, the total number of domestic wells increased during the simulation period (Maurer and others, 2009). Municipal, irrigation, and industrial and commercial wells remained relatively unchanged in Eagle Valley during the first decade of the 2000s. Small increases in the number of industrial and commercial wells in Dayton Valley were likely the result of water-right transfers from irrigation wells. In Churchill Valley, the number of irrigation wells remained constant, whereas the total number of municipal, industrial, and commercial wells increased during 2000 to 2010 (Maurer and others, 2009).
The history, structure, and geometry of the geologic units in the MCR have been described in several previous investigations. Maurer (2011) provides a summary of these investigations and documents the geologic framework that was adopted for the numerical model described in this report. A summary of bedrock units in the MCR and discussion about their control on groundwater flow along the boundaries of the study area is provided in this report (Maurer, 2011).
Tributary watersheds on the west side of Eagle Valley, namely Clear Creek and Ash Canyon, consist of deeply weathered granite (weathered to depths of 100 ft) and contribute small amounts of recharge to the basin-fill sediments. The northern Pine Nut Mountains form the east boundary of Eagle Valley and southern boundary of Dayton Valley and consist of several blocks bounded by north-trending normal faults. The faults have exposed granitic and metamorphic basement rocks on the east side of the Pine Nut Mountains tilting the mountain range to the west (Moore and Archbold, 1969, p. 18). Although separated by the Carson River, the Pine Nut Mountains and Virginia Range are considered a single, complex structural unit. The Carson River crosses a structural and topographic low point between these ranges (Moore and Archbold, 1969, p. 19). Schaefer and Whitney (1992; fig. 2A) used gravity and aeromagnetic data to estimate the maximum thickness of basin-fill deposits at this low point region to be 2,900 ft near the northwestern edge of the valley floor. Using a regression relationship between gravity measurements and known depths to bedrock, Maurer (2011, fig. 9A) used an array of gravity measurements taken from within Dayton Valley (Carson Plains subbasin) to estimate a shallower maximum thickness of 1,000 ft under the Carson River near the center of the valley floor. A unit of consolidated volcanic rock and semi-consolidated sediments of Tertiary age lies under the basin-fill sediment (Maurer, 2011), followed by a unit consisting of consolidated pre-Cenozoic rock that forms the basement of the model.
The Virginia Range separates the Truckee River (to the north) and Carson River drainage basins and is composed of thick sequences of Cenozoic volcanic rock. The western Virginia Range forms the northern boundary of Eagle and Dayton Valleys, whereas the eastern half of the range forms the northern boundary of Churchill Valley (Moore and Archbold, 1969, p. 21).
The east-trending Desert Mountains form the southern boundary of Churchill Valley and separate the Carson River basin from the Walker River basin (Mason Valley) to the south (fig. 2B). The Desert Mountains are composed of late Tertiary or Quaternary basalts overlying semi-consolidated tertiary rock (Moore and Archbold, 1969, p. 21).
The Dead Camel Mountains, composed of tertiary volcanic rocks, are in eastern Churchill Valley. Despite these mountains forming a distinct topographic boundary, measured groundwater elevations indicate that these mountains do not coincide with a groundwater divide, but rather indicate west-to-east flow through bedrock (Maurer, 2011, p. 37). Fractures in the basalt flows underlying the Dead Camel Mountains (Tingley, 1990) may provide the conduits necessary to sustain groundwater flow through them that provides discharge to the Carson Sink to the east.
The Carson River Basin has been the focus of numerous hydrologic investigations. The efforts of Worts and Malmberg (1966), among the first to appraise water resources in Eagle Valley to provide water resource managers with information pertinent to planning for population growth. Worts and Malmberg (1966) describe hydrogeologic controls in the unconsolidated basin-fill sediments and estimate inflow and outflow components of the Eagle Valley hydrographic area (HA 104; the use of “Eagle Valley” in this report refers to HA 104) water budget, including sources of groundwater recharge and pumpage. In a subsequent study, Glancy and Katzer (1976) expanded the analysis by Worts and Malmberg (1966) spatially by including the Dayton and Churchill Valley hydrographic areas. The hydrologic framework described in these reports, in addition to several maps detailing local geology (Moore, 1961; Bingler, 1977; Trexler, 1977; Pease, 1980), supported the development of an early groundwater model of Eagle Valley (Arteaga and Durbin, 1979; Arteaga, 1986).
Hydrologic studies of Eagle Valley during the 1990s focused on groundwater inflows to the unconsolidated sediments from surrounding watersheds to quantify local water resources better. Using Darcy’s Law, Maurer and others (1996) estimated the total subsurface groundwater inflow from Vicee, Ash, and Kings Canyons ranged between 2,800 and 3,000 acre-ft/yr. This estimate was approximately 1,800–2,000 acre-ft/yr more than what was reported in Worts and Malmberg (1966) and Arteaga and Durbin (1979). Maurer and Berger (1997) expanded the analysis to include all watersheds tributary to Eagle Valley. Total subsurface groundwater inflow to Eagle Valley was estimated to range between 3,200 and 6,100 acre-ft/yr (Maurer and Berger, 1997). The estimated net annual recharge on the floor of Eagle Valley (Maurer and Thodal, 2000, fig. 4), including infiltration from precipitation, irrigation, and streamflow, ranged between 4,300 and 5,900 acre-ft/yr (Maurer and Thodal, 2000, table 9). Together, the estimates described here indicate a total annual groundwater recharge to the Eagle Valley alluvial aquifer between 7,500 and 12,000 acre-ft.
Under predevelopment conditions, groundwater typically moved away from higher elevation mountain-front recharge areas toward points of natural discharge to the east-southeast (Worts and Malmberg, 1966, fig. 4; Arteaga, 1986, p. 6). As a result of pumping for municipal water supply in the vicinity of Vicee Canyon, recent groundwater-level measurements indicated that hydraulic gradients had reversed, and that local groundwater was flowing westward across part of the valley toward these wells (Lopes and others, 2006, p. 18).
The southern boundary of Eagle Valley is next to the northern boundary of Carson Valley (HA 105; fig. 1), for which a model was constructed by Yager and others (2012), referred to here as the “Carson Valley” model. The Carson Valley model describes groundwater flow and exchanges with surface water in Carson Valley. Owing to the presence of Clear Creek on the Carson and Eagle Valley boundary—a perennial stream draining the Carson Range (fig. 2A)—a natural hydrologic divide exists between Carson and Eagle Valleys that is sustained by streamflow loss from Clear Creek.
Hydrologic investigations in the Carson Plains—a subbasin of the Dayton Valley HA (fig. 1)—began with the work of Glancy and Katzer (1976), who reported 7,900 acre-ft/yr of potential groundwater recharge from the entire Dayton Valley hydrographic area (HA 103). Maurer (1997) refined the recharge estimate for Dayton Valley HA to 11,000 acre-ft/yr (not including recharge from the river). Recharge from infiltration of precipitation was estimated by applying the Maxey-Eakin method (Watson and others, 1976) in conjunction with alternative precipitation maps produced by Hardman (1936) and Hardman and Mason (1949), as well as precipitation–elevation relations by Glancy and Katzer (1976) and distance–elevation relations combined with the 1996 Nevada precipitation map (Daly and others, 1994). Estimates of infiltration from precipitation and streamflow in the Dayton Valley subbasins ranged between 5,100 and 16,000 acre-ft/yr in the Carson Plains subbasin, between 560 and 1,500 acre-ft/yr in the Stagecoach subbasin, and between 1,200 and 3,700 acre-ft/yr in the Bull Canyon subbasin (Maurer, 1997). Additional groundwater recharge from seepage from septic systems, lawn watering, agricultural return flow, wastewater-effluent ponds, and losing reaches of the Carson River to these subbasins was also estimated (Maurer, 1997). Although minor in most instances, these additional sources raised the estimated ranges of groundwater recharge for the Carson Plains, Stagecoach, and Bull Canyon subbasins to 7,000–18,000, 1,300–2,200, and 3,600–6,100 acre-ft/yr, respectively (Maurer, 1997). The Carson Plains subbasin in the Dayton Valley HA is the most heavily irrigated region of the three areas included in the study (table 2). Infiltration of flood irrigation contributes substantial recharge to the alluvial aquifer. Explicit estimates of recharge from irrigation return flow are not mentioned in Glancy and Katzer (1976); however, Maurer (1997) reported that 20–40 percent of pumped irrigation return flow was assumed to return as recharge (that is, irrigation return flow), and although not explicitly stated in the report, a similar percentage could be used for surface-water irrigation. Despite irrigation recharge, diversions for irrigation coupled with groundwater extraction in the Dayton Valley hydrographic area have resulted in hydraulic gradients that direct groundwater flow away from the river.
Limited hydrologic investigations are available for the subbasins east of Carson Plains. Harrill and Preissler (1994; p. H12) indicate basin-fill sediment thicknesses of 3,000 ft are reached in an area southeast of Misfits Flat (fig. 2B), and thicknesses range between 500 and 2,000 ft in the rest of the valley. Two cross-sections by Schaefer and Whitney (1992, p. 9) describe basin-fill sediment as deep as 2,900 ft northeast of Churchill Butte and as deep as 2,100 ft southeast of the head of Lahontan Reservoir in Churchill Valley (HA 102; the use of “Churchill Valley” refers to HA 102). Later work by Maurer (2011) suggests thicknesses in these two areas of only 400 and 600 ft, respectively, and a maximum thickness of 1,000 ft in an area approximately 7 mi due west of Silver Springs (approximately 5 mi north of Churchill Butte) based on 736 individual gravity station measurements spread throughout the middle Carson River Basin (Maurer, 2011, fig. 9B; Maurer and Medina, 2020). Schaefer and Whitney (1992) provide maps of groundwater depth, groundwater elevation, and groundwater-flow direction, which in 1982 was eastward toward Lahontan Reservoir; however, water-level declines of approximately 10 ft in two areas within Stagecoach subbasin could be slowing and possibly reversing the direction of groundwater movement (Maurer, 2011, fig. 14C).
Maurer and others (2008) describe streamflow gains and losses between the headwaters and Lahontan Reservoir (figs. 1, 2). For example, the difference between the mean annual river inflow and outflow between the Carson River near Carson City, NV 10311000) and the Carson River at Deer Run Road near Carson City, NV (10311300; fig. 2A) streamgages was 9,200 acre-ft/yr, an amount that approximately equals the annual diversion to the Mexican Ditch plus average groundwater extraction from seven municipal wells in this river reach. The average annual loss of streamflow from the Carson River through Dayton Valley between 1987 and 1992 was 10,000 acre-ft/yr (Maurer, 2011, p. 44) and, from 2000 to 2004, was 9,800 acre-ft/yr. From 1996 to 1998, however, streamflow records indicate a net gain of 22,000 acre-ft annually (Maurer, 2011, p. 44). During years of low precipitation, it was not uncommon for the Carson River to go dry in some or all the reaches flowing through the Dayton Valley HA during summer months (Glancy and Katzer, 1976, p. 36). Tables 3 and 4 summarize many of these, as well as some additional, water-budget components for the Carson River and simulated HAs, respectively.
Table 3.Water-budget summary for select Carson River reaches (U.S. Geological Survey, 2018).
[All reported river budget components calculated using the gaged record time series for water years (WY) 2001–10. Abbreviations: NV, Nevada; USGS, U.S. Geological Survey; —, no data]
Reported values represent average water year volumes (acre-feet per year) for the simulation period WY 2001–10, unless otherwise noted.
Includes Clear Creek at Fuji Park (10310518), Kings Canyon Creek near Carson City (10311100), Ash Canyon Creek near Carson City (10311200), and Vicee Canyon (records obtained from Nevada State Engineer).
Table 4.Summary of groundwater-budget components for Carson River reaches from published studies.
[ET, evapotranspiration; MCR, Middle Carson River; WY, water year; —, no data]
|Pumpage||—||21,330||—||Nevada State Engineer, WY 2001–10|
|Mountain front recharge||31,000; 1,260; 3,800;
|—||—||Estimates from Worts and Malmberg (1966, p. 16); Arteaga (1986, p. 16); Maurer & Thodal (2000); Maurer and Berger (1997, table 8), respectively|
|Seepage from surface-water4||—||2,700||59,200; 67,200||Out: Arteaga and Durbin (1979, p. 32); Net: Maurer and others (2008, p. 53)|
|Phreatophytic ET||—||4,000; 3,000;
|—||Estimates are from Arteaga and Durbin (1979, p. 32); Maurer (1997, table 9), respectively|
|Pumpage||—||2432||—||Nevada State Engineer, WY 2001–10|
|Mountain front recharge||82,200||—||—||Maurer (1997, table 13)|
|Seepage from surface-water9||—||—||10−10,000||Maurer and others (2008, p. 44, 57)|
|Phreatophytic ET||—||11977–2,497||—||Maurer (1997, table 9)|
|Recharge from precipitation||12440–580||—||—||Harrill and Preissler (1994) and Dettinger (1989)|
|Pumpage||—||2122||—||Nevada State Engineer, WY 2001–10|
|Natural recharge||1,300||—||—||Glancy and Katzer (1976, table 17)|
|Seepage from surface-water||—||—||1320,000||Maurer and others (2008, p. 61)|
Estimate based solely on records from Nevada State Engineer and does not include domestic pumping.
Estimates for Riverview subbasin located on the east side of Eagle Valley, but in the Dayton Valley hydrographic area.
Estimate is for entire Dayton Valley hydrographic area, including subbasins not explicitly modeled by MCR model.
From Carson River at Dayton, NV streamgage (10311700) to the Carson River near Fort Churchill streamgage (10312000).
Reported estimate factors in estimated diversions to ditches upstream from the Dayton streamgage (10311700) and Buckland Ditch about 0.1 mile upstream from the Carson River near Fort Churchill, NV streamgage (10312000) and is based on the period 1995–2006.
Among the earliest groundwater models constructed for the middle Carson River study reach, Arteaga (1986) built a two-layer steady-state finite-element groundwater model of the unconsolidated basin-fill sediments in the Eagle Valley HA to assist the Nevada Division of Water Resources evaluate the groundwater system. The model successfully simulated observed groundwater declines due to groundwater development on the western side of Eagle Valley. Model results also showed that groundwater discharge to the Carson River, on the eastern side of Eagle Valley, had remained unaffected through 1978. When projecting conditions by extending the pumping rates for 1978 through to 2000, groundwater declines of 150 and 15 ft were predicted in the western and south-central regions of Eagle Valley, respectively (Arteaga, 1986). Estimated decreases in evapotranspiration and discharge to the Carson River from 1978 to 2000 were 40 and 6 percent, respectively.
Building upon the work of Arteaga (1986), Schaefer and others (2007) built a two-layer steady-state finite-difference model using MODFLOW-2000 (Harbaugh and others, 2000) and used aquifer property values from the Arteaga model to estimate the extent of land in Eagle Valley contributing to recharge of public supply wells.
In a similar study of the Stagecoach subbasin, two valleys to the east, Harrill and Preissler (1994) developed a numerical groundwater-flow model for simulating pre-development (steady state) conditions at an estimated 920 acre-ft/yr of flow through each year: 550 acre-ft/yr of recharge from precipitation plus 280 and 90 acre-ft/yr of recharge from the upper and lower reaches of the Carson River in the Stagecoach subbasin, respectively, balanced by 630 acre-ft/yr of evaporation, 170 and 120 acre-ft/yr of groundwater discharge to neighboring Churchill Valley and the lower reach of the Carson River, respectively. Under post-development (transient pumping) conditions, the model was calibrated to observed changes in water levels from 1971 to 1981. Authors estimated that of the 10,000 acre-ft of net water pumped from the alluvial aquifer between 1971 and 1981, 7,000 acre-ft was derived from storage. Results from 11 alternative pumping scenarios demonstrated that pumping more than the natural discharge (920 acre-ft/yr) led to substantial capture of flow from the Carson River. Thus, pumping rates exceeding natural discharge are sustainable to the extent that reduced surface-water flows in the Carson River are tolerable. Between 1919 and 1969, Brown and Caldwell (2004, p. 7) estimated recharge from seepage loss along the Carson River and Truckee Canal at 268,000 and 170,000 acre-ft, respectively, and groundwater inflow from the Walker Basin totaled 1,000 acre-ft during this period. The uncertainty of this estimate is high because the estimated recharge often exceeds the annual flow in the Carson River and Truckee Canal combined. Maurer and others (2008, p. 60) provided a more realistic estimate of annual streamflow losses at approximately 55,000 acre-ft/yr between 1967 and 2006, which ranged from 13 to 17 percent of the of the annual streamflow measured at the Carson River near Fort Churchill, NV streamgage (10312000).
In additional to numerical modeling focused on simulating groundwater flow in the basin-fill sediments, previous modeling studies also have been developed to simulate streamflow. Hess (1996) and Hess and Taylor (1999) describe an operations model using the Hydrological Simulation Program-FORTRAN (HSPF; Bicknell and others, 1993) built for the upper Carson River Basin (UCRB). The HSPF model simulated streamflow, tributary inflow, and diversions using different flow regimes (that is, storm runoff, low flow) and a daily time step while accounting for the complex operating rules necessary for adherence to the Alpine Decree (U.S. District Court, 1980a, b). Although calibrated against historical streamflow observations, it was designed to facilitate an investigation of the relative effects of various alternative management strategies or water allocations on streamflow and reservoir storages (Hess and Taylor, 1999, p. 36). In a similar river-operation modeling study, Kennedy-Jenks Consultants (1998) used MODSIM (Labadie, 2006) to simulate the legal and physical components of the UCRB and evaluate the effects on annual inflow to Lahontan Reservoir from alternative management strategies that replaced agricultural consumptive-use demands with municipal and industrial consumptive-use demands.
MODFLOW-NWT (Niswonger and others, 2011) was selected for simulating groundwater flow through the basin-fill sediments of Eagle, Dayton, and Churchill Valleys. As do all versions of MODFLOW, MODFLOW-NWT supports the use of independent, modular components called “packages” to simulate pertinent aspects of a groundwater-flow system. In addition to this, MODFLOW-NWT overcomes model-cell wetting and drying problems that troubled unconfined groundwater-flow simulations in earlier versions of MODFLOW (for example, MODFLOW-2000, MODFLOW-2005, and others). By keeping dewatered cells active, the non-linearity associated with cell drying and rewetting are avoided in simulations including a falling and rising water table, which is common in the MCR.
Using MODFLOW packages, MODFLOW-NWT accounts for (1) spatially variable aquifer properties, including horizontal hydraulic conductivity, ratios of vertical to horizontal hydraulic conductivity, and specific yield; (2) recharge along the lateral boundaries of the model; (3) infiltration into the unsaturated zone with subsequent partition of unsaturated flow to recharge, evapotranspiration, and changes in unsaturated-zone storage; (4) pumping from wells spanning multiple model layers; (5) groundwater and surface-water interaction between the alluvial aquifer and the mainstem Carson River, perennial and ephemeral tributaries, irrigation ditches, and a large surface-water reservoir; (6) calculation of hydraulic heads at discrete points; and (7) calculation of surface-water flows for specific reaches of the surface-water network.
In many groundwater modeling applications, especially involving unconfined aquifers and groundwater–surface-water interaction, achieving numerical model convergence depends on finding a solution to a system of nonlinear equations. MODFLOW-NWT uses the Newton method for iterating toward a solution; previous solvers available with MODFLOW-2005 (Harbaugh, 2005) used Picard methods. Advantages of MODFLOW-NWT include robust handling of nonlinear boundary conditions created by the Multi-Node Well (MNW2; Konikow and others, 2009), Streamflow-Routing (SFR2; Prudic and others, 2004; Niswonger and others, 2006), and Lake (LAK; Merritt and Konikow, 2000) packages, each of which is active in the present model application. Thus, the use of MODFLOW-NWT, and specifically the χMD solver (Niswonger and others, 2011, appendix C), in the MCR model application greatly improves model convergence and computational efficiency compared to the solvers available with previous versions of MODFLOW.
Upstream Weighting Package
MODFLOW-NWT relies on the Upstream Weighting (UPW) package for calculation of inter-cell conductance terms for solving the discretized groundwater-flow equation (Niswonger and others, 2011). Calculation of the horizontal conductance is smoothed in MODFLOW-NWT by using combined quadratic and linear functions for increased stability as wetting and drying of cells were simulated. In this way, the UPW package eliminates the nonlinearities often responsible for non-convergence problems when using the Block-Centered Flow (BCF), Layer Property Flow (LPF), and Hydrogeologic-Unit Flow (HUF) packages (Anderman and Hill, 2000; Harbaugh, 2005). A more detailed description of the UPW package can be found in Niswonger and others (2011, p. 4).
The Unsaturated-Zone Flow (UZF1; Niswonger and others, 2006) package, available for use with MODFLOW-NWT, was used for simulating variably saturated flow conditions above the water table. The UZF1 package provides “an efficient means of simulating recharge” while also accounting for evapotranspiration (ET), storage, and flow effects (for example, drainage) in the unsaturated zone (Niswonger and Prudic, 2004; Niswonger and others, 2006). The UZF1 accomplishes this by using the method of characteristics to solve a kinematic wave equation for unsaturated flow by neglecting the diffusive term necessary with Richards-based solutions of variably saturated flow. Simulating vertical heterogeneity is not possible with the UZF1 package; only one value for each hydraulic property (for example, vertical hydraulic conductivity, Brooks-Corey epsilon) can be specified for each vertical column of cells (row and column combination).
The UZF1 package facilitates specification of stresses at land surface (for example, precipitation, applied irrigation, potential ET), where stresses are more readily observed and quantified, thereby alleviating the need to specify recharge (that is, recharge is calculated rather than specified). To calculate the recharge that results from infiltration, appropriate UZF1 input arrays must be specified (table 5). The input array IRUNBND directs the routing of groundwater discharge to land surface and of groundwater infiltration in excess of saturated vertical hydraulic conductivity (Niswonger and others, 2006). In the MCR simulation, groundwater discharge and infiltration at rates in excess of saturated vertical hydraulic conductivity values are routed to the nearest stream segment (or lake, for areas surrounding Lahontan Reservoir). The UZF1 package calculates the residual water content internally by subtracting the specific yield from the saturated water content (Niswonger and others, 2006, p. 6), thereby alleviating the need to specify residual water content.
Table 5.Required input arrays for the UZF1 package and their range of values. (Niswonger and others, 2006)
[—, no data; L, Length; T, Time; in., inch; ET, evapotranspiration; in/day, inch per day; ft, foot]
Construction of the Groundwater Flow Model
The MCR groundwater and surface-water model was built using the MODFLOW-NWT (Niswonger and others, 2011) simulation software, a standalone program that is based on the MODFLOW-2005 (Harbaugh, 2005) groundwater-flow simulator. MODFLOW-NWT uses the Newton solution method to overcome nonlinearities associated with the wetting and drying of computational grid cells while solving the unconfined groundwater-flow equation. MODFLOW-NWT relies on several input datasets that describe the groundwater and surface-water system, including the hydrogeologic units, unsaturated zone, surface-water features, initial conditions, transient boundary conditions, and water use.
Two primary hydrogeologic units are accounted for in the MCR model: (1) unconsolidated basin-fill sediments and (2) semi-consolidated permeable volcanic rock of tertiary age that underlies (and bounds) the basin-fill sediment. Pre-Cenozoic consolidated rock that underlies these units forms a no-flow boundary at the bottom of the model.
Spatial and Temporal Discretization
MODFLOW-NWT simulates three-dimensional movement of groundwater of constant density through porous media using a finite-difference approximation of equation 2–1 in Harbaugh (2005). To solve this equation, the spatial domain was horizontally and vertically discretized into rectilinear blocks commonly referred to as cells. Aquifer properties assigned to a cell are assumed to be homogenous throughout each cell. The grid used to discretize the MCR model domain is directly related to the model-grid discretization applied by Yager and others (2012, p. 29) to simulate groundwater flow in the Carson Valley, next to and south of the area of the MCR model (fig. 8). Identical grid discretization was applied to the MCR model to facilitate merging of the MCR with the Carson Valley models (Yager and others, 2012), should a unified model prove to be more useful in the future. Grid cells are 550 ft on each side throughout the modeled region and spatially aligned on a north–south axis with the grid documented in the Carson Valley model (Yager and others, 2012). Cell sizes are small enough for adequate simulation of groundwater and surface-water exchange, yet large enough to keep numerical computational times reasonable. The model grid consists of 376 rows and 546 columns. The model domain covers an area of approximately 389.7 square miles (mi2, or 249,400 acres), 278.1 mi2 (178,000 acres) of which is unconsolidated basin-fill sediment. The remaining area, 111.6 mi2 (71,500 acres), includes a part of the Dead Camel Mountains.
The lateral and vertical extent of the basin-fill sediment is based on the location of the permeable bedrock units that bound (mountains surrounding the valleys) and underly the unconsolidated basin-fill sediment. Vertically, the finite-difference grid (model domain) was divided into six layers. Layer 1 was used to represent Lahontan Reservoir. Layers 2–5 represent the unconsolidated basin-fill sediments, and each layer contains 25,627 active cells. Layer 6 represents the semi-consolidated rock underlying the basin-fill sediments. Two slug tests were done where the semi-consolidated rock unit outcrops, and as expected, hydraulic conductivities of less than 1 ft/day in this unit were less than those in the unconsolidated basin-fill sediments, where hydraulic conductivities are typically at least an order of magnitude greater (Maurer, 2011; table 1). Owing to the possibility of fracturing within this unit, however, the existence of highly transmissive zones within this unit remains an open question (Maurer, 2011, p. 27). In addition to underlying the unconsolidated basin-fill sediment, layer 6 also was used to represent the volcanic outcropping of the Dead Camel Mountains, where the hydraulic conductivity is considerably less than the neighboring basin-fill sediment, although still permeable owing to fractures, as discussed previously. Altogether, layer 6 hosts 35,916 active cells. This model characterization results in a total of 138,424 active cells in the model domain. The areal extent of active cells can be seen in figure 8.
Land-surface elevation was calculated using the zonal statistics tool available with Spatial Analyst for ArcGIS (Environmental Systems Research Institute, 2012) 10- and a 30-meter digital elevation model (DEM; U.S. Geological Survey, 1999). Except for the cells overlying the Dead Camel Mountains outcrop, the mean land-surface elevation was assigned to the top of layer 2. In the Dead Camel Mountains, land-surface elevation was assigned to the top of layer 6. Layer 2, where most of the surface-water network is simulated, was assigned a constant thickness of 20 ft to minimize the elevation difference between cell centroids and subgrid surface-water features found in many of the layer 2 cells. Through finer vertical discretization near surface-water features where vertical flow is concentrated, model error related to head and groundwater–surface water interaction can be reduced. In addition to vertical discretization, streambed hydraulic conductivities were kept low to provide resistance to flow as a safeguard against systemic bias of groundwater–surface-water exchange (Anderson, 2005; Markstrom and others, 2008).
The thickness of layers 3 through 5 varied throughout the model domain. Near the perimeter of the active cells, where the unconsolidated basin-fill sediments taper-off at the surface contact with consolidated rock, layers 3 through 5 were assigned a minimum thickness of 13.3 ft to avoid unnecessarily thin cells. Moving toward the center of the modeled subbasins, the thickness of basin-fill sediments increased with increasing distance from the outcropped consolidated rock. Likewise, layer thickness was increased to capture the full vertical extent of the basin-fill sediments estimated in Arteaga (1986, fig. 11) and Maurer (2011, fig. 9). The thicknesses of layers 3 through 5 were uniformly increased until the predetermined maximum thicknesses were reached: 40 ft for layer 3; 60 ft for layer 4; and the remaining thickness necessary to span the vertical extent of the basin-fill sediments (up to 1,880 ft) was assigned to layer 5. Thus, four progressively thicker layers were used to represent the unconsolidated basin-fill sediments. A similar approach is documented in Allander and others (2014), which served as the framework for this modeling effort. Layers 2 through 6 were simulated as convertible, meaning every cell in the active domain could dewater or contain the water table. Vertical elevations of the cell (layer) interfaces are recorded in the North American Vertical Datum of 1988 (NAVD 88).
Groundwater interaction between the unconsolidated basin-fill sediments and the semi-consolidated rock that underlies the modeled region is allowed in the model, but layer 6 hydraulic conductivities are several orders of magnitude less than those in the layers representing the unconsolidated basin-fill sediment. Hence, fluxes between layers 5 and 6 are small. Consolidated rock of pre-Cenozoic age that underlies layer 6 is assumed to be a no-flow boundary. The bottom of layer 6 was set equal to 3,000 ft above mean sea level. This elevation roughly corresponds to the average estimated elevation of the top of the pre-Cenozoic rock unit below the semi-consolidated rock units represented by layer 6 of the MCR model (Maurer, 2011, fig. 6).
The 11-year simulation, from January 1, 2000, to December 31, 2010, is divided into 574 weekly stress periods, with 1 time-step for each stress period. A steady-state stress period precedes the transient part of the simulation, however, to calculate the initial conditions priming the transient simulation. For the MCR model, steady-state stresses represent the average conditions during the transient part of the simulation and set model values (for example, groundwater elevations) at the beginning of the transient period roughly equal to observed conditions in early 2000. For example, the steady-state Carson River inflow is equal to the weekly inflow averaged across the full transient simulation. Although week-long stress periods were selected, a number of the stresses simulated by the model reflect shorter or longer lengths of time. For example, pumping rates used in the model change on a monthly basis because this was the frequency of data collection historically for pumped amounts supplied by the Nevada State Engineer’s Office. In contrast, surface-water inflows entering the model (for example, the Carson River on the south side of Carson City) vary by day and are averaged by MODFLOW according to the specified stress periods (weekly in the case of the MCR model).
The 11-year period selected for transient simulation and calibration facilitates assessment of model performance throughout a range of simulated groundwater elevations and streamflow representing drought as well as flood conditions. During the simulated period, the annual streamflow measured at the Carson River near Fort Churchill streamgage (10312000) ranged from as little as 99,000 and 102,000 acre-ft, in 2007 and 2008, respectively to as much as 535,000 acre-ft in 2006 (fig. 9). The 2006 and 2007 runoff totals represent about half of the average annual streamflow runoff for the period of record. In 2005, the annual streamflow was approximately 324,000 acre-ft, or more than 2.5 times larger than the average annual streamflow runoff for the period of record. In the remaining years, runoff totals ranged between −106,000 (2001) and −228,000 acre-ft (2010). It is worth noting that the average annual streamflow at the Carson River at Fort Churchill streamgage (10312000) during the period used in the transient simulation was 25 percent less than the average annual streamflow for the entire period of record (fig. 9).
Two distributed parameter fields necessary for groundwater modeling are hydraulic conductivity and specific yield. The hydraulic conductivity of geologic material describes the rate of groundwater movement through a given cross-sectional area of geologic material (basin-fill sediment or semi-consolidated rock). Specific yield is the ratio of the volume of water that gravity drains from a unit volume of aquifer as the groundwater level drops (Domenico and Schwartz, 1998).
For the MCR model, initial assignment of hydrologic property values was based on values reported in Maurer (2011). Computed hydraulic conductivities from 28 slug tests completed in the unconsolidated basin-fill sediment ranged from as little as 1 ft/day to as much as 300 ft/day, and values closer to the river were generally greater than those outside the flood plain. The lowest estimates of hydraulic conductivity, both less than 1 ft/day, were from two slug tests in wells drilled into the semi-consolidated hydrogeologic unit (sandwiched between the lower pre-Cenozoic and upper unconsolidated basin-fill sediment units) between the Carson Plains and Stagecoach subbasins of Dayton Valley.
Initial estimates of specific yield used in the MCR model were taken directly from delineations shown in figure 13 of Maurer (2011) and originally mapped by Stewart (1999). Unconsolidated basin-fill sediments were sub-divided into four general classifications: (1) Tertiary sediment; (2) alluvial fans; (3) lake deposits; and (4) basaltic rocks. In these four classifications, specific yield is expected to range between 1 and 40 percent (see Maurer, 2011, for a more detailed breakdown). In the semi-consolidated (that is, fractured) volcanic rock unit, specific yields range between 1 and 15 percent.
Whereas hydraulic conductivity and specific yield were varied during the calibration process, values of vertical anisotropy and specific storage remained fixed at their initial values of 0.5 and 1.0e–08, respectively (Domenico and Schwartz, 1998). A value of 0.5 for vertical anisotropy means that the vertical hydraulic conductivity is half of the horizontal hydraulic conductivity, which was the condition specified to represent basin-fill sediment free of inter-bedded confining units. Specific storage describes the amount of water that is released from storage because of lowering the groundwater head without dewatering. In effect, the amount of water released from storage as pressures are reduced is equivalent to the volumetric expansion of groundwater and simultaneous contraction of pore space (Domenico and Schwartz, 1998). Although this form of storage release can occur in unconfined aquifers, it is negligibly small relative to the volumes of water that are recovered from drainage of pore space.
Boundary conditions are the mathematical representations of where and how water enters and leaves the groundwater system (Anderson and others, 2015). Boundary conditions often are divided into three principal categories for describing how water may enter or exit the groundwater system: (1) specified-head (constant heads where data are sparse); (2) specified-flux (pumpage, recharge from the mountain block, stream inflows, and surface-water diversions); and (3) head-dependent flux. Head-dependent fluxes dynamically account for flow direction (into or out of the aquifer) and magnitude depending on relative differences in groundwater head and stream stage, lake stage, ET, and recharge from irrigation and precipitation. For example, in addition to water entering and exiting the model along its perimeter, water can also enter and exit the groundwater system through different stresses that cause hydraulic-head changes within its borders, like pumping.
In Eagle, Dayton, and Churchill Valleys, different combinations of sources and sinks (boundary conditions) are required for each HA to adequately simulate the cumulative effects of locally important stresses. Figures 10–15 highlight the boundary conditions described next, using simplified block diagrams (figs. 10, 12, 14) and detailed two-dimensional maps of the numerical grid showing the various boundary conditions depicted in different colors (figs. 11, 13, 15).
Three types of specified-flux boundaries were simulated in the MCR model: (1) mountain-front recharge; (2) groundwater pumping; and (3) streamflow across the boundary of the active model domain (figs. 11, 13, 15).
The amount of groundwater inflow from mountain blocks surrounding the model domain and specified in the model is based on the findings of previous studies (Maurer and Berger, 1997; Jeton and Maurer, 2011). Mountain-front recharge includes groundwater inflow from two sources—the groundwater inflow entering the active model domain through the alluvial ‘fingers’ extending up the perennial and ephemeral drainages and the groundwater flow across the contact between the mountain block and basin-fill sediments. Each of these two sources was accounted for separately but appear in the same input file (RCH; Morway and others, 2023b). Model setup ensured that active UZF1 cells (IUZFBND=1) did not coincide with active RCH cells to avoid conflicts in the numerical solver. In Eagle Valley, specified recharge along the mountain block and unconsolidated basin-fill sediment interface was specified on the west, north, and east sides of the valley (fig. 11). In Dayton (specifically, Carson Plains subbasin) and Churchill Valleys, recharge was specified along the north and south sides of the active model domain (figs. 13, 15).
Values of recharge from groundwater inflow originating in the alluvial ‘fingers’ extending up the perennial and ephemeral drainages are from table 4 of Jeton and Maurer (2011). Groundwater inflows of this type ranged from 40 acre-ft/yr (Eureka Canyon) to nearly 5,000 acre-ft/yr (Churchill Canyon). All values of groundwater inflow from the various drainages were allowed to vary during PEST calibration (Doherty, 2003), but were “regularized” (that is, preferred value was specified in the PEST input file) using values determined in Jeton and Maurer (2011).
Recharge from the mountain block was derived from estimates of annual rainfall calculated by the Parameter-elevation Regressions on Independent Slopes Model, or PRISM (Daly and others, 1994; PRISM Group, 2012). The average annual precipitation calculated by PRISM was resolved to an average annual volume by integrating the spatially variable rainfall over each catchment next to the active model grid. Recharge was estimated by applying mean infiltration efficiencies (Jeton and Maurer, 2011, table 4) to the average annual precipitation volume, converting to a daily rate, and distributing these values evenly along the mountain-block and basin-fill sediment contact. The PEST (Doherty, 2005) provided the means to vary estimated rates for mountain-front recharge.
Groundwater pumping from the unconsolidated basin-fill sediments was simulated in wells for which a historical record of pumped amounts is maintained by the Nevada State Engineer (figs. 10–15). Pumping from domestic wells was not simulated. Domestic wells pump very small volumes of water compared to irrigation and municipal wells and are assumed to have negligible effect on the groundwater system. The monthly pumped volumes obtained from the Nevada State Engineer (Morway and others, 2023b) were converted to an equivalent daily rate for assimilation into the appropriate input file. Locations of the pumping wells are shown in figures 11, 13, and 15.
The multi-node well package (MNW2; Konikow and others, 2009) was used to simulate groundwater pumping from wells. The MNW2 package accepts the well-screen intervals as one of its optional inputs. Well-screen intervals were known for each of the simulated wells. Using this option alleviates the need to specify from which layer groundwater is withdrawn; MODFLOW automatically makes this determination by well screen and model layer elevations. After the final model run of the baseline simulation, pumped amounts were evaluated to ensure that the specified pumped amounts were fully realized. This step was necessary because MODFLOW-NWT automatically reduces pumping to avoid numerical (model) instability associated with pumping from dry grid cells.
Estimates of surface-water inflows at the active MCR model domain were specified based on records from the Truckee Canal near Hazen, NV streamgage (10351400), shown in figure 2B and from the following streamgages shown in figure 11: (1) Clear Creek near Carson City, NV (10310500); (2) Carson River near Carson City, NV (10311000); (3) Kings Canyon Creek near Carson City, NV (10311100); (4) Ash Canyon Creek near Carson City, NV (10311200). In addition to the gaged inflows, estimates of surface-water inflow from Vicee Canyon were specified according to the quarterly release records maintained by the Nevada State Engineers Office (Matt Dillon, Nevada State Engineers Office, written commun., May 8, 2012). All specified inflow entering Vicee Canyon seeps into the ground soon after entering the model domain, which corresponds well to the physical system. Surface-water inflows for the ephemeral watersheds of the MCR were assumed to be negligible based on considerable indirect evidence of limited intermittent flows (Jeton and Maurer, 2011, p. 28). Even during wet years, tributary inflow remains a fraction of the overall Dayton (Carson Plains subbasin) and Churchill Valleys water budget (Maurer and others, 2008; Jeton and Maurer, 2011). The only surface-water flow exiting the MCR model domain is the Carson River downstream from Lahontan Dam. Estimates of surface-water flow from Lahontan Reservoir are based on the USGS streamgage 10312150 about 1-mile downstream from Lahontan Reservoir plus estimates of flow into Rock Dam Ditch, which diverts water upstream from the gage.
Surface water diverted from the Carson River is represented as specified fluxes to 10 irrigation ditches simulated in the MCR model (table 2; figs. 13, 15A, 15B). Seasonal (April–October) diversions were specified according to historical records obtained from the Federal Water Master in Reno, Nevada (Dave Wathen, Federal Water Masters Office, written commun., July 6, 2011). The SFR2 package, used to simulate surface-water flow, calculates the reduction in main-channel flow by an amount equal to that diverted at the point of diversion. Unused surface water diverted to the irrigation ditches returns to the main channel of the Carson River and is simulated as such.
Head-Dependent Flux Boundaries
Estimates for two head-dependent model inputs, ET, and irrigation applications, are based on a spatial delineation of riparian and irrigated areas. The riparian area bounding the Carson River was visually interpreted from 1-meter resolution National Agricultural Imagery Program (NAIP) imagery from 2010 (U.S. Department of Agriculture, 2012) and is spatially contained in those areas modeled with the UZF boundary type (figs. 11, 13, 15). The riparian-area boundary was delineated from the southernmost extent of the model domain, south of Eagle Valley, to the confluence with Lahontan Reservoir as defined by the National Hydrography Dataset (McKay and others ; http://nhd.usgs.gov/, accessed June 2011). The riparian area consists of vegetation typical of the western Great Basin such as Freemont cottonwood (Populus fremontii), salt cedar (Tamarix ramosissima), and willow (Salix spp.) interspersed with other grasses and shrubs. For the purposes of this delineation, riparian areas were defined as the vegetated area that appears in the NAIP imagery as a band of moderately lush, green vegetation bounding the banks of the Carson River. In much of the study area, the relatively dense, healthy vegetation of the riparian zone visually contrasts in the imagery with the sparser and drier native sage and scrublands and with the urban lands along the riparian corridor. This contrast was used to guide digitization of the boundary in a Geographic Information Systems (GIS) dataset. At the southern end of Eagle Valley and northern Carson Valley, the boundary between agricultural field boundaries and the riparian area is less distinct, so the riparian area was generalized to the area where meander scars from the Carson River are clearly visible in the imagery. This area does not fall in the active model domain but was digitized to ensure all active model cells would be included in the riparian lands dataset. Irrigated lands on the boundary or interior of the riparian zone in the active model domain were classified as agricultural lands and not included in the riparian area classification.
The delineation of agricultural lands for the 11-year simulation period were modified from field boundaries determined from GIS data supplied by the Nevada State Engineer’s office (M. Dillon, Nevada State Engineers Office, written commun., August 4, 2011). The field boundaries were compared to the 2010 NAIP imagery and modified as needed to match the imagery. Where irrigated fields were urbanized during the study period, water application ceased in the model at the appropriate time. Fields were added where visible in the imagery, but none existed in the source data.
The modified field-boundary data were compared against a series of satellite images acquired by the Thematic Mapper (TM) instrument aboard the Landsat 5 satellite (not the same as the NAIP imagery). The TM instrument collects information in seven spectral bands at wavelengths ranging between the visible blue (0.45 micrometer [μm]) to the thermal infrared (12.5 μm). The TM imagery is delivered as “scenes” consisting of one image for each spectral band. Each scene is imaged by the sensor every 16 days at 30-meter spatial resolution and covers approximately 31,110 square kilometers (km2). Scenes are identified by a unique path and row number. The study area falls entirely within path 43, row 32. Available TM scenes were evaluated using the U.S. Geological Survey Earth explorer web service (https://earthexplorer.usgs.gov/). One scene was acquired for each month of June, July, and August for 2000 through 2010 (table 6) if cloud-free scenes were available. The satellite images were corrected to top of atmosphere reflectance using standard formulas and calibration constants (Chander and others, 2009) and converted to normalized difference vegetation-index images (NDVI; Rouse and others, 1974). The NDVI is calculated as the difference between reflectance in the near infrared and red wavelengths measured by the satellite instrument divided by the sum of the same wavelengths. The index values are unitless and range between –1 and 1 such that higher values are indicative of relatively healthy or vigorous vegetation. Summertime NDVI values for irrigated agricultural lands in the MCR Basin commonly range between 0.5 and 0.95. Each agricultural field in the modified field boundary dataset was classified irrigated, not irrigated, permanently dried, riparian transition, or urban for each year between 2000 and 2010 as based on visual inspection of all the summer scenes for that year. In the event a single scene within the suite of scenes collected for a given year showed NDVI values greater than 0.5 for more than 50 percent of the field, the field was classified as irrigated for that year. Permanently dried fields are those that did not appear to change land use, but where irrigation appeared to cease at some point of the simulation period. Fields were classified as urban if they were converted to urban land use during the investigation. Riparian transition denotes fields that appeared to be in transition from agricultural use back to a natural riparian state, including fields used for grazing (table 7).
Table 6.Dates of collected Landsat scenes used to determine crop-cover classification (U.S. Geological Survey Earth explorer web service. Available at https://earthexplorer.usgs.gov/).
[mm/dd/yyyy, month/day/year; TM, thematic mapper]
Table 7.Summary of land classification total areas (in acres) for each year of the simulation.
[Land areas were estimated using a combination of Landsat imagery (irrigated vs. not irrigated) and aerial photography (permanently dried vs. urbanized vs. riparian transition)]
The UZF1 package was used to simulate ET along the riparian corridor and in agricultural areas of the MCR model (figs. 2, 11, 13, 15). This approach allows for the full ET demand to be met either by soil moisture stored in the unsaturated zone or by ET directly from the water table, or by some combination of these two. Precipitation on the valley floors was determined using historical observations from the Carson City and Lahontan Reservoir gages maintained by the Western Regional Climate Center (2013). The total annual ET demand exceeds annual precipitation on the valley floors, with the net result that only minor volumes of infiltrated water, if any, contribute to valley-floor recharge (Maurer and others, 2008). When the ET demand exceeds the available soil moisture and the groundwater elevation is above the extinction depth, UZF1 attempts to satisfy the remaining ET demand directly from groundwater using the same method as in the MODFLOW ET package (McDonald and Harbaugh, 1988, p. 10-1).
Grass and alfalfa hay are the most planted crops on lands serviced by 10 ditches in the MCR study (figs. 11, 13, 15) area. Simulated irrigation practices in the study area, including historical diversions, were based on communication with the Federal Water Master in Reno, Nev. (Dave Wathen, Federal Water Masters Office, written commun., April 2, 2012). A customized algorithm was used to distribute the water historically diverted among the fields served by each ditch. This approach avoids simulation of irrigation applications during intervals when diversion records indicate a ditch was dry. The GIS routines resolved the differences between fields and cells to apply the irrigation water to the correct cells. Once irrigation water was specified in the model, the UZF1 package partitioned this water to model-calculated infiltration, runoff, unsaturated-zone storage, and recharge.
Recharge and groundwater discharge resulting from irrigation applications were assigned to the highest active model layer in each active vertical column. Because of the proximity of the irrigated areas to the Carson River, where a shallow water table is sustained, this generally resulted in recharge and discharge to and from model layer 2, the uppermost active layer the water table could rise to (model layer 1 represents Lahontan Reservoir).
As described in the “Surface-Water Network” section, inflow laterally crossing the active model domain is simulated using specified-flux boundary conditions. Upon entering the active model domain, simulated streamflow moves down-channel and is lost to seepage or replenished by groundwater. This stream–aquifer interaction was simulated using the SFR2 package for the entire MCR surface-water network, including unsaturated flow beneath intermittent (for example, irrigation ditches) and ephemeral channels (Niswonger and Prudic, 2005). The SFR2 network consists of 71 segments divided into a total of 1,706 reaches (that is, reaches of channels in a model grid cell). The streamflow network of the MCR model is shown in figures 11, 13, and 15. The streamflow network includes 10 diversions (table 2), 3 perennial tributaries (Clear Creek, Kings Canyon, Ash Canyon), 2 controlled releases into the active model domain (Vicee Canyon and Truckee Canal), 10 ephemeral tributaries (Jeton and Maurer, 2011), and the mainstem of the Carson River.
The quantity of groundwater–surface-water interaction (volume per time) in each streamflow routing cell is determined as the product of a streambed conductance (area per time) and the difference between the water-table level and stream stage (L) in that cell. Calculation of the conductance term C is the same as that used for the standard River Package (RIV; McDonald and Harbaugh, 1988, chap. 6):
is the hydraulic conductivity of the streambed material (length per time),
is the length of the stream reach,
is the width of the stream reach, and
is the thickness of the streambed.
Initial values of Kstr were kept within the ranges reported in Stonestrom and Constantz (2003, appendix B, table 1) and Conlon and others (2003, fig. 8), but were varied during the automated parameter-estimation calibration procedure. The L was determined using GIS routines that calculated the length of channel in each cell. The l was set to 1 ft for all reaches in the SFR2 network, an approach used in similar modeling studies (Ely and Kahle, 2012, p. 28; Yager and others, 2012, p. 37). The final parameter necessary for calculating conductance, W, was calculated by SFR2 according to streamflow for the current time step; hence, the conductance varies during the simulation as a function of stream wetted width, itself varying with streamflow.
For better approximation of the wetted stream width, a generic eight-point description of a cross-sectional channel geometry was used to represent W more realistically (Prudic and others, 2004, fig. 4). For the Carson River, a generic geometry with a linear width (not accounting for cross-section bed slopes) of up to 200 ft was used. In the other natural drainages (tributaries), a generic geometry with a linear width of up to 50 ft was used. Other available alternatives, such as using a fixed rectangular geometry, would not represent the non-linear nature of groundwater–surface-water interaction as wells. Conversely, rectangular geometries with fixed widths of between 7 and 10 ft were assumed for the 10 irrigation ditches because the diversion records indicated relatively steady flows. In addition to varying stream wetted width, the eight-point cross section provides information on variable stream depths throughout the simulation period for improved representation of spatial and temporal variation affecting groundwater–surface-water interactions. Stream depth was calculated using a constant value for Manning’s equation (Prudic and others, 2004) equal to 0.04 (Cimbala and Çengel, 2008).
Streambed elevations in the MCR model were determined using a combination of light detection and ranging (lidar; BAE SYSTEMS Advanced Technologies, 2004) and 10-meter digital elevation model (DEM; U.S. Geological Survey, 1999) products for the Carson River and all other (tributary and diversion) channels, respectively. A DEM is a representation of topographic elevation that uses a grid of square cells each with an associated value that is the average elevation of the area covered by the cell. The streambed elevation was determined by calculating the mean elevation for each active model grid cell and then applying a constant value to every channel that represents the depth of stream incision. The lidar data were available for the entire Carson River reach but did not provide a continuous dataset for calculating streambed elevation of tributaries and diversions. For tributaries and diversions, notable differences between the two elevation datasets were found, with up to 20 ft differences in parts of the Carson Plains area. Thus, some inaccuracy in groundwater–surface-water exchange is attributable to the uncertainty associated with estimating streambed elevations. Because the active model domain is restricted to valley floors (except for the Dead Camel Mountain) that have gentle topographic relief and hydraulic gradients, the magnitude of the resulting errors of this type is assumed to be small.
Surface flows enter Lahontan Reservoir through the Carson River and Truckee Canals. Managed releases from the reservoir are controlled at the dam site (fig. 15B), as described in the “Surface-Water Network” section. Lahontan Reservoir is an important feature affecting groundwater–surface-water interaction in Churchill Valley (Maurer and others, 2008, p. 83). Depending on the prevailing hydrologic conditions each year, 70-ft swings in reservoir stage are possible. Such large swings substantially alter the gradients between the reservoir stage and surrounding water-table levels. For example, Maurer (2011, fig. 20) documents groundwater seepage filling the thalweg of the reservoir (the reservoir was nearly empty) during an extended period of no inflow to the reservoir from the Carson River. During full stage, gradients are reversed, and the water table in monitoring wells approximately 1 to 1.5 mi west of the reservoir showed a corresponding rise that peaked approximately 6 weeks after reaching full stage (Maurer and others, 2008, fig. 27; Maurer, 2011, fig. 19). Thus, groundwater exchange with Lahontan Reservoir varies depending on the time of year, location within the boundaries of the reservoir, and the reservoir stage relative to the surrounding groundwater elevations.
The LAK package (Merritt and Konikow, 2000) was used to simulate the groundwater–surface-water interaction between Lahontan Reservoir and the surrounding area. The LAK package updates a water budget for Lahontan Reservoir that is independent of the water budget maintained for the aquifer. The LAK package allows the spatial extent of grid cells representing exposed or submerged areas of Lahontan Reservoir to vary as a function of lake stage. As lake stage falls, groundwater elevations may be left higher than the lake stage, in which case the MCR model simulates groundwater discharge to the lake. Conversely, as lake stage rises, the model simulates surface-water seepage from the lake to the surrounding aquifer. The flow of water through the lakebed material (seepage of groundwater from the aquifer to the lake or from the lake to the aquifer) is controlled in the model by differences in hydraulic heads in the aquifer and lake stage and by the lakebed leakance parameter (Merritt and Konikow, 2000). Furthermore, volumes of water lost to evaporation from the surface of the lake, or gained through precipitation falling on the lake, vary according to a stage–surface area relationship provided by the user (fig. 16). At capacity, the surface area of Lahontan Reservoir is approximately 13,000 acres. Evaporation rates from the surface of Lahontan Reservoir were specified according to the findings of Huntington and McEvoy (2011, table 5) and varied monthly, from 0.9 inches in February to a maximum of 8.1 inches in August, totaling 50.0 inches each year. As did the evaporation rates, precipitation inputs across the surface boundary of the reservoir changed with each season and were approximately 5.5 inches each water year. The reported precipitation and evaporation depths are independent of each other and, therefore, represent total annual depths. Depending upon the surface area of the reservoir in each simulated time step, total volumes of evaporation and precipitation vary across all months (that is, no 2 months would be expected to have the same volume of evaporation or precipitation).
To represent the geometry of Lahontan Reservoir in MODFLOW, GIS datasets of the 2004 Lahontan Reservoir Bathymetric Survey (Ferrari, 2005) were acquired from Reclamation (Jeff Reiker, written commun., October 2011). The contour lines contained in the acquired datasets represent the underwater (bathymetric) topography of the reservoir and above water (bank) topography derived from digital aerial data. Contour elevations in the data were reported in Lahontan Stage Datum (LSD). The LSD is a local datum that measures 3.75 feet above National Geodetic Vertical Datum of 1929 and 0.30 feet above NAVD 88. The data were imported to geodatabase format (Environmental Systems Research Institute, 2012) and projected from State Plane Nevada West (FIPS 2703) to Universal Transverse Mercator Zone 11 North. Elevations for both datasets were in feet. The NAVD 88 elevation values were added to the attribute table and calculated for the contour lines.
The bathymetry and bank contours were used to create a DEM representing the bathymetric topography of the reservoir (Morway and others, 2023a). The elevation 4,164.7 feet above NAVD 88 (4,165 feet above LSD) was selected as the outer boundary of the bathymetric DEM. Bank contours representing the bounding elevation were selected from the source bathymetry dataset to create the boundary. The selected contours were not continuous around the perimeter of the reservoir, so missing segments of the contour were filled using contours derived from the 1/3 arc-second National Elevation Dataset (NED, U.S. Geological Survey, 1999). A single raster cell in the 1/3 arc-second NED has a spatial resolution of approximately 10 x 10 meters. Bathymetric and bank contours that were completely within the final bounding contour were selected from the source data. Only bank contours that did not intersect the bathymetric contours were included in the selection. The selected data and the bounding contour were merged to a new dataset and used to create the bathymetric DEM. The ArcGIS Topo-to-Raster tool (Environmental Systems Research Institute, 2012) was used to interpolate a 10-meter spatial resolution bathymetric DEM from the merged contour data. The bounding contour was used to create a mask that was used to limit the extent of the bathymetric DEM. The 10-meter bathymetric DEM was merged with the 10-meter NED data for the study area to create a new elevation surface that included the bathymetric data. The merged DEM was then resampled to match the 550 x 550 ft computational grid used by MODFLOW.
A new bathymetric elevation raster based on lake cells from the computational grid was extracted from the resampled DEM. The lake DEM raster cells were converted to polygons using ArcGIS conversion tools, and the elevation values were edited so that water could flow through the reservoir to end at a single cell near the reservoir outlet. Additional edits were made to ensure that no sinks existed in the data. Sinks form if the elevation of a single cell is lower than the four surrounding orthogonal cells. Flow to diagonal cells was not considered. At various stages of editing, stage-capacity volumes were calculated for the lake DEM using the ArcGIS 3D Analyst surface volume tool (Environmental Systems Research Institute, 2013). Volumes were compared to stage-capacity volumes reported in the 2004 bathymetric survey report (Ferrari, 2005). The final estimated bathymetric DEM volumes were within 1,000 acre-ft (0.3 percent of Lahontan Reservoir’s capacity) of the Reclamation reported peak storage at the spillway. Figure 16 shows the final stage-capacity to surface-area relationship used in MODFLOW based on the modified to the numerical grid described previously.
Because Lahontan Reservoir is managed, releases from the reservoir were specified as the sum of the measured flow at the USGS Carson River below Lahontan Reservoir near Fallon, NV streamgage (10312150) about 1-mile downstream from the dam (figs. 2B, 15B) and flow of water diverted to Rock Dam Ditch (Allison Danner, Bureau of Reclamation, written commun., October 2014) at a diversion point between the dam and the USGS gage (fig. 15B).
Specified-head boundaries were applied along the northeast, eastern, and southeastern borders of Churchill Valley, where little or no data were available (fig. 15). Based on the interpretation of groundwater elevations measured in early 2009 from wells in and outside the model domain, groundwater flows in along the southern boundary (Maurer, 2011, fig. 14D), and flows in and out along the north side of Churchill Valley. Substantially higher groundwater elevations were observed in Ramsey Canyon (fig. 2B), sustaining gradients likely to contribute groundwater flow to the model domain. Conversely, groundwater-level measurements from wells north of Lahontan Reservoir showed that hydraulic heads decrease to the north, indicating groundwater movement away from Lahontan reservoir. Groundwater-level observations from in and around the Dead Camel Mountains indicated that groundwater flows east through the mountains. Specified-head cells were based on the interpreted groundwater elevation determined by Maurer (2011); additional field data could allow refinement of the current set of specified heads. Table 8 summarizes the boundary conditions packages invoked for Eagle Valley, Dayton Valley (Carson Plains subbasin), and Churchill Valley.
Table 8.Summary of active MODFLOW (Harbaugh, 2005) packages in the three simulated valleys.
[CHD, specified heads; ET, evapotranspiration; GW, groundwater; LAK, lake cells; MNW, multi-node well; RCH, mountain front recharge; SW, surface water; UZF, unsaturated-zone flow; —, not active; X, active]
The MCR model was calibrated using the automated parameter-estimation software (PEST; Doherty, 2010a, b). PEST attempts to minimize the differences (residuals) between simulated (modeled) and observed (measured or estimated) values of water levels and fluxes through automated adjustment of selected model parameters. To accomplish that goal, PEST uses an objective function that quantifies differences between simulated and observed target values. Applying a gradient search technique, PEST determines the set of “best” parameter values that result in a reasonable fit between simulated and observed quantities. For the MCR model, target observations include (1) water-level observations from 202 wells; (2) river flows and the gains and losses between successive streamgages in the study area; (3) stage in Lahontan reservoir; and (4) estimates of actual ET using satellite data. The use of a diverse set of calibration targets facilitated an expanded parameterization of the flow model and the use of a broader set of parameters for model calibration than would otherwise be possible (Hunt and others, 2006). Final estimated parameter values were then evaluated to ensure they were hydrologically reasonable. Table 9 provides a list of calibrated MODFLOW parameters, including the package they are in.
Table 9.List of MODFLOW parameters calibrated by parameter-estimation software (PEST; Doherty, 2010a, c) in the middle Carson River (MCR) model.
[ft/d, feet per day; ft, feet; ft2/d, square feet per day; Kh, horizontal hydraulic conductivity; Kv, vertical hydraulic conductivity; Ss, Specific storage; Sy, Specific yield; —, no data; acre-ft/yr, acre-feet per year]
An important goal for the MCR model is to assess the effect of alternative management strategies (for example, transferred water-rights, reclaiming wastewater effluent, and increased pumping from existing wells) on flow in the Carson River in the study area. To help ensure model-predicted aquifer responses are realistic, the calibration period included a wide range of hydrologic conditions. By calibrating to this wide range of conditions (drought, average, and flood years), the baseline simulation is better for evaluating the effects of alternative water-management scenarios on the hydrologic system. Similar to the wide range of ambient hydrologic conditions during the simulation period, the historical record of anthropogenic forcings, for example river diversions and pumping, during the modeled period also spanned a wide range of water usage.
Pilot points are arbitrarily selected points that facilitate estimation of spatially distributed hydraulic properties of an aquifer. Because cell-by-cell estimation of aquifer properties is not computationally possible, pilot points offer a compromise between strict, piecewise constant-zonal (that is, “zonation”) approaches and underdetermined cell-by-cell estimation of spatially distributed aquifer properties (Doherty, 2003; Doherty and others, 2011). The flexibility afforded by pilot points allows parameter heterogeneity to emerge during automated parameter-estimation routines in areas where observations support it and, at the same time, keeps the number of adjustable parameters within a reasonable range. As the parameter values assigned at pilot-points are perturbed, the associated spatially continuous parameter field is re-Kriged (Deutch and Journel, 1998) and used by the process model, in this case, MODFLOW-NWT. As pilot-point values are adjusted, a lower overall objective-function value results (the objective function value is the simulated-to-measurement difference), and heterogeneity in the parameter field is established.
Pilot points are spaced 3,300 ft from each other in layers 2 through 5 of the numerical model grid. To help reduce the total number of adjustable parameters, the estimated hydraulic conductivity and specific yield arrays in layer 3 are the same as in layer 2, meaning that both layers use the same set of pilot points. Layers 4 and 5 also use the same set of pilot points and therefore have equivalent hydraulic conductivity fields. Even using this approach, more than 5,000 adjustable parameters were used in the MCR PEST simulation. Given the number of adjustable parameters, the lengthy model-run times (approximately 2–3 hours resulting from 138,424 active cells and 574 stress periods) and that the gradient search methodology in PEST requires one forward model run for each adjustable parameter, a number of steps were followed to help keep the parameter-estimation problem tractable: (1) initial model runs only included a steady-state stress period, where steady-state observations were the average of the available transient observations; (2) only one time step per stress period in the MODFLOW model was allowed; (3) a select set of parameter values was tied to other parameters so that their values changed together, thereby reducing the total number of adjustable parameters; and (4) some of the adjustable parameter values were set to specific values to further reduce the total number of adjustable parameters.
In addition to reducing the number of adjustable parameters, the MCR calibration procedure used the regularized inversion technologies available in PEST, including Tikhonov regularization and single-value decomposition, to help stabilize the parameter-estimation process (Doherty, 2003).
Regularization helps to bring numerical stability to the inverse problem PEST is designed to solve (Doherty, 2003). Moreover, regularization allows the modeler to apply expert knowledge (sometimes referred to as “soft” knowledge) to the parameter estimation problem. That is, regularization provides a user-specified “fallback” value deemed reasonable for each parameter (or pilot point) if observations are lacking (Morway and others, 2023b). In regions of the model where historical observations provide enough information to override user-specified preferred values (that is, vastly improved model fits resulted from adjusting parameter values away from the regularized, or “preferred,” values), PEST introduces parameter heterogeneity supported by the data. Without regularization, PEST can over-fit the MODFLOW model, such that parameters take on widely varying and unreasonable values for only minor improvements to model fit.
In the MCR model, regularization also was used to enforce preferred homogeneity throughout the model domain. With preferred homogeneity, individual pilot points only deviate from their neighbors if a great enough improvement in model fit to the observation data is achieved. In this way, the introduction of unnecessary heterogeneity (“over-fitting”) to the spatially distributed parameter arrays (for example, hydraulic conductivity) is minimized.
The objective function minimized by PEST is the sum of squared weighted residuals, where residuals are calculated as the simulated value minus the observed value. There is no limit on the number of observations or observation types that can be incorporated in PEST’s objective function. Because the relative contribution of each residual to the overall objective function value depends on the weight assigned to each observation, however, observation weights must be chosen carefully. In addition, the selection of appropriate observation weights can limit the effect of uncertain observations and enables comparison of measurements with non-commensurate units in a single objective function because weighted residuals are dimensionless.
Observations for the MCR model were divided into groups by observation types (that is, groundwater elevations, differences between streamflow at two successive river gages, and so on). The objective function was balanced by adjusting observation weights such that the relative contribution of the sum of the squared weighted residuals was roughly equal among all groups. This prevents any one observation type from dominating the parameter-estimation process at the expense of the others, which can happen when an observation group has many more observations than other groups. Individual observation weights within groups were identical.
Time-Series Processing Approach
Several surface-water flow time-series datasets were used during model calibration. In this report, time-series were weekly, the length of stress periods and time steps used in the model, or in some instances observed and modeled equivalents were aggregated in monthly time series for subsequent comparison. With appropriate processing, the information content contained in surface-water flow time-series data can be distilled to useful observation targets from otherwise noisy time-series data (Walker and others, 2009).
To this end, eight alternative time series based on surface-water streamflow data from several streamgages were prepared using the surface-water model utility time series processor (TSPROC; Doherty, 2008; Westenbroek and others, 2012). The datasets as shown in table 10 are (1) weekly average flow rates (TSPROC datasets 1–3; Morway and others, 2023b); (2) weekly average Lahontan Reservoir stage (TSPROC dataset 4); (3) monthly mean flow rates (TSPROC dataset 5–7); (4) mean monthly flow rates (TSPROC dataset 8–10); (5) weekly change in average flow rates or, in other words, the difference in flow rates at a given location between successive model time steps (TSPROC datasets 11–13); (6) change in weekly average stage of Lahontan Reservoir (TSPROC dataset 14); (7) cumulative annual volume of water passing a USGS gage (TSPROC datasets 15–17); and finally, (8) cumulative groundwater gains and losses in surface-water flows between any two gages in the study area for the simulation period (TSPROC datasets 18 and 19). The TSPROC utility also was used to post-processes model results to be equivalent to the processed observations.
Table 10.Description of the observation group names appearing in the automated parameter-estimation (PEST) control file (Morway and others, 2023b).
[TSPROC, time series processor]
Assessment of Baseline Model Calibration
The calibration performance measures used to establish the reliability of the baseline model before applying the model to various alternative management scenarios are described in this section. The term “baseline” model refers to the calibrated transient model that simulates historical conditions spanning water-years 2001 to 2010. Comparisons between observations of streamflow, groundwater elevations, reservoir stage, and ET and the respective simulated equivalents are described.
The effect of pumping on gains and losses in the MCR is an important aspect of the baseline simulation; the variation in river gains and losses resulting from the investigated alternative management strategies is a key model prediction that area stakeholders are keenly interested in. Hence, to verify that the baseline model adequately simulated historical river flows, observed, and simulated streamflows are presented for three of the five gages in the model region. Two of the five streamgages, Carson River near Carson City, NV (10311000) and Carson River below Lahontan Dam near Fallon, NV (10312150), were used to specify inflow to the modeled region and releases from Lahontan Reservoir, respectively, and were therefore not used for calibration. Owing to the wide range over which streamflow varies in the Carson River system, a direct comparison between simulated and observed streamflow could potentially mask model bias. For example, model misfit of low flows could appear inconsequential relative to model misfit of high flows because of the relative magnitude of these two groups of residuals. Thus, to investigate whether bias was present, absolute, and cumulative streamflow were compared for observed and simulated river flow, which is presented here from upstream to downstream streamgages on the Carson River (at Deer Run Road near Carson City, 10311400; at Dayton 10311700; and near Fort Churchill, 10312000).
The Carson River at Deer Run Road streamgage was operational during the entire MCR simulated period. Figure 17A shows no meaningful differences between the simulated and gaged streamflow record, except for small differences during high flows. Simulated results compared well to the gaged streamflow during low- and mid-flow conditions, but projections of cumulative annual streamflow were slightly greater than that of gaged streamflow by about 0.7 percent, 15,900 acre-ft, of the cumulative streamflow of 2,354,000 acre-ft at the Deer Run Road streamgage by the end of the 11-year simulation, about 1,450 acre-ft annually (fig. 17B). To reduce this bias, fits to the other targeted observation groups, for example the simulated Lahontan stage values, might have to be downgraded (poorer fit). Thus, the slight overestimation of roughly 2 ft3/s for the simulation period was considered reasonable.
The Carson River at Dayton streamgage operated from October 1, 2002 (the start of water year 2003), through the end of the simulation period; hence, the period of observation constitutes roughly 75 percent of the simulation period. As with the Deer Run Road streamgage, simulated low- and mid-flow periods matched their respective gaged streamflows well, except during high flows, a difference likely caused by comparing daily observed values to average weekly simulated streamflows (fig. 18A). Overall, bias in simulated values at the Dayton streamgage was minimal, amounting to 0.86-percent (16,500 acre-ft) overestimation of streamflow during the period for which gaged streamflow records were available, or 2,060 acre-ft/yr (fig. 18B). The bias in simulated values at the Dayton streamgage was not independent of the bias at the Deer Run Road streamgage. That is, if groundwater–surface-water (GW–SW) upstream gains and losses were simulated accurately between the Deer Run Road and Dayton streamgages, the degree of bias for the Dayton streamgage would be expected to be very similar to that of the Deer Run Road streamgage. Because the bias for the Dayton streamgage was less than that for the Deer Run Road streamgage, some amount of compensatory GW–SW gains and losses likely took place between the two streamgages that counteracted the bias in simulated streamflow estimates between the Carson River near Carson City and Deer Run Road streamgages.
The Carson River at Fort Churchill streamgage was active during the entire simulated period. As with the previously discussed streamgages, low- and mid-flows were well simulated by the model (fig. 19A). Similar to the Dayton streamgage, the bias in simulated values at the Fort Churchill streamgage was minimal at about 37,300 acre-ft of cumulative streamflow during the 11-year simulated period, which was 1.6 percent of the total cumulative observed streamflow (fig. 19B), or about 3,390 acre-ft/yr. Less noticeable in figure 19B is that the variability of the year-by-year bias for some years was negative (simulated cumulative streamflow was less than observed cumulative streamflow). For example, in 2000, simulated cumulative streamflow at Fort Churchill was less than the observed streamflow by 7.6 percent or 14,570 acre-ft, and in 2002, simulated cumulative streamflow was 13.1 percent greater than observed streamflow, or 16,820 acre-ft more. This variability in the year-by-year bias is a result of the calibration and parameter-estimation process in which the objective function seeks the minimum error variance. Overall, the simulated mean error in simulated flows for the Carson River at Fort Churchill streamgage was 4.5 cubic feet per second, which was roughly 1.6 percent of the mean flow during the simulation period. The Nash-Sutcliffe efficiency (NSE; Krause and others, 2005), a unitless goodness-of-fit metric commonly used in the hydrologic sciences to evaluate model performance, was 0.93 at the Fort Churchill streamgage. The maximum NSE value possible of 1.0 indicates an exact match.
Many previous seepage studies along the MCR relied on differencing values recorded at successive USGS streamgages (Maurer, 1997; Maurer and others, 2008). Differencing of the gaged streamflows at the Carson River near Carson City and Fort Churchill streamgages, for example, showed a steady depletion of river flow (loss of river flow to the groundwater system) between the two sites (fig. 20), with the most notable exceptions in years 2000 and 2006. The slight gain in streamflow between these two sites in early 2000 reflects a year of below-average groundwater discharge (fig. 9) that followed a 5-year period of higher-than-average discharge. Under this type of streamflow-aquifer condition, gains are often caused by groundwater seepage to the river channel following wetter than average years (Maurer, 2011). In contrast, the total annual streamflow from the Carson River in 2006 was well above average (nearly twice the average annual flow, fig. 9) and generally followed a period of below-average streamflow indicating that streamflow to the river from the ephemeral tributaries along the MCR during this wet year was contributing flow in excess of the loss to groundwater.
Simulated and observed net loss of streamflow between the Carson River near Carson City and Carson River near Fort Churchill streamgages was similar despite small differences in the magnitude of annual gains or losses. Factors contributing to the differences between simulated and observed streamflow include uncertain diversion rates among the 10 diversions between these 2 gages, ungaged inflows from tributaries that were assumed to be dry during the entire simulated period, and inaccurate simulation of GW–SW gains and losses. Overall, the long-term simulated trend agreed well with cumulative observed streamflow, although variability was greater in simulated annual streamflow.
Simulated Groundwater Elevations
Observations of groundwater elevations were important to guiding the calibration of aquifer properties, including the effect of hydraulic conductivity, aquifer storage, and recharge distribution. Overall, 5,296 groundwater-level observations from 202 distinct wells in Eagle, Dayton, and Churchill Valleys were used to calibrate the model, results from which are summarized in table 11. In the parameter-estimation routine, every monitoring well was assigned to one of three groups, depending on the valley in which the well is located. In each of these groups, weights assigned to wells were uniformly adjusted up or down such that the contribution to the overall objective by each group was roughly equal.
Table 11.Groundwater monitoring wells used in calibration of the middle Carson River model (Morway and others, 2023b).
[Horizontal coordinates registered to North American Datum of 1983; —, no data; Elevation datum: NGVD 29, National Geodetic Vertical Datum of 1929; NAVD 88, North American Vertical Datum of 1988. Abbreviations: USGS, U.S. Geological Survey; ID, identification; mm/dd/yyyy, month/day/year]
|USGS site identifier||Well ID
|Pumping well ID||Latitude
|Land-surface elevation datum||Elevation accuracy
|Well construction date
|Number of water-level measurements, construction date through
Sept. 30, 2010
|Average water-level altitude, construction date through Sept. 30, 2010
|Number of water level measurements, Oct. 1, 2000, through Sept. 30, 2010||Average water-level altitude, Oct. 1, 2000, through
Sept. 30, 2010