History of Water Use

The basin has a long history of water use. Native Americans utilized surface-water resources in the basin for thousands of years prior to the arrival of non-Native settlers (Risser, 1982). Diversions of surface water for irrigation increased following the development of settlements in the area near the communities of Bluewater and Grants during the 1870s (Risser, 1982). Settlers began diverting springflow from Ojo del Gallo around 1870 (fig. 1), which had formerly discharged into a swampland and flowed into the Rio San Jose about 3 mi downstream from the community of Grants, to irrigate 1,000–2,000 acres of land in the area near the community of San Rafael (Risser, 1982; Frenzel, 1992). Ojo del Gallo springflow was utilized for irrigation of farms and gardens until the flow ceased in 1953 (Gordon, 1961). Prior to the construction of Bluewater Dam, settlers diverted surface water from Bluewater Creek for irrigation. Flows were regulated on Bluewater Creek from 1894 to 1909 when a series of earth dams were built near the current location of Bluewater Dam. Following the construction of Bluewater Dam in 1927, four main canals were constructed to deliver water diverted from Bluewater Creek to over 10,000 acres of land within the Bluewater-Toltec Irrigation District (Gordon, 1961). However, runoff to Bluewater Lake proved to be inconsistent from year to year and never large enough to support the intended acreage. By 1948 the amount of land within the Bluewater-Toltec Irrigation District had been permanently reduced to 5,488 acres (Gordon, 1961).

The inconsistency of interannual surface-water supply resulted in the introduction of large-capacity irrigation wells in 1944 and the rapid expansion of the number of wells during the 1940s (Risser, 1982). Most of these wells obtained water from the Permian San Andres Limestone and Glorieta Sandstone aquifer system (SAGA) (Gordon, 1961). On the Acoma and Laguna Pueblos, irrigation wells obtain water primarily from the Quaternary alluvium and alluvial-basalt sequences aquifer (Risser and Lyford, 1983).

Following the discovery of uranium in the basin in the early 1950s, the economy and resultant groundwater use began to shift from predominantly agricultural to industrial (Gordon, 1961; Risser, 1982). From 1951 to 1980, the Grants uranium district yielded more uranium than any other district in the United States (McLemore, 2011). Most of the uranium production was from the Westwater Canyon Member of Morrison Formation (McLemore, 2011), which at many mines required the dewatering of the ore-bearing and overlying units to access the ore. Numerous mills were constructed to process uranium, which also required large-capacity wells that pumped groundwater from the SAGA (Gordon, 1961). Untreated mine dewatering and mill-process waters were discharged to natural waterways and were a significant source of contamination of sediments, alluvial aquifers, and deeper aquifers in areas of faulting (Schoeppner, 2008). A variety of factors led to the cessation of mining activities in the 1980s (McLemore, 2011). Many springs ceased to flow during the mining period, and partial recovery of some springs was observed after mining operations ceased (Lorenzo and Watchempino, 2003).

Municipal development and water use expanded in the 1950s to meet the demands of the mining industry (Gordon, 1961; Frenzel, 1992). The population of the community of Grants grew from 2,251 in 1950 to 10,274 in 1960 (U.S. Census Bureau, 1960). The municipal water supply for the community of Grants was primarily obtained from the Quaternary alluvium and alluvial-basalt sequences aquifer prior to 1958 and from the SAGA afterward (Gordon, 1961; Cooper and West, 1967). Smaller communities around the community of Grants typically obtain municipal water supply from the SAGA and the Quaternary alluvium and alluvial-basalt sequences aquifer (Gordon, 1961; Cooper and West, 1967). Municipal water supply on the Acoma and Laguna Pueblos is obtained from a variety of aquifers, including the Quaternary alluvium, Dakota Formation, Zuni Sandstone, Entrada Sandstone, and sandstones in the Chinle Formation (Risser and Lyford, 1983).

Water for domestic and livestock use in the basin is obtained from springs and wells (Gordon, 1961). These domestic and livestock springs and wells obtain water from nearly all of the aquifers in the basin, but most wells obtain water from the SAGA or the Quaternary alluvium and alluvial-basalt sequences aquifer (Gordon, 1961). Many springs utilized for domestic and livestock use issue from the Cenozoic basalts aquifer on mesas near the community of Grants and Mount Taylor and in the Zuni uplift (Gordon, 1961; Risser and Lyford, 1983).

Model Boundary Delineation and Horizontal Discretization

The active PRMS area conforms to the basin, encompassing a total area of about 3,600 mi² (fig. 1) that was discretized into 500- by 500-meter (m) grid cells, each representing an HRU. The land surface elevation of each HRU was derived from a National Elevation Dataset 10-m horizontal resolution DEM (USGS, 2017) at the center of each grid cell. A total of 36,774 HRUs were active in PRMS, with 13 HRUs set as lake type to simulate Bluewater Lake, 2 HRUs set as swale type to simulate closed surface-water subbasins, and the remaining HRUs set as land type. The two swale HRUs were defined as the lowest DEM-derived land surface elevation within two closed subbasins (subbasins 14 and 15 in fig. 2) identified from National Hydrography Dataset (NHD) Hydrologic Unit Code 12 (USGS, 2019a) in the southwestern part of the study area.

Stream Network Delineation and Streamflow Routing

The HRU land surface elevations were used by the GSFLOW-Arcpy toolkit to develop the stream network, which includes the main-stem Rio San Jose and its major tributaries, and delineate subbasins in the active PRMS area (fig. 2). The level of detail of the stream network (density of tributaries to the Rio San Jose) was adjusted with the user defined values for flow accumulation and flow length thresholds until a level of detail that represented the Rio San Jose and its major tributaries was obtained. The stream network was modified to better match the location of NHD flow lines (USGS, 2019a) and stream networks visible in aerial imagery by manually adjusting select HRU land surface elevations. Fifteen subbasins were delineated using the 12 streamgages that were used for model calibration and verification (see “Precipitation-Runoff Modeling System, Calibration Approach and Datasets” section), the 2 swale HRUs (low areas that capture flow and only contribute to groundwater flow), and 1 outlet location where the Rio San Jose exits the basin.

Local depressions in the HRU land surface elevations that were not identified as swales were filled by the GSFLOW-Arcpy toolkit using the Cascade Routing Tool (version 1.3.1) (Henson and others, 2013). These local depressions (or undeclared swales) were filled, because they become HRUs with no downslope flow path and can cause numerical issues if a swale is not physically present (Henson and others, 2013). After the adjustments to the HRU land surface elevations, the toolkit used the Cascade Routing Tool with the stream network to compute the surface and shallow subsurface lateral flow paths (called cascades) among HRUs. The Cascade Routing Tool was also used to develop the associated cascade parameters required by PRMS to route simulated flows from upslope HRUs to downslope HRUs and cascade-discharge features such as streams, swales, or lakes. This study used the srunoff_smidx module to control surface runoff processes and the muskingum_lake module to control streamflow routing (Markstrom and others, 2015; Regan and LaFontaine, 2017).

