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.

Geographic Setting

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).

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 (ft³/s), in 2007 (15th smallest year on record out of 81 years) and 550,600 acre-ft/yr (760 ft³/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).

Hydrogeologic Setting

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.

Newton Solver

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).

Unsaturated-Zone Flow

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.

Hydrogeologic Units

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 (mi², or 249,400 acres), 278.1 mi² (178,000 acres) of which is unconsolidated basin-fill sediment. The remaining area, 111.6 mi² (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).

Hydrologic Properties

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

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).

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 [2012]; 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 (km²). 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 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/).

Landsat scene identifier	Sensor	Date (mm/dd/yyyy)
LT50430332000176XXX02	Landsat 5 TM	06/24/2000
LT50430332000208XXX02	Landsat 5 TM	07/26/2000
LT50430332000240XXX02	Landsat 5 TM	08/27/2000
LT50430332001162XXX04	Landsat 5 TM	06/11/2001
LT50430332001194LGS01	Landsat 5 TM	07/13/2001
LT50430332001226LGS02	Landsat 5 TM	08/14/2001
LT50430332002165LGS01	Landsat 5 TM	06/14/2002
LT50430332002197LGS01	Landsat 5 TM	07/16/2002
LT50430332002213LGS01	Landsat 5 TM	08/01/2002
LT50430332003152LGS01	Landsat 5 TM	06/01/2003
LT50430332003184LGS03	Landsat 5 TM	07/03/2003
LT50430332003216PAC02	Landsat 5 TM	08/04/2003
LT50430332004171PAC02	Landsat 5 TM	06/19/2004
LT50430332004203PAC01	Landsat 5 TM	07/21/2004
LT50430332004219PAC01	Landsat 5 TM	08/06/2004
LT50430332005173PAC01	Landsat 5 TM	06/22/2005
LT50430332005205PAC01	Landsat 5 TM	07/24/2005
LT50430332005237PAC01	Landsat 5 TM	08/25/2005
LT50430332006160PAC01	Landsat 5 TM	06/09/2006
LT50430332006192PAC01	Landsat 5 TM	07/11/2006
LT50430332006224PAC01	Landsat 5 TM	08/12/2006
LT50430332007163PAC01	Landsat 5 TM	06/12/2007
LT50430332007195PAC02	Landsat 5 TM	07/14/2007
LT50430332007227PAC01	Landsat 5 TM	08/15/2007
LT50430332008166PAC01	Landsat 5 TM	06/14/2008
LT50430332008198PAC01	Landsat 5 TM	07/16/2008
LT50430332008230PAC02	Landsat 5 TM	08/17/2008
LT50430332009200PAC01	Landsat 5 TM	07/19/2009
LT50430332009232PAC01	Landsat 5 TM	08/20/2009
LT50430332010171PAC01	Landsat 5 TM	06/20/2010
LT50430332010203EDC00	Landsat 5 TM	07/22/2010
LT50430332010235EDC00	Landsat 5 TM	08/23/2010

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)]

Table 7.    Summary of land classification total areas (in acres) for each year of the simulation.

Land classification	Year
	2000	2001	2002	2003	2004	2005	2006	2007	2008	2009	2010
Eagle Valley hydrographic area
Mexican Ditch
Irrigated	439.7	447.7	487.3	438.0	463.7	463.7	455.8	453.7	491.2	401.9	471.4
Not irrigated	57.5	49.4	9.9	59.1	33.4	33.4	41.3	43.4	5.9	95.2	25.7
Permanently dried	31.7	31.7	31.7	31.7	31.7	31.7	31.7	31.7	31.7	31.7	31.7
Urbanized	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Riparian transition	180.7	180.7	180.7	180.7	180.7	180.7	180.7	180.7	180.7	180.7	180.7
Dayton Valley hydrographic area: Carson Plains subbasin
Rose Ditch1
Irrigated	134.8	134.8	134.8	118.9	0.0	11.1	22.9	0.0	0.0	0.0	0.0
Not irrigated	0.0	0.0	0.0	15.9	46.5	35.4	23.7	46.5	46.5	46.5	134.8
Permanently dried	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Urbanized	0.0	0.0	0.0	0.0	88.3	88.3	88.3	88.3	88.3	88.3	0.0
Riparian transition	8.1	8.1	8.1	8.1	8.1	8.1	8.1	8.1	8.1	8.1	8.1
Fish Ditch
Irrigated	104.7	104.7	104.7	104.7	104.7	104.7	104.7	104.7	104.7	104.7	104.7
Not irrigated	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Permanently dried	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Urbanized	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Riparian transition	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Baroni Ditch
Irrigated	117.7	117.7	117.7	117.7	92.5	99.9	75.3	75.3	84.8	75.3	84.8
Not irrigated	0.0	0.0	0.0	0.0	25.2	17.8	9.4	9.4	0.0	9.4	0.0
Permanently dried	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Urbanized	0.0	0.0	0.0	0.0	0.0	0.0	32.9	32.9	32.9	32.9	32.9
Riparian transition	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Cardelli (Rocky Point) Ditch
Irrigated	628.4	589.8	558.0	499.1	386.1	355.9	382.4	138.3	170.9	250.8	305.0
Not irrigated	133.6	152.7	184.5	165.1	223.3	253.5	227.0	437.5	404.9	325.0	457.0
Permanently dried	19.2	19.2	19.2	19.2	19.2	19.2	19.2	19.2	19.2	19.2	19.2
Urbanized	0.0	19.5	19.5	97.8	152.6	152.6	152.6	186.3	186.3	186.3	0.0
Riparian transition	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Quilici Ditch
Irrigated	219.5	219.5	219.5	219.5	219.5	219.5	219.5	219.5	219.5	219.5	219.5
Not irrigated	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Permanently dried	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Urbanized	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Riparian transition	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Gee Ditch
Irrigated	202.9	202.9	202.9	202.9	197.6	204.0	209.4	197.6	197.6	197.6	202.9
Not irrigated	6.5	6.5	6.5	6.5	11.9	5.4	0.0	11.9	11.9	11.9	6.5
Permanently dried	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Urbanized	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Riparian transition	14.7	14.7	14.7	14.7	14.7	14.7	14.7	14.7	14.7	14.7	14.7
Carson Plains subbasin totals
Irrigated	1,408.0	1,369.4	1,337.6	1,262.7	1,000.3	995.1	1,014.2	735.3	777.4	847.8	916.9
Not irrigated	140.1	159.2	191.0	187.5	306.9	312.1	260.1	505.3	463.2	392.8	598.3
Permanently dried	19.2	19.2	19.2	19.2	19.2	19.2	19.2	19.2	19.2	19.2	19.2
Urbanized	0.0	19.5	19.5	97.8	240.9	240.9	273.8	307.4	307.4	307.4	32.9
Riparian transition	22.8	22.8	22.8	22.8	22.8	22.8	22.8	22.8	22.8	22.8	22.8
Dayton Valley hydrographic area: Bull Canyon subbasin
Koch Ditch2
Irrigated	194.0	176.5	75.7	200.2	75.8	218.4	218.4	194.0	218.4	166.3	218.4
Not irrigated	24.4	41.9	142.7	18.2	142.6	0.0	0.0	24.4	0.0	52.1	0.0
Permanently dried	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Urbanized	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Riparian transition	9.9	9.9	9.9	9.9	9.9	9.9	9.9	9.9	9.9	9.9	9.9
Houghman-Howard Ditch
Irrigated	231.7	231.7	231.7	231.7	231.7	249.0	237.0	254.2	254.2	231.8	254.2
Not irrigated	22.5	22.5	22.5	22.5	22.5	5.3	17.3	0.0	0.0	22.5	0.0
Permanently dried	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Urbanized	151.8	151.8	151.8	151.8	151.8	151.8	151.8	151.8	151.8	151.8	151.8
Riparian transition	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Bull Canyon subbasin totals
Irrigated	425.7	408.2	307.4	431.9	307.5	467.4	455.4	448.3	472.6	398.1	472.6
Not irrigated	46.9	64.4	165.2	40.7	165.1	5.3	17.3	24.4	0.0	74.5	0.0
Permanently dried	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Urbanized	151.8	151.8	151.8	151.8	151.8	151.8	151.8	151.8	151.8	151.8	151.8
Riparian transition	9.9	9.9	9.9	9.9	9.9	9.9	9.9	9.9	9.9	9.9	9.9
Churchill Valley hydrographic area
Buckland Ditch
Irrigated	555.7	532.4	532.4	404.0	306.0	469.2	555.7	452.4	555.7	555.7	532.4
Not irrigated	17.2	23.3	23.3	151.7	249.7	86.4	0.0	103.3	0.0	0.0	23.3
Permanently dried	14.7	31.9	31.9	31.9	31.9	31.9	31.9	31.9	31.9	31.9	31.9
Urbanized	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Riparian transition	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
Middle Carson River study area
Study reach totals
Irrigated	2,829.0	2,757.7	2,664.7	2,536.6	2,077.6	2,395.4	2,481.1	2,089.7	2,296.9	2,203.5	2,393.4
Not irrigated	261.66	296.28	389.31	438.99	755.05	437.2	318.66	676.33	469.17	562.54	647.21
Permanently dried	65.58	82.8	82.8	82.8	82.8	82.8	82.8	82.8	82.8	82.8	82.8
Urbanized	151.75	171.25	171.25	249.59	392.62	392.62	425.5	459.17	459.17	459.17	184.63
Riparian transition	213.4	213.4	213.4	213.4	213.4	213.4	213.4	213.4	213.4	213.4	213.4

