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.