Thirteen lake-type HRUs were used to simulate Bluewater Lake. The lake-type HRUs representing Bluewater Lake were defined as the grid cells with at least 30 percent of the cell area containing the NHD waterbody feature (USGS, 2019a). Managed releases from Bluewater Lake have altered the natural hydrologic system that is simulated by PRMS. Therefore, releases from Bluewater Lake were simulated in PRMS using the muskingum_lake module, with the outflow from Bluewater Lake specified in the PRMS data file. This flow replacement method does not maintain a complete water balance because the state (calculated water-content storage) of the lake is computed but not included in the water balance (Regan and LaFontaine, 2017).

Daily releases from Bluewater Lake were estimated for the transient simulation period without gaged flows using measured daily mean streamflow obtained from the USGS National Water Information System (NWIS) (USGS, 2019b) at the Bluewater Creek below Bluewater Dam, N. Mex. streamgage (referred to hereafter as “Bluewater Creek below Bluewater Dam”; USGS site number 08341500), which is about 2 mi downstream from Bluewater Dam. Using the streamflow at this streamgage to estimate releases from Bluewater Lake relied on the assumption that streamflow gains from groundwater and streamflow losses to groundwater were negligible between Bluewater Dam and the streamgage, which was reasonable because of the proximity of the streamgage to Bluewater Dam. Measured streamflow was available at the Bluewater Creek below Bluewater Dam streamgage for only a portion (March 17, 1951, through September 29, 1963, and July 21, 1989, through February 22, 2001) of the transient simulation period of the RSJIHM. Therefore, streamflow at the Bluewater Creek below Bluewater Dam streamgage needed to be estimated during the intervening transient simulation periods of the RSJIHM not covered by the measured streamflow data. About 8 mi downstream from Bluewater Dam is the Bluewater Creek near Bluewater, N. Mex. streamgage (referred to hereafter as “Bluewater Creek near Bluewater”; USGS site number 08342000). Measured daily mean streamflow obtained from the USGS NWIS (USGS, 2019b) was available for the Bluewater Creek near Bluewater streamgage from October 1, 1960, through January 4, 1973. Because of its proximity to Bluewater Dam, the Bluewater Creek below Bluewater Dam streamgage is likely most representative of releases from Bluewater Lake, but the Bluewater Creek near Bluewater streamgage provides an additional period (September 30, 1963, through January 4, 1973) of measured streamflow that was used to estimate streamflow at the Bluewater Creek below Bluewater Dam streamgage during this period.

The Bluewater Creek below Bluewater Dam and Bluewater Creek near Bluewater streamgages had an overlapping period of streamflow measurements from October 1, 1960, through September 29, 1963. Linear regressions (fig. 3) were developed for streamflow measurements at Bluewater Creek below Bluewater Dam versus two groups of streamflow measurements (daily mean streamflow less than or equal to 6 ft³/s and greater than 6 ft³/s) at Bluewater Creek near Bluewater during the overlapping period. To check the ability of these linear regressions to predict daily mean streamflow at Bluewater Creek below Bluewater Dam, the regressions were used to estimate daily mean streamflow at Bluewater Creek below Bluewater Dam from the daily mean streamflow at Bluewater Creek near Bluewater during the overlapping period. The estimated daily mean streamflow closely matched the measured daily mean streamflow at Bluewater Creek below Bluewater Dam, with a sum of daily differences (estimated minus measured) of 5 ft³/s and a root-mean-square error (RMSE) of daily differences of about 1 ft³/s, calculated as

Therefore, the linear regressions were deemed suitable to estimate daily mean streamflow at Bluewater Creek below Bluewater Dam during the period from September 30, 1963, through January 4, 1973, when measured daily mean streamflow was only available at Bluewater Creek near Bluewater. The linear regressions had nonzero vertical axis intercepts (estimated daily mean streamflow at Bluewater Creek below Bluewater Dam), and therefore negative values of estimated daily mean streamflow at Bluewater Creek below Bluewater Dam were set to 0 ft³/s.

During periods when measured daily mean streamflow was not available at either Bluewater Creek below Bluewater Dam or Bluewater Creek near Bluewater, daily mean streamflow was estimated using the annual mean streamflows (each being the arithmetic mean of the daily mean streamflow derived from Bluewater Creek below Bluewater Dam and Bluewater Creek near Bluewater, as described previously, for each year) and annual estimates of surface-water diversions for irrigation within the Bluewater-Toltec Irrigation District (data are available on request from the NMOSE). When annual surface-water diversions were estimated to be zero, daily mean streamflow was held constant during the year at 0.21 ft³/s (computed as the arithmetic mean of the annual mean streamflows for years with no surface-water diversions). When annual surface-water diversions were estimated to be greater than zero, daily mean streamflow was specified for growing season (28 ft³/s) and nongrowing season (2.6 ft³/s) periods (computed as the arithmetic mean of the annual mean streamflows for growing and nongrowing season subdivisions of each year). The growing and nongrowing season subdivisions of each year were estimated at a monthly frequency, with the growing season corresponding to months with a daily minimum air temperature greater than 32 °F for at least two-thirds of the month and the nongrowing season corresponding to the other months that did not meet these criteria. The daily minimum air temperature was obtained from the gridded, hydrometeorological datasets described in the “Climate Distribution Parameterization” section and analyzed for the growing and nongrowing season criteria at the HRU where surface water was diverted from Bluewater Creek below Bluewater Dam for irrigation within the Bluewater-Toltec Irrigation District.

Land Surface and Soil Zone Parameterization

Publicly available geospatial datasets along with the GSFLOW-Arcpy toolkit were used to produce spatially distributed land surface and soil zone parameters needed to run PRMS. The National Land Cover Dataset 2011 (USGS, 2014) was used to derive initial values for the fraction of each HRU area that is impervious (hru_percent_imperv), which were later adjusted during calibration (see “Precipitation-Runoff Modeling System, Calibration Approach and Datasets” section). Calibration of hru_percent_imperv was necessary, because the toolkit generated hru_percent_imperv values generally contradicted the classifications in the National Land Cover Dataset. For example, areas of the National Land Cover Dataset classified as “Developed, Low Intensity” to “Developed, High Intensity” (representing urban areas) had the lowest hru_percent_imperv values while areas of the National Land Cover Dataset classified as “Herbaceous” to “Shrub/Scrub” had the highest hru_percent_imperv values. The LANDFIRE 2014 (lf_1.4.0) Existing Vegetation Cover (Wildland Fire Science, Earth Resources Observation and Science Center, USGS, 2016a) and Existing Vegetation Type (Wildland Fire Science, Earth Resources Observation and Science Center, USGS, 2016b) datasets were used to generate parameters related to vegetation. The cov_type parameter defines the dominant vegetation type for each HRU (0=bare soil, 1=grasses, 2=shrubs, 3=trees, 4=coniferous), and with each vegetation type, an associated rooting depth is assigned. In addition, the summer vegetation cover density and winter vegetation cover density (covden_sum and covden_win) and summer rain interception storage capacity and winter rain interception storage capacity (srain_intcp and wrain_intcp) are determined from the vegetation type assigned to each HRU. The summer rain interception storage capacity (srain_intcp) derived from the GSFLOW-Arcpy toolkit was later adjusted during calibration (see “Precipitation-Runoff Modeling System, Calibration Approach and Datasets” section).