¹

The Rose Ditch was abandoned in 2004.

²

Alternately named Chaves/Koch Ditch.

Evapotranspiration

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).

Irrigation Applications

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).

Streams

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):

$$C=\frac{K_{str}LW}{l}$$

(1)

where

K_str

is the hydraulic conductivity of the streambed material (length per time),

L

is the length of the stream reach,

W

is the width of the stream reach, and

l

is the thickness of the streambed.

Initial values of K_str 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.

Lahontan Reservoir

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).

Pilot Points

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

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.

Observations Weighting

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.

Simulated Streamflow

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 ft³/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.

The difference between observed and simulated hydraulic head, termed the groundwater-level residual, is summarized by valley in figure 21. Residuals are calculated by subtracting the observed values from the simulated values; hence, a negative residual indicates a simulated groundwater level that is below the observed level, and a positive residual indicates a higher simulated than observed groundwater level. A count of simulated groundwater elevations that were within ±5 ft of their corresponding observation in Eagle, Dayton, and Churchill Valleys is 543, 2,363, and 208, respectively, out of 1,022, 3,459, and 810 total observations for each valley. This represents 53, 68, and 26 percent of the residuals, respectively, for each valley. Total counts of simulated groundwater elevations within ±10 ft are 840, 2,698, and 497 in Eagle, Dayton, and Churchill Valleys, respectively. Simulated heads in the MCR model are considered accurate and a reasonable representation of observed conditions because (1) the observed groundwater elevations were generally well-distributed over the modeled regions (fig. 22); (2) most of the groundwater-level residuals were less than 10 ft (fig. 21); and (3) the overall bias in simulated values, as indicated by the mean error, was less than +1 ft in Eagle and Churchill Valleys and roughly +3 ft in Dayton Valley (fig. 21). In each valley, residuals of groundwater elevations were generally least near the center of the valley where topographic relief is lowest. Conversely, monitoring wells near the higher relief boundaries of the model, where steeper hydraulic gradients persisted, generally had larger average residuals (fig. 22). Considering all 5,296 simulated and observed groundwater level pairs (i.e., residuals), the mean error was 1.42 ft; the mean absolute error was 7.71 ft, and the percent bias was −0.1 percent.

In addition to the evaluation of groundwater-level residuals, simulated trends in groundwater elevation also were assessed for each valley. Areas identified as showing long-term drawdown, corresponding to areas highlighted in figure 14A of Maurer (2011), are drawn in each respective valley in figure 22 and are hereinafter referred to as “recognized drawdown areas.” The average simulated groundwater level was calculated for these regions for the model period, and it is discussed in the sections that follow. In addition to the average groundwater level in the recognized drawdown areas, simulated and observed groundwater-level hydrographs were compared for select wells in Eagle, Dayton, and Churchill Valleys to help assess the accuracy of model results.

Eagle Valley

Groundwater elevation hydrographs for observation wells W-5, W-8, W-9, and W-14 (fig. 23) in Eagle Valley, in the recognized drawdown areas, were compared to simulated groundwater elevations. Although the model appeared to capture overall groundwater elevation trends in Eagle Valley, simulated water levels for wells W-9 and W-14, near the northern extent of recognized drawdown area (fig. 22A), did not capture the annual variability in the observed record that was likely due to annual pumping patterns. Thus, model results should be evaluated on a multi-year basis and not used to evaluate within-year variability.

Dayton Valley

In Dayton Valley, observation wells W-48, W-50, W-52, and W-54 are near an area of long-term groundwater elevation decline (fig. 22B) of about 10 ft from 1995 to 2009 (fig. 24). The simulated groundwater elevations in three of the wells (W-50, W-52, and W-54) were higher than observed groundwater elevations by more than 15 ft (fig. 24B–D). In contrast, the simulated groundwater elevation for well 48 was lower than the observed groundwater elevation by more than 20 ft (fig. 24A). Although the groundwater-elevation hydrograph for well W-50 indicated long-term drawdown in the vicinity of the well, wells W-52 and W-54, and even well W-50, indicated a slow rebound in the groundwater elevation in this area starting in 2008 that coincided with the start of the economic recession and slowdown in home building (Coleman and others, 2008). Thus, for this location, the effectiveness of the parameter-estimation procedure was limited by the conflicting water-table residuals (underestimation and overestimation) among neighboring observation wells. Relatively large and variable residuals are due, in part, to simulated groundwater elevations that were relatively constant compared to the greater within-year variability of observed groundwater elevations in wells W-48, W-50, and W-54. The observed groundwater elevation drawdown was approximately 20 ft in well W-50 (fig. 24B) between 2000 and 2008, with considerable within-year variability. No observations were available from 2000 to 2008 to use to deduce groundwater elevation trends in wells W-52 and W-54. Of the three wells—W-50, W-52, and W-54—well W-50 is the closest to a large supply well. Model discretization coupled with a gridded pilot-point spacing (approximately 1 pilot point per 36 cells) led to a final calibration for this part of the model that was not able to fully capture all the localized effects of pumping.

Churchill Valley

Owing to the large area of western Churchill Valley (fig. 22C) and Lahontan Reservoir (fig. 22D), 13 observation wells distributed throughout the region were selected to represent and illustrate the model performance—4 wells (W-92, W-93, W-98, and W-103), for which hydrographs are shown in figures 25A–D, are in recognized drawdown areas west of Churchill Butte and Misfit Flats (fig. 22C); 4 wells (W-107, W-112, W-122, and W-124), for which hydrographs are shown in figures 26A–D, are in or near a recognized drawdown area west to northwest of Churchill Butte (fig. 22C); and 5 observation wells (W-182, W-185, W-188, W-190, and W-192), for which hydrographs are shown in figures 27A–E, surround Lahontan Reservoir (fig. 22D). The simulated groundwater elevations generally agreed with the four observation-well hydrographs shown in figures 25A–D, where the largest (less than 4 ft) difference in simulated and observed groundwater elevation was in wells W-98 and W-103 (figs. 25C, 25D, respectively). Simulated groundwater elevations in all four observation wells shown in figure 25 are lower at the end of the 10 yr simulation period, similar to the observations.

Similarly, groundwater elevations were simulated reasonably well for the four observation wells closer to Churchill Butte, but still in the recognized drawdown area in western Churchill Valley (figs. 26A–D); differences between simulated and observed water levels ranged between approximately 3 (fig. 26A) and –4 ft (figs. 26C–D). Moreover, the downward trend evident in the observed groundwater-levels through time was reproduced by the steadily declining simulated groundwater levels for all four wells shown in figure 26.

Differences between simulated and observed groundwater elevations for the five wells around the perimeter of Lahontan Reservoir (fig. 27) ranged from approximately −5 (fig. 27B) to −32 ft (fig. 27A), although the difference was approximately 8 ft for well W-190 (fig. 27D). In general, groundwater-level residuals for nearly all the wells on the southwest side of Lahontan Reservoir (west of the Dead Camel Mountains) were within ±10 ft (fig. 22D), whereas groundwater-level residuals for wells near the northwest side of Lahontan Reservoir were nearly +30 ft (well W-159, fig. 22D), but generally were +20 ft or less (for example, wells W-153, W-155, and W-176; fig. 22D). All the wells on the north side of the model area also are near the boundary of the model, where groundwater-level gradients were steep and uncertain boundary conditions were fixed (specified by the user). A modeled groundwater level within 40 ft, and much less in some of the wells (for example, wells W-187 and W-190), was deemed adequate because groundwater-level contours generally matched those interpreted by Maurer (2011, fig. 14D). Because each of the simulated alternative management scenarios was evaluated for effects on Lahontan Reservoir storage levels (among other effects), and not the simulated heads around Lahontan Reservoir, accurate simulation of Lahontan Reservoir stage was weighted more heavily (that is, given priority) than close fits to observed groundwater elevations north and east of the reservoir.

Simulated groundwater elevations in the basalt aquifer in the Dead Camel Mountains (figs. 27A, 27C) were higher and lower than observed groundwater elevations; differences ranged from about −30 to +40 ft, which is similar to fits achieved among observation wells completed in unconsolidated basin-fill sediments. Higher residuals for wells completed in the basalt aquifer of Churchill Valley are likely due to complex and uncertain geology (fractured semi-consolidated rock), less well-known land-surface elevations in areas of varying topographic relief from which groundwater elevations are calculated, and a limited number of observations (typically, 2–10). Simulated hydraulic head residuals for wells around the northeastern part of Lahontan Reservoir reveal that the model is poorly calibrated to observed water levels in this area (fig. 22D). Some simulated water levels were lower than observed values by nearly 60 ft (well W-200, fig. 22D), whereas other simulated groundwater elevations were higher than observed values by more than 100 ft (well W-193, fig. 22D). This indicates steeper than simulated groundwater-level gradients, possibly the result of complex geology owing to the presence of the volcanic Dead Camel Mountains. Moreover, the proximity of the model boundary—represented with specified heads—to a hydrologically dominant feature like Lahontan Reservoir makes it difficult to accurately simulate groundwater elevations in this region of the model. Nevertheless, efforts were made during the calibration process to adjust model parameters (for example, lake conductance) to reduce head residuals in this area; however, smaller head residuals could not be achieved without causing greater error in the water budget for Lahontan Reservoir and the surrounding groundwater system. For the purpose of alternative management scenario evaluation, it is more important for the model to accurately simulate Lahontan Reservoir stage than groundwater elevations in this region of the model. That is, groundwater elevations near the reservoir are not evaluated under alternative management scenarios, but reservoir stage and storage are.

Simulated Evapotranspiration

Remotely sensed ET, as calculated by Mapping EvapoTranspiration at high Resolution with Internal Calibration (METRIC; Allen and others, 2007), was compared to MODFLOW-simulated values, an increasingly common practice to constrain the overall water budget better for the simulated region (Morway and others, 2013; Carroll and others, 2015; Doble and Crosbie, 2017). When comparing MODFLOW-simulated values with METRIC-derived values, it is important to keep in mind that METRIC values are themselves “modeled” estimates, although they provide better estimates of spatially and temporally varying ET than do regionally constant estimates. As with observations of groundwater elevations or streamflow, matching simulated and observed ET helps constrain an otherwise uncertain component of the overall water budget and, as a result, can improve confidence in the simulation. Although METRIC results are themselves a modeled output, the model-predicted values are constrained by land-based climate stations combined with Landsat imagery. Lastly, with models like METRIC or RESET (remote sensing of evapotranspiration; Elhaddad and Garcia, 2011), greater accuracy has been achieved estimating ET for irrigated lands than for dry-land areas (Allen and others, 2013); therefore, simulated ET was compared with METRIC results only for irrigated areas, where METRIC performs best.