Soil zone parameters were derived from soils data obtained from the Digital General Soil Map of the United States or STATSGO2 (U.S. Department of Agriculture, Natural Resources Conservation Service, 2006). A dominant soil type (soil_type parameter; 1=sand, 2=loam, 3=clay), soil depth, and parameters influencing water movement through the soil were derived from the soils data. Soil depth is defined by the greater of the soil depth values provided by the soils data or rooting depth (Gardner and others, 2018). The maximum water holding capacity of the soil is determined by the parameter soil_moist_max, which influences processes such as evapotranspiration, surface runoff, and recharge (Markstrom and others, 2015). Parameters ssr2gw_rate and slowcoef_lin route water to the groundwater reservoir or to downslope HRUs (Markstrom and others, 2015).

Climate Distribution Parameterization

The Climate-by-HRU Distribution Module, climate_hru, was used to input daily precipitation and minimum and maximum air temperature by HRU to PRMS (Markstrom and others, 2015). This study assigned daily precipitation and minimum and maximum air temperature to each HRU by extracting data from gridded, hydrometeorological datasets (Livneh and others, 2015; Thornton and others, 2016) at the centroid of each HRU. The Livneh and others (2015) dataset was gridded at 1/16-degree horizontal resolution (approximately 6 kilometers [km]) and covered the period from 1950 through 2013. The Thornton and others (2016) dataset, known as Daymet, was gridded at a 1-km horizontal resolution and covered the period from 1980 through 2019. This study used the Livneh and others (2015) dataset for the RSJIHM transient simulation period from 1950 through 2013 and the Daymet dataset for the RSJIHM transient simulation period from 2014 through September 30, 2018. This study used the ddsolrad module to control the distribution of solar radiation to each HRU, the transp_tindex module to determine whether the dominant vegetation in each HRU is actively transpiring, and the potet_jh module to calculate the amount of potential evapotranspiration in an HRU.

Snow-Covered Area Parameterization

This study assigned snow area depletion curve values (snarea_curve parameter) and the maximum threshold snowpack water equivalent below which the snow-covered-area curve is applied (snarea_thresh parameter) to each HRU by extracting data from the USGS National Hydrologic Model (Regan and others, 2018) HRU occupying the majority of each RSJIHM HRU. In total, values for 184 snow area depletion curves from the USGS National Hydrologic Model were assigned to HRUs in the RSJIHM. Related parameters that were modified to accommodate these edits include ndepl (the number of snow-depletion curves), ndeplval (the number of values in all snow-depletion curves, equal to ndepl times 11), and hru_deplcrv (the index number for the snowpack areal depletion curve associated with each HRU).

Water Use Estimation and Distribution

Surface-water use in the basin has been primarily for irrigation and has varied in quantity throughout the simulation periods of the RSJIHM (Risser, 1982). Diversions of surface water from Bluewater Creek below Bluewater Dam for irrigation within the Bluewater-Toltec Irrigation District were simulated in PRMS using the water_use_read module as an external transfer of water from the stream network outside the active PRMS area (dest_type variable equal to 4) (Regan and LaFontaine, 2017). Water applied to the land surface for irrigation from those same diversions was simulated to infiltrate below the soil zone to recharge the groundwater system (hereafter referred to as “deep percolation”) in MODFLOW-NWT (see “MODFLOW-NWT, Water Use Estimation and Distribution” section) only, because PRMS has a limited representation of the groundwater system. Because of a lack of reported measurements or observations of the occurrence of surface runoff from irrigation water within the Bluewater-Toltec Irrigation District, this runoff was assumed to be a minor component of the hydrologic budget and was not included in the RSJIHM.

Daily surface-water diversions from Bluewater Creek below Bluewater Dam for irrigation within the Bluewater-Toltec Irrigation District were estimated using the estimated daily releases from Bluewater Lake (see “Precipitation-Runoff Modeling System, Stream Network Delineation and Streamflow Routing” section) and annual estimates of surface-water diversions for irrigation within the Bluewater-Toltec Irrigation District (data are available on request from the NMOSE). Daily surface-water diversions were estimated by scaling the estimated daily releases from Bluewater Lake so that the total volume of daily surface-water diversions during the growing season of each year equaled the annual estimate of surface-water diversion. The scale factors were calculated as the annual estimate of surface-water diversion divided by the total volume of daily releases during the growing season of each year.

Calibration Approach and Datasets

Calibration of PRMS was performed using PEST++ (version 4.3.20), specifically the iterative ensemble smoother algorithm (White and others, 2020). PEST++ was run using the USGS Yeti supercomputer (U.S. Geological Survey Advanced Research Computing, 2021). Parameters adjusted during calibration were those that control the (1) solar radiation; (2) potential evapotranspiration distribution; (3) evapotranspiration and sublimation; (4) air temperature and precipitation distribution; (5) streamflow; (6) Hortonian surface runoff, infiltration, and impervious storage; (7) interception; (8) soil zone storage, interflow, gravity drainage, and Dunnian surface runoff; (9) groundwater flow; and (10) snow computations (table 1), and generally included parameters identified as sensitive by Markstrom and others (2016). Parameter dimensions in the RSJIHM are identified in the “Parameter Dimension in the Rio San Jose Integrated Hydrologic Model” column of table 1; “One” indicates parameters that are constant in space and time, “HRU” indicates parameters that can vary spatially, “Monthly” indicates parameters that can vary temporally, and “HRU and monthly” indicate parameters that can vary spatially and temporally.

The calibration of parameters in PRMS that can vary spatially was achieved by using pilot points following the guidelines described in Doherty and others (2010). Pilot points were distributed uniformly across the PRMS grid at a spacing of 12.5 km, located at the centroid of every 25th HRU in the general east-west and north-south directions, resulting in 57 pilot points within the active PRMS area (fig. 2). PEST++ adjusted parameter values at each pilot point during calibration, with the parameter values for the remaining HRUs calculated by interpolating between pilot points using ordinary kriging. The pilot point approach allowed for the calibration of spatially variable parameters while minimizing computational resources (57 adjustable values per spatially variable parameter), as opposed to adjusting values at each of the 36,774 HRUs for each spatially variable parameter.

Adjustments to parameter values during calibration were based on the ability of model simulated values to reproduce measured or estimated values for several observation groups: (1) solar radiation, (2) potential evapotranspiration, (3) actual evapotranspiration, (4) precipitation and minimum and maximum air temperature, (5) snow water equivalent, (6) snow-covered area, and (7) streamflow. Observation weights within PEST++ were applied to give relatively equal importance to each observation group within the objective function. Within observation groups, values were given the same PEST++ observation weights. The ability of PRMS to reproduce measured or estimated values is discussed in the “Model Performance” section.

Model Boundary Delineation and Discretization

Like PRMS, MODFLOW-NWT was discretized horizontally into 500- by 500-m grid cells, with each grid cell directly underlying the corresponding PRMS HRU where the active areas overlap. The MODFLOW-NWT grid consists of 322 rows, 318 columns, and 8 layers with a total of 405,344 active grid cells in the horizontal and vertical dimensions and a no-flow boundary at the vertical base of the model. MODFLOW-NWT row numbering increases in a generally north-to-south direction and column numbering increases in a generally west-to-east direction. Layer 1 of MODFLOW-NWT was specified as convertible in the Upstream Weighting package (Niswonger and others, 2011) to allow for the simulation of a water table condition in the uppermost hydrogeologic unit. All other layers of MODFLOW-NWT were simulated as confined in the Upstream Weighting package.