To improve MODFLOW-simulated ET estimations, adjustments were made to a global multiplier (across space and time) applied to the extinction-depth term in the UZF input file. Extinction depth is the depth at which groundwater ET ceases (Harbaugh, 2005; Niswonger and others, 2006). Through adjustment of the extinction depth, over- and underestimated ET could be removed from the model until a good fit with METRIC results was achieved. Because of the hydrologic connection with other components of the water budget, adjustments that effectively increased or decreased estimates of ET could lead to a poorer fit among other simulated factors, such as total streamflow.

Model-simulated monthly ET rates compared well to METRIC-derived (“observed”) rates of ET (Huntington and others, 2018). For irrigated acreage supplied by each of the selected ditches in Eagle, Dayton, and Churchill Valleys, simulated ET for the growing season (April–October) was less than observed ET, on average, by –10, –12, and –7 percent, respectively (fig. 28). Pursuit of improved simulation of ET for better model fits was compromised by two main factors: (1) model grid-cell boundaries did not correspond with field boundaries, where field boundaries define strong differences in ET rates between irrigated lands and surrounding desert lands; and (2) the “observations” were themselves modeled estimates and therefore had an associated error of their own. The total simulated volume of ET for the 11 water years of the period of analysis was 305,900 acre-ft (approximately 27,805 acre-ft/yr), 18 percent of which was from the unsaturated zone and supplied through irrigation water, and the remaining 82 percent originated from groundwater.

Simulated Lahontan Reservoir Stage, Inflows, and Losses

In general, the calibrated MODFLOW model does an excellent job of simulating the rise and fall of the reservoir stage as volumes are captured and released before and during each irrigation season, (fig. 29A). Two of the three major inflows and outflows from the reservoir are specified. Inflow from the Truckee Canal is specified based on the streamflow from the Truckee Canal near Hazen, NV streamgage (10351400) measured near the outlet of the canal. Reservoir releases are specified according to streamflow from the Carson River below Lahontan Reservoir near Fallon, NV (10312150) measured downstream from the dam. The third major water budget component, inflow from the Carson River, is the primary source of flow entering the reservoir and is a function primarily of the specified inflow entering the model at the Carson River at Carson City streamgage (10311000) modified somewhat by approximately 65 miles of simulated GW–SW gains and losses affected by irrigation diversions, irrigation return flows, groundwater pumping, and tributary inflows from the perennial streams in Eagle Valley (fig. 11). As with the groundwater-level residuals, weekly stage residuals were determined by subtracting the observed Lahontan Reservoir stage from the simulated reservoir stage. In general, a long-term (multi-year) oscillation and an annual oscillation in the residuals are evident (fig. 29B). At the end of each irrigation season, the residual trend changed markedly—through 2006 and resuming in 2008, the upward and downward trends in the residuals reverse coincident with the start and end of the irrigation season, respectively. In other words, stage residuals grow increasingly positive during the irrigation season, and this trend tends to reverse at the end of the irrigation season. This behavior indicates that when simulated releases from Lahontan Reservoir to the Newlands project ceased, model fit with respect to simulated reservoir stage improved, and that the within-year stage residuals were affected by uncertain, yet specified, diversions both upstream from and directly from Lahontan Reservoir. When irrigation activity stopped each year and the only fluxes to and from the river upstream from Lahontan Reservoir were groundwater and surface-water exchange, model performance improved. This result might also be partially explained by the cessation of riparian and phreatophyte ET along the river corridor and around the delta region where the Carson River enters the reservoir. With less ET being simulated, specified inflow at the upstream end of the model is transmitted to the reservoir with one less sink groundwater evapotranspiration (GWET) to satisfy. Finally, errors related to the specified monthly reservoir evaporation rate likely contribute to the within-year variability of stage residuals. Thus, on an annual basis, less simulated activity (riparian ET, diversions, tail-water returns, among others) in the model leads to reduced residuals.

The long-term variation in simulated stage generally ranged within 6 ft of the observed stage, with one exception. Early in water year (WY) 2010, the reservoir reached the lowest stage on record up to that point in the simulation. During this time, simulated reservoir stage was consistently higher than observed stage (more than 6 ft), but the stage residual was reduced as WY 2010 progressed and returned to no more than ±5 ft of the observed values toward the end of the irrigation season. Owing to how little water is in the reservoir at low stage, the poorly simulated stage had minimal effect on the overall reservoir water budget. Overall, the mean reservoir-stage residual was 2.25 ft, and the root mean square error was 3.54 ft.

Lahontan Reservoir Groundwater Surface-Water Interaction

There are three primary sinks from Lahontan Reservoir: (1) scheduled releases from the reservoir; (2) groundwater–surface-water interaction resulting in a net seepage loss from Lahontan Reservoir to the underlying groundwater system; and (3) evaporation losses. Whereas scheduled, or managed, releases are considered a beneficial use of water by the State of Nevada, seepage and evaporation losses are considered non-beneficial uses of water.

Groundwater–surface-water interactions with Lahontan Reservoir were evaluated for three geographic areas: (1) west of the semi-consolidated rock unit (Dead Camel Mountains); (2) overlying the semi-consolidated rock unit; and (3) northeast of the semi-consolidated rock unit (fig. 15B). In the western and eastern parts of Lahontan Reservoir, groundwater and surface water interact in the unconsolidated basin-fill sediments west and east of the Dead Camel Mountains. Along the central part of the reservoir, water is lost to (or gained from) the groundwater system through basalt rock in the Dead Camel Mountains.

The part of the reservoir west of the Dead Camel Mountains is the only area in the reservoir with appreciable inflow from groundwater (fig. 30A). Simulated monthly groundwater inflow reached a maximum of approximately 1,000 acre-ft in 2000, 2001, and 2007–09 and a maximum closer to 500 acre-ft between 2002 and 2006 and again in 2010. Simulated groundwater discharge to the reservoir generally starts in mid- to late-summer and extends through the winter when the reservoir stage typically reaches the lowest levels of the year and bank storage is draining back to the reservoir. Conversely, when the reservoir stage is high, typically during the spring and early summer months, gradients are reversed and simulated seepage from the reservoir recharges the adjacent groundwater system. The month in which maximum seepage losses were greatest was in the spring of 2005 (fig. 30A), when more than 2,000 acre-ft was lost to groundwater. In the middle (fig. 30B) and eastern (fig. 30C) parts of Lahontan Reservoir, simulated groundwater–surface-water interaction was almost entirely reservoir seepage losses. The losses ranged from zero in the second half of 2008 to approximately 1,000 acre-ft/month in various years. Only during the fall months of 2007, 2008, and 2009, when simulated stages were lowest, did groundwater enter the central part of Lahontan Reservoir in the canyon cutting through the Dead Camel Mountains (fig. 30B) and only in small (less than 200 acre-ft/month) amounts.

Monthly total groundwater inflow for the reservoir was roughly the same as for the west side of the Dead Camel Mountains (fig. 30D), less than 1,000 acre-ft for most years. The monthly total seepage loss from the reservoir reached a maximum of approximately 4,000 acre-ft in the second half of water year 2005 when the reservoir stage was at its second highest stage (June 2006 reached a higher stage). The monthly net monthly difference between groundwater inflow and seepage loss indicated that during low stages, there was a large contribution to the reservoir storage from groundwater (fig. 30D).

Seepage losses in Lahontan Reservoir were greatest during periods of high stage (near or at maximum storage capacity); for example, simulated seepage loss rates were high near pumping wells in the northwest part of the reservoir during the week of June 11, 2005, near pumping wells and in the northeast part of the reservoir where the reservoir stage is much higher than surrounding groundwater (fig. 31A). At high stage, nearly the entire reservoir area contributes to seepage. Integrating over the lake model-grid cells, the total seepage loss for the week depicted in fig. 30A (June 11, 2005) was 996 acre-ft (or 72 ft³/s), which was approximately 0.4 percent of the total storage in the reservoir. This seepage rate is equivalent to a flux of 0.005 ft³/s per wetted lakebed acre. Much of this water is gained back in fall and winter when the reservoir is at low stage and the bank storage drains back to the reservoir. Note that groundwater discharge to the reservoir was only considered for cells that were inundated; the LAK package does not accumulate discharge for cells above the lake stage.

Simulated groundwater inflow to the reservoir during the 7 days starting on September 20, 2008, a period of low stage, was approximately 18.6 ft³/s. A simultaneous seepage loss of roughly 0.6 ft³/s was simulated near the dam, making for a net accrual of groundwater in the reservoir roughly equal to 18.0 ft³/s. As a demonstration of the model’s ability to simulate Lahontan Reservoir inflows and outflows, a water budget based on modeled values was completed for the week of September 20, 2008, and compared to a water budget using gaged values. During this period, the average gaged release from the reservoir was 7 ft³/s, and the combined inflow from the Carson River near Fort Churchill, NV and Truckee Canal near Hazen, NV streamgages averaged 121 ft³/s. Accounting for the observed reservoir storage increase during this period required more inflow than could be explained by the difference between surface inflow and reservoir release, and the difference was matched by the net groundwater discharge to the reservoir. For the observed increase in storage, an additional 22.8 ft³/s of inflow was needed, of which 18.6 ft³/s was simulated by the model, a reasonable difference of 4.2 ft³/s. Table 13 summarizes average annual water-budget components for Lahontan Reservoir.

Table 13.    

Simulated and observed water budget components for Lahontan Reservoir, Nevada, 2000–10.

[ft³/s, cubic foot per second; NA, not available; WY, water year; —, no data]

Table 13.    Simulated and observed water budget components for Lahontan Reservoir, Nevada, 2000–10.

Site	Simulated (ft3/s)	Observed (ft3/s)
Average annual inflow, WY 2001–10
Carson River1	200,150	200,300
Truckee Canal	124,800	124,800
Precipitation	2,670	NA
Groundwater	4,210	NA
—	331,830	—
Average annual outflow, WY 2001–10
Lahontan Dam managed release	284,000	279,600
Rock Dam Ditch	3—	4,200
Evaporation	23,900	NA
Groundwater	15,380	NA
Change in storage2	280	2,660
—	323,560	—

¹

Additional gains or losses between the reservoir and the Carson River near Fort Churchill streamgage may be measured at the gage.

²

Estimated increase in storage (difference) between October 1, 2000, and September 30, 2010, was 72 acre-feet.

³

Simulated releases into the Carson River include the measured release at the Carson River below Lahontan Reservoir stream gage in addition to the Rock Dam Ditch release, which diverts water upstream of the Carson River below Lahontan Reservoir stream gage.

Cumulative Losses from Lahontan Reservoir

Lahontan Reservoir, like any reservoir in an arid environment, contributes to water losses from the Carson River system through evaporation and seepage. Total losses from the reservoir during the simulated period were roughly 400,000 acre-ft (fig. 32). Figure 32A shows that a little over 30 percent of the total loss from the reservoir was caused by seepage from the reservoir, or about 125,000 acre-ft. Annually, this equates to approximately 11,000 acre-ft/yr, although figure 32A shows that cumulative losses due to seepage slowed in 2007. Consistently lower average stage in the reservoir and a reduction in pumping in the Silver Springs area after 2007 help explain the slowdown in cumulative losses from the reservoir. Figure 32B presents another view of the seepage where weekly simulated seepage amounts are normalized by the wetted-lakebed area for that week. Viewed in this way, the per acre seepage loss rate never exceeded 0.04 ft³/s. Only when the reservoir stage was very low in late 2008 and 2009 did groundwater inflow exceed 0.05 ft³/s/acre.

The slowdown in reservoir losses was also evident in the cumulative evaporation losses that equaled approximately 275,000 acre-ft for the period of the simulation, or 23,900 acre-ft/yr (fig. 32A). Considering the large surface area covered by the reservoir, nearly 16,000 acres at full capacity, and that the annual rate for total evaporation, about 4.25 ft/yr, was within the normal range for reservoirs in northern Nevada (Huntington and McEvoy, 2011), the estimated cumulative evaporative losses for Lahontan Reservoir were reasonable.

Alternative Management Scenario 1: 40–40–20

As municipalities search for additional water to bolster supplies for meeting the needs of a growing population (fig. 3), transfer of water currently used for agricultural production to municipal use is one possibility for meeting that need. This scenario explores the impact of a potential management scenario that includes the exercise of a surface-water right(s) sourced from an “induction well(s).” In Nevada, induction wells are located close to the river and require a 100 ft sanitary seal (Ed James, Carson Water Subconservancy District, written commun., February 2023).

When a surface-water right is sourced from groundwater, the constraints regulating the maximum amount of monthly diversion for irrigation remain intact. A description of the rules governing the timing and maximum diversion amounts are given in the Alpine Decree (U.S. District Court, 1980b). Under the Alpine Decree, maximum surface-water diversions may not exceed 40 percent of the total water right during April, 40 percent of the total water right during May, and the remaining 20 percent of the total water right is limited to June (fig. 33). Irrigators could choose to irrigate less than the maximum amount each month and, in so doing, irrigate beyond these 3 months (without exceeding the total water rights). For example, irrigators could choose to irrigate with 20 percent of their total allotment in April, again in May, and yet again in June, leaving 40 percent for the coming months (assuming the river does not go into regulation, in which case irrigators only get their water when in priority). In the scenario simulated for this strategy, however, the maximum monthly limits for April, May, and June were exercised by the induction wells for each year. Monthly allotments of 40, 40, and 20 percent in April, May, and June, respectively, gave this scenario its “40–40–20” name.