The active MODFLOW-NWT area conforms to PRMS in the southwestern portion of the basin (plate 1), where faults, absent units from erosional processes, and reported groundwater divides were assumed to impede groundwater flow to the southwest (no-flow boundary). However, the active MODFLOW-NWT area extends beyond PRMS in the western, northern, eastern, and southern portions of the study area, encompassing a total area of about 6,300 mi². The active MODFLOW-NWT area was extended to the west to generally correspond to the location of a groundwater divide in the SAGA as contoured by Baldwin and Anderholm (1992). The northern boundary generally corresponds to the extent of the San Andres Limestone from McLaughlin and Johnson (1987). The eastern boundary corresponds to the western boundary of the Middle Rio Grande Basin groundwater-flow model (McAda and Barroll, 2002). The southern boundary is generally perpendicular to groundwater-level contours in the SAGA from Baldwin and Anderholm (1992).

In the steady-state simulation period, groundwater flow across the eastern boundary of MODFLOW-NWT into the Middle Rio Grande Basin was specified as the estimated flow from the calibration of the Middle Rio Grande Basin groundwater-flow model (McAda and Barroll, 2002) and was simulated using the Flow and Head Boundary package (Leake and Lilly, 1997); all other horizontal boundaries of MODFLOW-NWT were simulated as no-flow. In the transient simulation period, the possibility of groundwater flow across the western, northern, and eastern boundaries was simulated using the General-Head Boundary package (Harbaugh and others, 2000), with the head on the boundary specified as constant at the simulated head from the end of the steady-state simulation period. The southwestern and southern boundaries of MODFLOW-NWT were simulated as no-flow (plate 1).

Stream Network Delineation and Streamflow Routing

Several parameters required for the SFR2 package were generated by the GSFLOW-Arcpy toolkit. The stream network within the SFR2 package is simulated as groups of grid cells that are assumed to have uniform or linearly varying characteristics, referred to as “segments,” with individual grid cells referred to as “reaches.” The toolkit generated parameters to define the model row and column number of the grid cell containing each reach, the number assigned to each segment, the downstream sequential numbering of each reach within a segment, the length of the stream channel within each reach, the elevation of the top of the streambed within each reach (set as 1 ft below the adjusted HRU land surface elevations described in the “Precipitation-Runoff Modeling System” section), and the routing of streamflow from segment to segment.

The stream slope across each reach in the SFR2 package was calculated as the difference between the elevation of the top of the streambed at the reach and the next downstream reach divided by the length of the stream channel within each reach. The thickness of the streambed in each reach was set as 1 ft. The hydraulic conductivity of the streambed in each reach was adjusted during model calibration. Stream depth within each segment was set in the SFR2 package to be calculated using Manning’s equation and assuming a 50-ft-wide rectangular channel. Manning’s roughness coefficient for each segment was assumed to be equal to 0.02 (Arcement and Schneider, 1989).

Components of the stream network outside the active PRMS area (the Rio Puerco and the Rio San Jose from its outlet from PRMS to where it exited the active MODFLOW-NWT area) also were simulated in MODFLOW-NWT using the SFR2 package (plate 1). These stream network features were not delineated by the GSFLOW-Arcpy toolkit because they are outside the active PRMS area; instead, they were manually delineated using NHD flow lines (USGS, 2019a) and aerial imagery. The Rio Puerco was simulated from the Rio Puerco above Arroyo Chico near Guadalupe, N. Mex. streamgage (referred to hereafter as “Rio Puerco above Arroyo Chico”; USGS site number 08334000) to its confluence with the Arroyo Chico and then downstream to where it exited the active MODFLOW-NWT area. For the components of the stream network outside PRMS, the length of the stream channel was estimated as the length of the NHD flow line (USGS, 2019a) within each reach. The elevation of the top of the streambed within each reach was set to linearly decrease from the upstream to downstream reach on the basis of the stream channel lengths and to be below the HRU land surface elevations derived from the National Elevation Dataset DEM described in the “Precipitation-Runoff Modeling System” section. The elevation of the top of the streambed in these reaches ranged from 1 to 216 ft below the HRU land surface elevations, because the Cascade Routing Tool and resultant adjustments to the HRU land surface elevations were not applied to the stream network outside PRMS.

Inflow to the Rio Puerco and Arroyo Chico were simulated in the SFR2 package using daily mean streamflow data obtained from the USGS NWIS (USGS, 2019b) at Rio Puerco above Arroyo Chico and the Arroyo Chico near Guadalupe, N. Mex. streamgage (referred to hereafter as “Arroyo Chico near Guadalupe”; USGS site number 08340500), respectively (plate 1). Daily mean streamflow data were available at Rio Puerco above Arroyo Chico from October 1, 1951, through September 30, 2018, and at Arroyo Chico near Guadalupe from October 1, 1943, through September 29, 1986, and from October 1, 2005, through March 29, 2018. Missing daily mean values were estimated as the mean daily streamflow statistic obtained from NWIS (USGS, 2019b), which are calculated as the arithmetic mean of the daily mean streamflow for each calendar day over the period of record. A steady-state streamflow was estimated by calculating the average of the annual mean streamflow statistic obtained from NWIS (USGS, 2019b), which was calculated as the arithmetic mean of the 365 (or 366 for leap years) individual daily mean streamflows for each year.

Bluewater Lake was simulated in MODFLOW-NWT using the Lake package (Merritt and Konikow, 2000). MODFLOW-NWT grid cells in layer 1 directly underlying the 13 PRMS lake-type HRUs were used to simulate Bluewater Lake. In the steady-state simulation period, surface-water releases from Bluewater Lake were computed by the SFR2 package on the basis of lake stage, as simulated by the Lake package, relative to the top of the streambed for the first reach in the segment (ISEG variable equal to 94) downstream from the lake, by setting the FLOW variable in ISEG 94 equal to 0. In the transient simulation period, surface-water releases from Bluewater Lake were specified as the estimated daily releases from Bluewater Lake (see “Precipitation-Runoff Modeling System, Stream Network Delineation and Streamflow Routing” section) to represent the anthropogenic alteration to the natural hydrologic system. The estimated daily releases from Bluewater Lake were used to calibrate both the steady-state and transient periods (see “MODFLOW-NWT, Calibration Approach and Datasets” section). Monthly precipitation rates on the surface of Bluewater Lake were estimated using the gridded, hydrometeorological datasets described in the “Climate Distribution Parameterization” section by computing the average monthly mean precipitation for the 13 grid cells (plate 1) used to simulate the lake. A steady-state precipitation rate was estimated as the average of the monthly precipitation rates for the transient simulation period. Monthly evaporation rates from the surface of Bluewater Lake were estimated using the “Monthly Average Pan Evaporation” data at the Laguna station (Western Regional Climate Center, 2021). A steady-state evaporation rate was estimated as the average of the monthly evaporation rates for the transient simulation period.