The general intent of the 40–40–20 management scenario was to decrease municipal pumping withdrawals in April, May, and June in areas of long-term groundwater-level drawdown in Eagle and Dayton Valleys (figs. 22A–B), as well as provide a reliable source of water for new residential development. Thus, reductions in pumping from municipal supply wells were offset with increased pumping from the induction wells using the transferred surface-water rights. At the beginning of July, when the water rights transferred under this scenario are exhausted, simulated pumping from existing municipal supply wells are increased to levels above their baseline pumping amount to supply not only the original demand, but the additional demand associated with the newly permitted residential development based on the transferred water rights. Because the transfer of water rights from surface water to groundwater for augmented induction well pumping alters the timing and amount of return flows from the previously flood-irrigated fields (the prior use of the transferred water right), model results from this scenario were compared to those of the calibrated baseline model to assess downstream effects that arise from the change in water use.

Selected ditches for which the transfer of water rights was used in the 40–40–20 scenario included a partial reduction of irrigated acreage for the Mexican Ditch in Eagle Valley and a complete reduction of irrigated acreage for Fish Ditch in Dayton Valley. Acreages associated with the water rights for Mexican Ditch (328 acres) and Fish Ditch (133 acres) diversions were based on the Federal Water Master records (table 14).

The total volume of water available for transfer to the induction wells was calculated as the original acreage associated with the transferred water right multiplied by 2.5 ft (table 14). The 2.5 ft multiplier is from the Alpine Decree and is generally less than the consumptive use of the irrigated crop. Hence, more water would be left in the river in this scenario in average to wet years because, historically, much more than 2.5 ft was diverted from the river when water was available. Furthermore, inefficiencies associated with flood irrigation, including seepage loss along the delivery ditch, recharge from over-irrigation, and tail-water runoff from the end of the field, would remain in the river in this scenario. Although these inefficiencies were returned to the river in the baseline simulation (minus any additional riparian ET losses associated with a higher water table), the temporal shift in river flows was a key difference of this scenario in that irrigation return flows were diminished, decreasing later season flows, and springtime flows were higher because water was not removed for irrigation. Hence, the 40–40–20 scenario resulted in a net savings of water, meaning more water remained in the river than was withdrawn by the induction wells. In Eagle Valley, the 40–40–20 scenario reduced the irrigated land supplied by Mexican Ditch by 328 acres. Multiplying this area by the 2.5 ft duty brought the total volume associated with the transferred water right for Eagle Valley to 820 acre-ft. Thus, under 40–40–20 scenario, 40 percent of 820 acre-ft, or 328 acre-ft was pumped from the Eagle Valley induction well in April and May, leaving 164 acre-ft (20 percent) to be pumped in June. Similarly, in Dayton Valley, 133 acres irrigated from Fish Ditch were dried by transferring the associated water-rights to the Dayton Valley induction well; An example of the estimated monthly pumping amounts for Dayton and Eagle Valleys are given in table 15.

Simulation of the 40–40–20 scenario was accomplished by modifying model-input files defining pumping rates and ditch diversion amounts. During April, May, and June, pumping from municipal wells in and around the recognized drawdown areas were reduced in an amount corresponding to the increased induction well pumping that captured water associated with the transferred water rights. Conversely, during the remaining months of the year (July through March), pumping from the municipal wells increased above baseline pumping to supply the expanded residential development. Moreover, pumping from the induction wells dropped to its baseline level from July through March. Because the water associated with the transferred water rights recovered by the induction wells was available during the 3-month 40–40–20 period (April–June), the way monthly pumping rates were determined for the 40–40–20 alternative management scenario is described next.

The monthly estimates of water usage for each new home were estimated, and the total number of new homes associated with each transferred water right (Mexican Ditch in Eagle Valley and Fish Ditch in Dayton Valley) was determined. Using these values, the amount of the transferred water right consumed by new development was calculated. The remaining volume of water available for established homes during the 40–40–20 simulation period was the amount that can be relaxed from the municipal supply wells in the recognized drawdown areas during the same period. To make these calculations, it was assumed that the total annual water usage by a new home was 0.6 acre-ft/yr (Ed James, Carson Water Subconservancy District, personal commun., May, 2014), and that roughly half the amount (0.3 acre-ft/yr) was for indoor purposes whereas the other half was for water requirements for landscaping between April and September (table 15). Further, the number of new lots supported in Eagle and Dayton Valleys, respectively was calculated as follows:

Utilizing the monthly breakdowns of transferred water-right usage by new development (table 15), the amount of water available compared to historical pumpage in the recognized drawdown areas was calculated in table 16. In summary, more than 22 and 9 million ft³ were used for existing homes. Water remaining after satisfying the total new home usage was 17,286 and 4,356 ft³ (that is, the “leftover” after supplying the new residential development) in Eagle and Dayton Valleys, respectively, during the 40–40–20 simulation period.

A summary of the simulated monthly increases in induction well pumping during the 40–40–20 period to capture the full allotment of water transferred from Mexican and Fish Ditches is provided in table 17. During the non-40–40–20 period, pumping from the municipal supply wells (wells other than the induction wells) needed to be increased to provide water to the established and new homes after the full duty associated with the transferred water rights was met. Note that the simulated totals for induction wells in the last row of table 17 verified the values in table 16. Because Eagle and Dayton Valleys have only one induction well each, all of the water transferred from Mexican and Fish Ditches was simulated as pumped from the induction well in the respective valley.

Monthly adjustments to existing municipal supply wells were such that the sum of pumping reductions from multiple wells equaled the values shown in table 18. Note that the total reduction in pumping from the production wells was nearly equal to the “water available to satisfy other homes during 40–40–20” in table 16.

During dry years, when water availability is limited, the full duty may not be realized by induction wells after the river goes into regulation during irrigation season and available surface-water supplies are depleted. The MCR model, however, does not distribute the available water supply in priority; instead, water use (diversions) is specified before running the model. As a result, the pumping increases associated with the transferred water right are exercised to their full allotment, even during the dry years of the simulation period. Because of this, the MCR model estimates resulting from the transferred water rights within the MCR system are conservative (that is, likely maximums) with respect to effects on groundwater; the actual pumping would likely be somewhat less than the simulated amount during dry years when the river system is in regulation.