According to annual groundwater monitoring reports for Homestake mill (Hydro-Engineering, LLC, 1999, 2000, 2001, 2002, 2003; Homestake Mining Company of California and Hydro-Engineering, LLC, 2004, 2009, and 2013), groundwater was extracted from wells located upgradient of the uranium mill to prevent groundwater flow into the mill site where contaminated groundwater was being remediated and discharged to a nearby surface-water channel from approximately 1993 through 2011. Discharge of extracted groundwater to the surface-water channel adjacent to Homestake mill was simulated in MODFLOW-NWT using the SFR2 package. Monthly volumes of groundwater extraction from the wells located upgradient of Homestake mill were estimated from the annual groundwater monitoring reports (Hydro-Engineering, LLC, 1999 and 2001, fig. 2.1–2; Homestake Mining Company of California and Hydro-Engineering, LLC, 2014 and 2016) and converted to a daily volumetric rate by dividing by the number of days per month. The estimated daily volumetric rates of groundwater extraction were assumed to be equal to the daily surface-water discharge and were specified as an inflow to the stream network in the first reach of the segment adjacent to Homestake mill (ISEG variable equal to 113) (plate 1).

According to Risser (1982), the community of Grants waste-water treatment plant began discharging treated effluent to the Rio San Jose around 1957. Around 1978 a flow of about 1 to 2 ft³/s of treated effluent became perennial from the waste-water treatment plant downstream to the area around Horace Springs (fig. 1; Risser, 1982). The community of Grants stopped regular discharge from the waste-water treatment plant to the Rio San Jose in 1992 (Wolf, 2016). Treated effluent is now applied to ponds on the community of Grants golf course (Roca Honda Resources, LLC, 2013). Discharge of the treated effluent was simulated in MODFLOW-NWT using the SFR2 package. From 1957 through 1977, discharge of treated effluent to the Rio San Jose was specified at a volumetric rate of 0.5 ft³/s (half of the minimum of the range of flow reported by Risser [1982] beginning in 1978). From 1978 through 1992, discharge of treated effluent to the Rio San Jose was specified at a volumetric rate of 1.5 ft³/s (average of the minimum and maximum of the range of flow estimated by Risser [1982]). Discharge of treated effluent to the Rio San Jose was specified as an inflow to the stream network in the first reach of the segment downstream from the community of Grants (ISEG variable equal to 179) (plate 1).

Spring Delineation and Spring Discharge Routing

Spring locations were identified at 168 locations within the active MODFLOW-NWT area from the USGS NWIS (USGS, 2019b). Discharge from these springs was simulated using the SFR2 package. Thirty springs were in stream reaches and were simulated as the groundwater discharge to the existing reaches. Thus, discharge from these 30 springs could directly contribute to streamflow in the stream network. Discharge from Horace Springs was simulated at two existing reaches on the Rio San Jose (ISEG variable equal to 207 and IREACH variable equal to 3 and 4). The remaining 138 springs were added to the SFR2 package as a segment with one reach, with the length of the stream channel within each reach set equal to the length of the side of a model grid cell (500 m) and the elevation of the top of the springbed (STRTOP variable) within each reach was set as the adjusted HRU land surface elevations described in the “Precipitation-Runoff Modeling System” section. Four grid cells contained two springs each and the springs in each grid cell were simulated with the same segment and reach in the SFR2 package. Except for Ojo del Gallo, discharge from these 138 springs was assumed to be entirely removed from the hydrologic system by not specifying a downstream segment (OUTSEG variable equal to 0). During the nongrowing season (October through March), discharge from Ojo del Gallo was routed to an existing segment (ISEG variable equal to 193) in the SFR2 package that is a tributary to the Rio San Jose to simulate that the spring was allowed to flow along its natural course (discharging to a marsh and then to a channel that flowed east to the Rio San Jose) (Gordon, 1961; Baldwin and Anderholm, 1992). During the growing season (April through September) and in the steady-state simulation period, discharge from Ojo del Gallo was removed from the model to simulate the historical use of the springflow for irrigation in the region of the community of San Rafael (Gordon, 1961; Baldwin and Anderholm, 1992). Deep percolation of Ojo del Gallo springflow that was applied to the land surface for irrigation in the region of the community of San Rafael was simulated using injection wells in the Well package (Harbaugh and others, 2000). The amount of deep percolation of Ojo del Gallo springflow was based on data from a preliminary historical surface-water and groundwater use compilation for the Bluewater Declared Underground Water Basin (data are available on request from the NMOSE) (see “MODFLOW-NWT, Water Use Estimation and Distribution” section).

Model Layering and Parameterization from Geologic Framework Model

A three-dimensional geologic framework model of the basin was constructed and included 18 geologic units, the locations of faults that truncate and offset the geologic units, and the locations of volcanic vents and igneous dikes that penetrate all units (Sweetkind and others, 2020;Sweetkind and Galanter, in press) (plate 1). For the purposes of this study, the 18 geologic units from the geologic framework model were grouped into 9 hydrogeologic units (numbered 1 to 9 from youngest to oldest) on the basis of similar age, lithology, and reported hydraulic properties (hydraulic conductivity and storage) (table 4). MODFLOW-NWT was discretized vertically into eight layers (numbered 1 to 8 from highest to lowest), with the top and bottom elevations of each layer assigned using the top and bottom elevations of the 9 hydrogeologic units defined in this study.

Model Layering

The top elevation of layer 1 corresponded to the land surface elevation developed for PRMS (see “Precipitation-Runoff Modeling System” section), except for the 13 grid cells used to simulate Bluewater Lake where the top elevation of layer 1 was set equal to the maximum lake stage of 7,406 ft above NAVD 88 for “excellent boating” reported by the New Mexico Energy, Minerals, and Natural Resources Department (2022). The bottom elevation of layer 1 corresponded to the bottom elevation of the uppermost hydrogeologic unit at each model row and column location, except for the 13 grid cells used to simulate Bluewater Lake or where the Precambrian hydrogeologic unit was the uppermost hydrogeologic unit. For the 13 grid cells used to simulate Bluewater Lake, the bottom elevation of layer 1 corresponded to the elevation of the dry lakebed (7,358 ft above NAVD 88), which was estimated from lake surface areas and stages for “excellent boating” and “fair boating” reported by the New Mexico Energy, Minerals, and Natural Resources Department (2022). Where the Precambrian hydrogeologic unit was the uppermost hydrogeologic unit, the bottom elevation of layer 1 was defined as 150 ft below the top elevation of layer 1. This approach for defining the bottom elevation of layer 1 created a continuous layer 1 across the entire active MODFLOW-NWT area (fig. 1) and allowed the Precambrian hydrogeologic unit to be incorporated in MODFLOW-NWT.

The top and bottom elevations of MODFLOW-NWT layers 2 through 8 were assigned on the basis of the presence or absence of successively older hydrogeologic units. The Precambrian hydrogeologic unit was only used to define the top and bottom elevations of MODFLOW-NWT layers where it was the uppermost hydrogeologic unit, in which case layers 1 and 2 were each defined as 150 ft thick. The 150 ft thickness of MODFLOW-NWT layers associated with the Precambrian hydrogeologic unit were assigned to minimize potential numerical instabilities caused by thin layers and to produce a total Precambrian hydrogeologic unit thickness of 300 ft, which corresponds to the maximum thickness of weathered Precambrian rocks in the study area that transmit water (Baldwin and Rankin, 1995). The vertical discretization of MODFLOW-NWT is demonstrated in two cross sections: A–A′ (fig. 4) oriented north to south and B–B′ (fig. 5) oriented west to east (locations of cross-section lines are shown on plate 1). Cross-section A–A′ extends from the northwestern portion of the study area in the Chaco slope subdivision of the San Juan structural basin, crosses the Continental Divide, passes through Bluewater Lake, the Zuni uplift, and El Malpais National Monument, and ends in the southwestern portion of the study area. Cross-section B–B′ extends from the Continental Divide in the western portion of the study area, passes through the Zuni uplift, Ojo del Gallo, the Rio San Jose, the southern flank of Mount Taylor and associated tributaries to the Rio San Jose, the Rio Puerco fault zone, and ends in the western portion of the Middle Rio Grande Basin. The cross sections illustrate the large range of MODFLOW-NWT layer elevations, the offsets of layers associated with faulting, and the discontinuous nature of hydrogeologic units throughout the study area.

The number of active layers at a model row and column location in MODFLOW-NWT was equal to the number of hydrogeologic units 1 to 8 (excluding the Precambrian hydrogeologic unit) that had a nonzero thickness in the geologic framework model (figs. 4 and 5). An arbitrary thickness value of 5 m was assigned where a geologic framework model unit was absent in the subsurface because of an erosional unconformity and between two units that were present (Sweetkind and others, 2020). In MODFLOW-NWT, the thickness of these absent hydrogeologic units was set as an active layer and had hydraulic properties assigned with the adjacent hydrogeologic unit (below) that was present. Because of the discontinuous nature of hydrogeologic units throughout the study area (table 4) (Sweetkind and others, 2020), MODFLOW-NWT layers represent multiple hydrogeologic units.

Model Parameterization

The locations of the faults, volcanic vents, and igneous dikes were incorporated in MODFLOW-NWT as areas of potential resistance to horizontal flow using the Horizontal Flow Barrier package (Hsieh and Freckleton, 1993) (plate 1). Eight calibration parameters were defined for the hydraulic characteristic within the Horizontal Flow Barrier package, one each for volcanic vents and igneous dikes and six for faults. The six fault parameters were created by grouping the faults on the basis of age (offsets pre-Quaternary units or Quaternary units) and fault strike (northwest, north, and northeast), with the groupings designed to combine faults with similar age and tectonic setting on the assumption that these similarities might result in similar fault behavior in restricting horizontal flow.

Multiple hydrogeologic unit zones were defined for 5 of the 9 hydrogeologic units by geologic framework model unit and previous studies that delineated distinct zones with similar hydraulic properties (horizontal hydraulic conductivity, vertical anisotropy, specific storage, and specific yield) within units (table 4 and figs. 6–10). The remaining 4 hydrogeologic units (Mancos Shale, Dakota Formation, Chinle Formation, Precambrian) were each assigned a single hydrogeologic unit zone (homogeneous hydraulic properties within the entire extent of the hydrogeologic unit) largely because of a lack of data from previous studies to inform spatial variability of hydraulic properties.

The Cenozoic and Lower Permian-Pennsylvanian hydrogeologic units were subdivided into 4 and 3 hydrogeologic unit zones, respectively (table 4, figs. 6 and 10), on the basis of the geologic framework model unit that had the greatest thickness in each hydrogeologic unit. The Mesaverde Group and Morrison-Entrada hydrogeologic units were subdivided into 6 and 7 hydrogeologic unit zones, respectively (figs. 7 and 8), on the basis of the geologic framework model unit that had the greatest thickness in each hydrogeologic unit and reported transmissivity zones (Stone, 1981; Frenzel and Lyford, 1982). Two transmissivity zones (50–100 and 100–350 feet squared per day [ft²/d]) for the Gallup Sandstone in the San Juan structural basin that were developed by Stone (1981) were used to subdivide the Mesaverde Group hydrogeologic unit. Five transmissivity zones (less than 50, 50–150, 150–250, 250–400, and 400–500 ft²/d) for the Morrison Formation in the San Juan structural basin that were developed by Frenzel and Lyford (1982) were used to subdivide the Morrison-Entrada hydrogeologic unit.

The SAGA hydrogeologic unit was subdivided into five hydrogeologic unit zones (fig. 9) on the basis of previous studies by Kelley and Wood (1946), Baars (1962), and Baldwin and Anderholm (1992). Hydrogeologic unit zone number 3 represents the typical San Andres Limestone, with the lower half to two-thirds usually a thin- to medium-bedded alternation of dolomites, shales, siltstones, and sandstones, overlain by a middle white sandstone and an overlying typically massive limestone. The San Andres Limestone in hydrogeologic unit zone number 3 was never exposed to subaerial erosion, so there are few solution-enhanced openings or caverns. Hydrogeologic unit zone number 3 is relatively distant from the faulting in the area of the community of Grants, so there is relatively limited secondary fracturing. Baldwin and Anderholm (1992) estimated the transmissivity of hydrogeologic unit zone number 3 to be 800 ft²/d on the basis of aquifer-test results from a single well. Hydrogeologic unit zone number 1 is like the typical San Andres Limestone in hydrogeologic unit zone number 3, except the San Andres Limestone thins such that more of the SAGA hydrogeologic unit hydraulic properties represent the Glorieta Sandstone. Baldwin and Anderholm (1992) reported transmissivity values ranging from 30 to 280 ft²/d and estimated a representative transmissivity of 140 ft²/d in hydrogeologic unit zone number 1. In hydrogeologic unit zone number 2, cavernous zones, faults, and fractures are common, and the large range in reported transmissivity (100–450,000 ft²/d) reflects well completions in various configurations of the faulted, fractured, and cavernous carbonate terrain. Baldwin and Anderholm (1992) estimated an average transmissivity of 50,000 ft²/d for hydrogeologic unit zone number 2. The San Andres Limestone changes abruptly to interbedded dark red shale and sandstone toward the north and northeast in hydrogeologic unit zone number 4. Although thin carbonates are usually present in hydrogeologic unit zone number 4, they are a minor feature. Baldwin and Anderholm (1992) estimated a representative transmissivity of 800 ft²/d in hydrogeologic unit zone number 4, but no wells were analyzed. The San Andres Limestone changes abruptly to evaporites toward the southeast in hydrogeologic unit zone number 5, thickening markedly into an evaporite basin (Baars, 1962). Kelley and Wood (1946) described a thick series of evaporites in the Sierra Lucero uplift of hydrogeologic unit zone number 5 that may be equivalent to the lower dolomites and clastics in the Zuni uplift. On the basis of lithologic and hydrologic data from seven wells, the transmissivity for hydrogeologic unit zone number 5 was estimated by Baldwin and Anderholm (1992) to range from 10 to 150 ft²/d.

Water Use Estimation and Distribution

Water use in the study area has varied both by category of use and quantity throughout the simulation periods of the RSJIHM. Although available data are sparse for some areas and years, effort was taken to accurately represent the water use changes that took place within the active MODFLOW-NWT area. Water use, including pumping, deep percolation, injection of water and seepage to groundwater from tailing piles related to uranium mill operations, and seepage to groundwater from water supply treatment plants, was simulated in MODFLOW-NWT using the Well package. Pumping was simulated using extraction wells, and deep percolation, injection, and seepage were simulated using injection wells (fig. 11). In addition, diversions of surface water from Bluewater Creek below Bluewater Dam for irrigation within the Bluewater-Toltec Irrigation District were simulated in MODFLOW-NWT as diversions from the SFR2 package.

Calibration Approach and Datasets