An important limitation of the 40–40–20 scenario is that under the Alpine Decree, unused amounts of a transferred water right remain in the river and become available to the next junior water-right holder, a process the current model is not equipped to simulate. Therefore, the reduction of agricultural lands along the middle Carson River is likely to result in additional water being delivered to Lahontan Reservoir rather than the next junior water right. Future model modifications, using an integrated hydrologic and river-operations code like MODSIM-MODFLOW (Morway and others, 2016), could facilitate simulation of water-right priorities with the ability to route augmented river flows after a water-right transfer to the next in-priority user. In average and above-average water years, the effect of this limitation is expected to be minimal owing to the fulfillment of junior water rights under these conditions.

Effects on Regional Water Levels and Stream and Canal Seepage

Alternative management scenario results showed a range of aquifer responses across the three simulated valleys. Results of each scenario for Eagle, Dayton, and Churchill Valleys were evaluated for changes from baseline conditions in groundwater elevations, Carson River streamflow, and Lahontan Reservoir storage. For all three comparisons, differences between baseline conditions and each respective scenario were calculated using the model stress periods between January 2000 and December 2010.

In Eagle Valley, the MCR model results indicate the average groundwater-level changes in the recognized drawdown areas were roughly the same as for the baseline (model calibration) scenario by the end of the 11-year simulation period (fig. 34A). In general, the MCR model results indicate that adopting the 40–40–20 management scenario did not reverse the relatively substantial groundwater-level decline of nearly 30 ft in Eagle Valley (fig. 34A). Differencing the average groundwater elevations between the 40–40–20 and baseline scenarios through time shows a maximum of roughly 2 ft less drawdown than the baseline simulation in 2006 (fig. 34B). Along with the long-term declines in average groundwater elevations, seasonal head varied by 1 to 2 ft during within-year periods relative to the baseline simulation (fig. 34B).

The extent of the roughly 2 ft of higher water levels under 40–40–20 management was evident on the west side of Eagle Valley, roughly half of which is inside the recognized drawdown area. Higher water levels of 0.1 to 0.5 ft outside the recognized drawdown area were indicated for much of Eagle Valley under 40–40–20 management. A small pocket of drawdown in northwestern Eagle Valley resulted from additional water being pumped from a municipal well during the non-irrigation season (fig. 35). No effect to the water table was observed in the vicinity of the induction well, indicating that the Carson River was the main source of water for induction-well pumping rather than groundwater storage.

In addition to effects on the water table, model results were post-processed to determine the effect of the 40–40–20 scenario on GW–SW gains and losses. Simulated differences in seepage between the river and shallow alluvial aquifer near the induction wells increased on the east side of Eagle Valley (fig. 36). The largest difference in river seepage was in the cell containing the northern induction well, which increased approximately 0.5 ft³/s, or about 1 acre-ft/day on average, over baseline conditions during the growing season. Note that in figure 36, positive values indicate a net decrease in river flow either from increased losses from the river to the shallow alluvial aquifer or diminished groundwater return flow to the river.

Simulated average groundwater elevations in the recognized drawdown area in Dayton Valley (fig. 22B) for the 40–40–20 management scenario were roughly equivalent to those in the baseline scenario (fig. 37A). There were brief periods of slightly higher average groundwater elevations; however, differences typically were less than a foot (fig. 37B). That is, sustained long-term drawdown in Dayton Valley was not apparent, neither under the baseline, nor the 40–40–20 scenario.

Although changes in groundwater elevations from baseline conditions were minimal near the induction well in central Dayton Valley (fig. 38), there was a pronounced drawdown near Fish Ditch, where irrigation diversions were reduced, in southwestern Dayton Valley. Without canal seepage loss and recharge from flood irrigation to sustain groundwater elevations in this area, water levels declined up to about 5 ft along some areas next to Fish Ditch.

Simulated surface-water (seepage) losses to the shallow alluvial aquifer increased in Dayton Valley because of the 40–40–20 scenario, or alternatively, groundwater return flows to the river diminished. Regardless of whether seepage losses from the river or diminished groundwater return flow was the cause, the change of roughly 0.01–0.03 ft³/s per model cell, for a total of approximately 0.7 ft³/s as simulated at the Fort Churchill streamgage location (located just downgradient of figure 39), was minor.

The largest change in seepage rates along the Dayton Valley reach of the Carson River was in the cell containing the induction well. In the baseline scenario, groundwater discharge to the river in this model cell was about 0.1 ft³/s. In the 40–40–20 scenario, there was a slight seepage loss from the river to the alluvial aquifer in the same model cell, indicating a slight reversal from a seepage gain to a seepage loss condition; however, the change from gain to loss was small, and when rounded to the nearest whole ft³/s, resulted in no net gain or loss from the river during the irrigation season. Thus, a net reduction of 0.1 ft³/s in simulated river flow through this reach of river was the result of a reduction in groundwater discharge to the river and was within the uncertainty associated with the model as well.

A time series of the average simulated groundwater elevations in the recognized drawdown area in Churchill Valley (fig. 22C) shows no response to the 40–40–20 scenario (fig. 40). A lack of groundwater-level response in Churchill Valley was expected because there was no physical implementation of the 40–40–20 scenario in this valley, and the recognized drawdown area in Churchill Valley was spatially removed (more than 4 mi) from simulated changes in upstream river conditions.

Effects on Flow at Fort Churchill and Storage in Lahontan Reservoir

Total annual flow at the Fort Churchill streamgage site because of 40–40–20 management was less than the baseline simulation for all but one of the simulated years (figs. 41A–B). The largest difference in flow was for 2001 at roughly three-quarters of a percent, which equates to about 800 acre-ft/yr (table 19). For all other simulated years, the total annual flow at the Fort Churchill streamgage generally was less than the baseline simulation by an average of 0.43 percent (table 19), except for 2006, which showed a slight positive change in the total annual flow. At the end of the 11-year simulation period, Lahontan Reservoir stage was 0.35 feet lower, corresponding to a reduction in storage of 1,810 acre-ft relative to the baseline simulation. The total annual flow reduction at the Carson River near Fort Churchill streamgage (fig. 41) was not expected to exactly match the overall reduction in Lahontan Reservoir storage (fig. 42) owing to changes in GW–SW gains and losses in the reservoir as estimated by simulation of the alternative management scenarios. The differences in figure 42 were calculated as the management scenario time series minus the baseline time series.

The original focus of the 40–40–20 scenario was to explore how changes in water management affect other water users of the river. The scenario did not explicitly consider land-use change along the river that may accompany the 40–40–20 water-management practices. Owing to the unexpected reduced flows simulated at the Carson River near Fort Churchill streamgage, a more detailed analysis of the water budget was performed for the fallowed lands. The results revealed that although surface-water delivery to the fallowed land served by Mexican Ditch had ceased, the model continued simulating ET from the shallow groundwater at nearly the same rate as when the fields were being irrigated. Because of the proximity of the fallowed land to the Carson River, as well as its location downgradient from the Mexican Ditch, which still conveys some irrigation water in the scenario, the water table remained shallow and, as a result, sustained ET directly from the water table. With the total ET from the fallowed lands remaining roughly the same and the associated increase in induction well pumping simulated by this scenario, the discharge at the Carson River near Fort Churchill, NV streamgage was less for most years in the 40–40–20 scenario relative to the baseline scenario (fig. 41B).