Calibration of MODFLOW-NWT was performed using BeoPEST (version 13.6) (Schreüder, 2021). BeoPEST was run using the USGS Yeti supercomputer (U.S. Geological Survey Advanced Research Computing, 2021). The parameters adjusted during calibration were

Horizontal hydraulic conductivity, vertical anisotropy, specific yield, and specific storage parameters were adjusted for groups of grid cells corresponding to the hydrogeologic unit zones in each model layer. General-head boundary hydraulic conductance was adjusted for 7 parameter groups of grid cells, consisting of 5 parameter groups that corresponded to the five regions of groundwater flow into the Middle Rio Grande Basin used in McAda and Barroll (2002) and 1 parameter group each for groundwater flow across the northern and western boundaries of MODFLOW-NWT (plate 1). Eight parameters were defined for the hydraulic characteristic within the Horizontal Flow Barrier package (see “Model Layering and Parameterization from Geologic Framework Model” section). Parameters within the SFR2 package were defined separately for all stream network grid cells within each of the 15 subbasins, the Rio Puerco upstream from its confluence with the Arroyo Chico, the Rio Puerco downstream from its confluence with the Arroyo Chico, Horace Springs, and all other springs. A single lakebed leakance parameter was defined for the 13 grid cells representing Bluewater Lake.

Adjustments to parameter values during calibration were based on the ability of model simulated values to reproduce measured, estimated, or reported values for several observation groups: (1) streamflow, (2) head, (3) springflow at Ojo del Gallo, (4) springflow at Horace Springs, (5) surface-water releases from Bluewater Lake, and (6) surface-water diversions for irrigation within the Bluewater-Toltec Irrigation District. Observation weights within PEST were applied to give relatively equal importance to each observation group. Within observation groups, values were given the same PEST observation weights, except for head values, which were given variable weights based on the assumed quality of the data. Higher weights were assigned to head data assumed to be of higher quality, such as data from NWIS that have documented accuracies; lower weights were assigned to head data assumed to be of lower quality, such as data from NMOSE POD records that do not have documented accuracies and may not represent static heads. The ability of MODFLOW-NWT to reproduce measured, estimated, or reported values is discussed in the “Model Performance” section.

Hydraulic Properties of Model Layers

The horizontal hydraulic conductivity, vertical anisotropy, specific yield, and specific storage values by hydrogeologic unit and hydrogeologic unit zone that resulted from the calibration of the transient RSJIHM are presented in table 7. During calibration, the hydraulic conductivity and storage property values were allowed to vary within the range of reported hydraulic property values obtained from previous studies (table 7). The calibrated horizontal hydraulic conductivity ranged from 0.0065 to 80 ft/d. The lowest horizontal hydraulic conductivity value was associated with the Cenozoic hydrogeologic unit zone number 3 where Miocene to Pliocene basaltic to andesitic lava flows with local dacites near Mount Taylor of the Neogene volcanics geologic framework model unit are predominant. The highest horizontal hydraulic conductivity value was associated with SAGA hydrogeologic unit zone number 2 where cavernous zones, faults, and fractures are common. The calibrated vertical anisotropy (ratio of horizontal hydraulic conductivity divided by vertical hydraulic conductivity) ranged from 8.0 to 61,000. The lowest vertical anisotropy value was associated with Cenozoic hydrogeologic unit zone numbers 2 and 4 where the basaltic to andesitic lava flows and alluvial deposits of the Quaternary basalts geologic framework model unit and the Oligocene to early Pleistocene sedimentary and volcanic partly consolidated conglomerate sandstone, shale, volcanics of the Santa Fe Group basin fill geologic framework model unit are predominant, respectively. The highest vertical anisotropy value was associated with the Mesaverde Group hydrogeologic unit zone number 3 where shale, siltstone, and coal of the Crevasse Canyon Formation geologic framework model unit are predominant. Specific yield was calibrated to range from 0.020 to 0.29, and specific storage was calibrated to range from 1.0×10^–08 to 4.0×10^–05 per foot.

General-Head Boundary Hydraulic Conductance

The General-Head Boundary package hydraulic conductance values that resulted from the calibration of the transient RSJIHM are presented in table 8. During calibration, the hydraulic conductance values were allowed to vary within a range of 10 orders of magnitude because of a lack of data to constrain the values. The calibrated hydraulic conductance ranged from 0.000036 ft²/d along the western boundary to 68 ft²/d along the Rio Salado region of the Middle Rio Grande Basin boundary used in McAda and Barroll (2002) (plate 1).

Horizontal Flow Barrier Hydraulic Characteristic

The Horizontal Flow Barrier package hydraulic characteristic values that resulted from the calibration of the transient RSJIHM are presented in table 9. During calibration, the hydraulic characteristic values were allowed to vary within a range of 10 orders of magnitude because of a lack of data to constrain the values. The calibrated hydraulic characteristic ranged from 0.00000000010 per day for faults that offset pre-Quaternary units and have a general northwest strike to 0.070 per day for volcanic vents (plate 1). The faults with the highest calibrated hydraulic characteristic were those that offset pre-Quaternary units and have a general northeast strike and extend from El Malpais National Monument to the north and northeast of the community of Grants. The faults with the lowest calibrated hydraulic characteristic tended to be those that are along and near the boundaries of the active MODFLOW-NWT area.

Streambed and Springbed Hydraulic Conductivity

The streambed and springbed hydraulic conductivity values that resulted from the calibration of the transient RSJIHM are presented in table 10. During calibration, the streambed and springbed hydraulic conductivity values were allowed to vary within the range of reported horizontal hydraulic conductivity values obtained from previous studies (table 7). The calibrated streambed hydraulic conductivity ranged from 0.000010 ft/d in subbasin 11 to 77 ft/d in subbasin 2. The calibrated springbed hydraulic conductivity ranged from 10 ft/d at Horace Spring to 27 ft/d at all other springs.

Streambed and Springbed Elevation

The elevation of the top of the streambed and springbed as derived from the GSFLOW-Arcpy toolkit (set as 1 ft below the adjusted HRU land surface elevations described in the “Precipitation-Runoff Modeling System” section) was calibrated by optimizing reductions to the elevation (table 11). During calibration, the elevation reductions were allowed to vary from 1 ft to the reduction that would cause the bed of at least one reach in a parameter group to fall below the uppermost model layer. The calibrated elevation reductions ranged from 1.0 ft in subbasin 12 to 82 ft in subbasin 1. The highest calibrated elevation reductions tended to be for reaches in higher elevation subbasins where greater slopes result in a larger range of elevations within the 500- by 500-m grid cells (USGS, 2017). Therefore, the horizontal discretization of the RSJIHM was too coarse in these areas to accurately represent the streambed elevation without adjustment.

Lakebed Leakance

The calibrated lakebed leakance within the Lake package representing Bluewater Lake was 0.000000024 day⁻¹. During calibration, the lakebed leakance value was allowed to vary within a range of 10 orders of magnitude because of a lack of data to constrain the value.

Mean Monthly