To further explore the effects land fallowing may have on the simulated discharge at the Carson River near Fort Churchill, NV streamgage, a variation of the original 40–40–20 scenario was created. Under this 40–40–20 variation scenario, the potential ET (pET) rate for the fallowed lands was set to half of the pET rate assigned to the same lands under the original 40–40–20 scenario. The pET represents the maximum ET rate that can be achieved under well-watered conditions and is generally greater than actual ET during water-stressed conditions. The revised pET rate, although arbitrarily chosen, represents a reduced ET rate that may accompany a land-use change away from irrigated pasture. An example land-use change could include the conversion to a land-cover with less consumptive water use or to a partial removal of phreatophytic vegetation. Note that a 50-percent reduction in pET may not correspond to a 50-percent reduction in actual simulated ET because of the non-linear controls affecting crop consumptive use in the model.

With pET of the fallowed lands cut in half, the simulated ET (that is, consumptive use) was reduced enough to result in an overall annual increase of flow at the Carson River near Fort Churchill streamgage (fig. 41C). Moreover, the augmented flows entering Lahontan Reservoir slightly increased the storage of Lahontan Reservoir relative to the baseline simulation (fig. 42). Thus, the efficacy of the 40–40–20 scenario depended on the land-use changes that accompanied its implementation. Without a corresponding reduction of ET from the fallowed irrigated lands next to the river, increased induction well pumping in April, May, and June negatively affected total annual flows. If a land-use change that accompanies the 40–40–20 scenario successfully reduced ET from the retired lands, however, model results indicated it is possible to offset the effects of increased induction well pumping simulated by this scenario.

Table 19.    

Change in total annual flow at Fort Churchill, Nevada, expressed as a percentage, under the 40–40–20 and 40–40–20 with reduced potential ET scenarios, 2000–10.

[pET, potential evapotranspiration]

Table 19.    Change in total annual flow at Fort Churchill, Nevada, expressed as a percentage, under the 40–40–20 and 40–40–20 with reduced potential ET scenarios, 2000–10.

Year	Change (percent)
	40–40–20	40–40–20 with reduced pET
2001	–0.77	–0.35
2002	–0.27	0.25
2003	–0.30	0.06
2004	–0.37	0.11
2005	–0.11	0.14
2006	0.02	0.18
2007	–0.73	–0.17
2008	–0.75	–0.27
2009	–0.30	0.10
2010	–0.24	0.11
Average	1–0.43	0.02

¹

Average does not include 2006.

Alternative Management Scenario 2: Reclaimed Effluent for Aquifer Storage and Recovery

Reclaiming treated-wastewater effluent offers water managers an additional and potentially reliable source of water to support the region’s projected population and home growth. Lyon County Utilities (in Dayton Valley) produced about 600,000 gallons per day (gal/d) of treated effluent in 2017 (Mike Workman, Lyon County Utilities, written commun., 2017). From April through October every year, this water is contracted for delivery to a local golf course. During the rest of the year (November through March), this water is sent to rapid infiltration basins (RIBs) and currently is considered “un-reclaimed” water. Given the relative proximity of the RIBs to the Carson River, treated effluent sent to the RIBs recharges the underlying alluvial aquifer and eventually flows back to the river as groundwater discharge (fig. 43).

The general intent of alternative management scenario 2 was to use induction wells to reclaim treated effluent discharged to RIBs in Dayton Valley and then returned to the Carson River as groundwater return flow. Alternative management scenario 2 preserved wastewater-effluent sent to the RIBs in the baseline simulation but increased Dayton Valley induction-well pumping by an amount equivalent to the RIB infiltration. Moreover, an equal volume of induction-well pumping was subtracted, or “relaxed,” from that of municipal wells in the recognized drawdown area south of the town of Dayton. Unlike alternative management scenario 1, induction-well pumping and relaxed municipal pumping were simulated for the winter and spring months (November–March).

Alternative management scenario 2 was simulated in two parts to evaluate the different rates of treated-wastewater effluent sent to the RIBs and the associated municipal and induction well pumping adjustments during the winter and spring period. In scenario 2A, generation of wastewater effluent continued at the baseline rate of 600,000 gal/d; however, the baseline induction-well pumping rate was augmented by an additional 600,000 gal/d from November through March (fig. 43). Wastewater-effluent was simulated as specified recharge in the cells containing the RIBs from November through March and was specified as 0 gal/d during the summer period (April through October) when the treated wastewater effluent was contracted for delivery to a local golf course. For April through October, when induction-well pumping was restored to its baseline level, municipal-well pumping also returned to its baseline levels. For November through March, simulated municipal-well pumping rates in the recognized drawdown area (southwest region of Dayton Valley) were reduced by a total of 600,000 gal/d in accord with the enhanced induction-well pumping.

Alternative management scenario 2B simulated a “potential maximum build-out capacity” of wastewater effluent to the Dayton Valley RIBs (Mike Workman, Lyon County Utilities, written commun., 2017; (fig. 44). Under this scenario, total specified recharge rates in cells containing the RIBs (wastewater effluent) were increased to 1,500,000 gal/d more than the baseline simulation (900,000 gal/d more than scenario 2A) during the November through March period and increased to 900,000 gal/d (300,000 gal/d more than scenario 2A) during the summer period (April through October) when 600,000 of the 1,500,000 gal/d were under contract for delivery to the local golf course. A wastewater-effluent discharge rate of 1,500,000 gal/d represented a 2.5-fold increase in the wastewater-effluent rate simulated in scenario 2A. To achieve such a large increase in wastewater-effluent, alternative management scenario 2B assumed an accompanying 2.5-fold increase in municipal-well pumping relative to the baseline simulation. An increase in groundwater pumping of this magnitude was justified because of existing groundwater permits issued by the Nevada State Engineer that had not yet been exercised. Thus, municipal- and induction-well pumping amounts varied depending on the season and are summarized in table 20. As in alternative management scenario 2A, the wastewater-effluent in scenario 2B was assumed to return to the Carson River by groundwater flow that originated in the RIBs. Under scenario 2B, pumping from the Dayton Valley municipal supply wells (excluding the induction well) increased by 220 percent compared to the baseline scenario for the months of April through October. Between November 1 and March 31, when water demand slackens, average pumping among the Dayton Valley municipal supply wells (excluding the induction well) was reduced by 8 percent relative to the baseline simulation because the induction well—recovering water released to the RIBs—was able to supply the total water demand.