The Nash-Sutcliffe model efficiency coefficient between measured and simulated mean monthly streamflow ranged from −1.9 (“not satisfactory”) at the outlet of subbasin 5 to 0.99 (“very good”) at the outlet of subbasin 3 and was 0.53 (“satisfactory”) for all streamgages (table 15). The percent bias between measured and simulated mean monthly streamflow ranged from −100 at the outlets of subbasins 5 and 6 to 51 at the outlet of subbasin 2 (both “not satisfactory”) and was 2.8 (“very good”) for all streamgages. The R² between measured and simulated mean monthly streamflow ranged from 0.0025 (“not satisfactory”) at the outlet of subbasin 7 to 1.0 (“very good”) at the outlet of subbasin 3 and was 0.53 (“not satisfactory”) for all streamgages. The R² could not be calculated at the outlets of subbasins 5 and 6 because simulated mean monthly streamflow was zero. Overall, these quantitative measures indicate spatial variability in the performance of the RSJIHM at simulating measured mean monthly streamflow throughout the active MODFLOW-NWT area.

Hydrographs of measured and simulated mean monthly streamflow at the 12 streamgages are shown in fig. 13. In the Zuni uplift, the RSJIHM was unable to simulate the pronounced peak in measured mean monthly streamflow during March and April at the outlet of subbasin 1 but did a much better job of simulating the peak in streamflow at the outlet of subbasin 2. The RSJIHM closely simulated the measured mean monthly streamflow at the outlet of subbasin 3, which is about 2 mi downstream from Bluewater Dam and thus provides an indication of the performance of the RSJIHM at simulating surface-water releases from Bluewater Lake. At the outlet of subbasin 4, which is about 8 mi downstream from Bluewater Dam, the RSJIHM simulated the overall patterns of streamflow corresponding to surface-water releases from Bluewater Lake from May through September, but oversimulated streamflow from December to March. The RSJIHM was unable to simulate the relatively low (less than 2.5 ft³/s) measured mean monthly streamflow at the outlets of subbasins 5 and 6, which drain the Chaco slope subdivision of the San Juan structural basin (fig. 2).

Measured mean monthly streamflow at the outlet of subbasin 7, which drains an area extending from the Continental Divide to near the community of Grants (fig. 2), shows a bimodal pattern with pronounced peaks during April and May and July through September (fig. 13). The RSJIHM simulated a bimodal pattern of mean monthly streamflow at the outlet of subbasin 7 that differed from the measured streamflow in the timing and magnitude of the peaks. The RSJIHM simulated a relatively steady mean monthly streamflow at the outlet of subbasin 8, which drains a relatively small area near the community of Grants, that differed from the pronounced measured peak during July through September. Measured mean monthly streamflow at the outlet of subbasin 9 shows a bimodal pattern with pronounced peaks during April and May and July through September. The RSJIHM simulated a bimodal pattern of mean monthly streamflow at the outlet of subbasin 9 that differed from the measured streamflow in the timing and magnitude of the peaks. The RSJIHM simulated the overall patterns of measured mean monthly streamflow at the outlet of subbasin 10, which drains the eastern side of Mount Taylor, with a pronounced peak during July through September, but simulated about half of the maximum measured streamflow peak during August. The RSJIHM was generally able to simulate the bimodal patterns of measured mean monthly streamflow at the outlets of subbasins 11 and 12 near the terminus of the basin, with lower streamflow during May and June and a pronounced peak during July through October.

Monthly Mean

The Nash-Sutcliffe model efficiency coefficient between measured and simulated monthly mean streamflow ranged from −0.37 (“not satisfactory”) at the outlet of subbasin 7 to 0.98 (“very good”) at the outlet of subbasin 3 and was 0.43 (“not satisfactory”) for all streamgages (table 15). The percent bias between measured and simulated monthly mean streamflow ranged from −100 at the outlets of subbasins 5 and 6 to 54 at the outlet of subbasin 2 (both “not satisfactory”) and was 8.7 (“good”) for all streamgages. The R² between measured and simulated monthly mean streamflow ranged from 0.00010 (“not satisfactory”) at the outlet of subbasin 9 to 0.98 (“very good”) at the outlet of subbasin 3 and was 0.42 (“not satisfactory”) for all streamgages. The R² could not be calculated at the outlets of subbasins 5 and 6 because simulated monthly mean streamflow was zero. Overall, these quantitative measures indicate spatial variability in the performance of the RSJIHM at simulating measured monthly mean streamflow throughout the active MODFLOW-NWT area.

Annual Mean

The Nash-Sutcliffe model efficiency coefficient between measured and simulated annual mean streamflow ranged from −3.6 (“not satisfactory”) at the outlet of subbasin 6 to 0.99 (“very good”) at the outlet of subbasin 3 and was 0.55 (“satisfactory”) for all streamgages (table 15). The percent bias between measured and simulated annual mean streamflow ranged from −100 at the outlets of subbasins 5 and 6 to 52 at the outlet of subbasin 2 (both “not satisfactory”) and was 7.9 (“good”) for all streamgages. The R² between measured and simulated annual mean streamflow ranged from 0.00010 (“not satisfactory”) at the outlet of subbasin 9 to 1.0 (“very good”) at the outlet of subbasins 3 and 11 and was 0.58 (“not satisfactory”) for all streamgages. The R² could not be calculated at the outlets of subbasins 5 and 6 because simulated annual mean streamflow was zero. Overall, these quantitative measures indicate spatial variability in the performance of the RSJIHM at simulating measured annual mean streamflow throughout the active MODFLOW-NWT area.

Hydrographs of measured and simulated annual mean streamflow at the 12 streamgages are shown in fig. 14. In the Zuni uplift, the RSJIHM simulated the overall patterns of measured annual mean streamflow peaks and reductions, with undersimulation of the magnitudes of the peaks at the outlet of subbasin 1 and a mixture of undersimulation and oversimulation of the magnitudes of the peaks at the outlet of subbasin 2. The RSJIHM closely simulated the measured annual mean streamflow at the outlet of subbasin 3. The RSJIHM generally simulated the overall patterns and magnitudes of measured annual mean streamflow at the outlet of subbasin 4, but oversimulated streamflow by greater than 3 ft³/s during 1968, 1971 and 1972. The RSJIHM was unable to simulate the relatively low (less than 2 ft³/s) measured annual mean streamflow at the outlets of subbasins 5 and 6.

The RSJIHM was generally unable to simulate the patterns and magnitudes of measured annual mean streamflow at the outlet of subbasin 7 (fig. 14). The RSJIHM generally simulated the trends of measured annual mean streamflow at the outlet of subbasin 8 but was unable to simulate the peakiness of measured streamflow. The RSJIHM generally simulated the trends of measured annual mean streamflow at the outlet of subbasin 9 but was unable to simulate the magnitudes of the streamflow peaks. At the outlet of subbasin 10, the RSJIHM simulated the overall patterns of measured annual mean streamflow, but only simulated about a quarter of the measured streamflow peak during 1988. The RSJIHM was generally able to simulate the patterns and magnitudes of measured annual mean streamflow at the outlets of subbasins 11 and 12.

During the high precipitation verification year, 1983, the RSJIHM undersimulated measured annual mean streamflow by greater than 6 ft³/s at the outlets of subbasins 7 and 9 (fig. 14). The RSJIHM performed better at simulating measured annual mean streamflow during 1983 at the outlets of subbasins 8, 10, and 12. During the low precipitation verification year, 2018, the RSJIHM oversimulated measured monthly mean streamflow at the outlet of subbasin 9, with the average of the monthly mean streamflow during January through September simulated at 4.1 ft³/s and measured at 3.1 ft³/s. Measured streamflow data were not available at any other streamgages during these verification periods.