Potentiometric-Surface Construction (Kriging)

The potentiometric surfaces of the Gallup Sandstone were developed using a kriging spatial interpolation algorithm (Ritchie and Pepin, 2020). This methodology uses a dataset that includes geographic coordinates and groundwater elevations from wells located within the study area, identified as the surface and subsurface extent of the Gallup Sandstone (fig. 1; Kernodle, 1996). These data are used to develop a semivariogram model of the correlation versus distance relation, which in turn is used by the kriging algorithm to estimate groundwater elevations throughout the study area. The kriging workflow was developed in the R programming language (R Core Team, 2019). This workflow used the “ObsNetwork” package, a spatial optimization algorithm (Fisher, 2013). ObsNetwork uses a geostatistical technique known as universal kriging to interpolate the groundwater elevation at unmeasured sites within the study area. The general implementation of this workflow involves (1) spatial trend modeling, (2) semivariogram development, and (3) kriged surface development and assessment.

Spatial trend modeling addresses the kriging assumption that the mean value of the data being estimated does not change when shifted in space. Groundwater elevations typically have spatial trends that violate this assumption. To correct for this, a trend model is subtracted from the groundwater elevations, and the semivariogram and kriged surface are developed for the residuals. The trend is added back to the kriged surface to acquire the final kriged groundwater-elevation estimates.

For the semivariogram development, empirical and theoretical semivariograms were developed for the new surface. Semivariance for all pairs of monitoring sites are calculated as one-half the variance of the difference between the residual groundwater elevations at the sites. The empirical variogram is computed by grouping the semivariance values into bins based on the distance (lag) between monitoring site pairs and averaging the semivariance values within each bin. Kriging relies on a continuous theoretical semivariogram that is fit to the empirical variogram. This fitting process requires the consideration of semivariogram shape, which is strongly influenced by three components called the nugget, sill, and range. The nugget is the semivariance value at a lag of zero, which typically is related to noise in the data and is user specified by visual inspection. The sill is the semivariance value at which the mean semivariance stops changing with increasing distance. The range is the distance at which the sill is reached and indicates the distance over which data are correlated. Circular, exponential, gaussian, linear, and spherical models were all considered in fitting the theoretical semivariogram.

For the kriged surface development and assessment, universal kriging was used to interpolate groundwater elevation and quantify uncertainty associated with the interpolated values. Interpolation was performed at a uniform grid size measuring 500 by 500 ft. The grid resolution was chosen to be smaller than the minimum distance between sites and the groundwater model grid size.

Leave-one-out cross validation was used to test the performance of the theoretical semivariogram model (Fisher, 2013). This method removes each measurement point from the dataset one by one and estimates the groundwater elevation at the removed site by using the theoretical semivariogram model developed for the entire dataset to krige with the remaining data. The estimation error from this method is calculated as the difference between the observed and estimated groundwater elevation at the omitted site.

Kriging is advantageous because it produces measurements of uncertainty associated with the estimates of water-level elevation. Uncertainty is presented as standard error, or the square root of the estimated variance (Fisher, 2013). Standard error is essentially a scaled version of the distance to the nearest measurement point, meaning standard error is zero at the water level observation points, given a nugget of zero, and increases with distance.

Potentiometric surfaces were kriged for both modern data (data range between 2016 and 2020, table 1) and for historical data (data range between 1954 and 1989, table 2). These new potentiometric surfaces were used to derive drawdown maps and calibrate the groundwater model. Drawdown was calculated as the difference between the two kriged potentiometric surfaces.

Description of Numerical Groundwater Flow Model

The numerical groundwater flow model (“model”) (Shephard and others, 2023) was constructed using MODFLOW-2005 version 1.12.0 (Harbaugh, 2005), a well-documented and commonly used three-dimensional finite-difference groundwater flow model developed by the USGS that simulates steady and nonsteady flow in irregularly shaped flow systems. The MODFLOW code was executed in part in Groundwater Vistas (GWV) version 6 (Rumbaugh and Rumbaugh, 2011), a graphical user interface that can run source code for several different groundwater models and particle tracking simulators. GWV also has built-in visualization functions as well as calibration and sensitivity analysis functions. The model was also run directly from the MODFLOW-2005 program executable for various runs and simulations. MODPATH version 7.2.001 (Pollock, 2016), a USGS particle tracking model for MODFLOW also executed in GWV and directly from the program executable, was used to assess relative groundwater flow paths and traveltimes within the constructed MODFLOW model.

Spatial and Temporal Discretization

The active model area is a rectangle measuring approximately 25 x 31 miles (mi), approximately centered around CCNHP (fig. 1). These model dimensions were chosen to ensure that pumping operations inside or outside of the proposed model domain produce minimal changes in the potentiometric surface at the model boundaries. The grid cells are 1,500 by 1,500 ft, and the model grid consists of 119 rows and 104 columns. Cell size was determined on the basis of layer thicknesses (described in the “Model Layer Description” methods section). Given the relatively small thicknesses of the Gallup Sandstone in the model domain, the Mancos Shale model layers were split into multiple layers of smaller cells to avoid the creation of relatively thin cells in the Gallup Sandstone in MODFLOW. The grid is rotated approximately 50.3 degrees counterclockwise from grid north to align with the general direction of groundwater flow and the geographic extent of the Gallup Sandstone (Kernodle, 1996).

The model simulates steady-state conditions relative to the modern potentiometric surface in the Gallup Sandstone (see the “Water-Level Measurements” methods section). Under steady-state conditions, the potentiometric surface does not change with time. The results presented in this study comparing the historical and modern potentiometric surfaces in the Gallup Sandstone indicate that the hydrologic system is not in steady state. However, construction of a model with time-varying inflows and outflows was beyond the scope of this study. A model of this type could change in response to well injection and extraction of water resources over time. Incorporation of time-varying flows into further versions of the model may be useful to understand the temporal dynamics of influences on the Chaco well.

Model Layer Descriptions

The model domain features one Gallup Sandstone layer and two main Mancos Shale units. The model specifically only simulates the deep groundwater system of these units below CCNHP and does not simulate flow in shallower aquifers or groundwater flow of the entire San Juan structural Basin or greater region. The upper and lower Mancos Shale units within the model domain each range from 600 to 1,300 ft in thickness, whereas the Gallup Sandstone ranges from 100 to 200 ft in thickness. Large differences in thickness between adjacent layers could cause potential issues when running MODFLOW. To account for this, the Mancos Shale units are split into multiple layers so that no layer is drastically thicker than the Gallup Sandstone layer. The upper and lower Mancos Shale model units consist of 5 layers each, which together with the Gallup Sandstone layer, compose a total of 11 model layers. The upper Mancos Shale bottom geometry was assumed to be the same as the Gallup Sandstone top geometry. The top geometry for the upper Mancos Shale unit was derived from available structural data (Kernodle and others, 1989) on the top elevation of the Point Lookout Sandstone, the stratigraphic unit that overlies the upper Mancos Shale in the model domain. Using well log data, a relative thickness of the Point Lookout Sandstone in the model domain was estimated. This thickness was subtracted from structural top data of the Point Lookout Sandstone to construct the elevation and geometry of the top of the Upper Mancos Shale. The Gallup Sandstone layer was constructed using available data on the elevation of the structural top of the unit from Kernodle and others (1989). These data were digitized using a geographic information system (GIS) and then brought into the model to produce the top geometry.

The thickness of the Gallup Sandstone was based on available well log data. This thickness was subtracted from the structural top of the Gallup Sandstone to produce the bottom geometry for the layer. The lower Mancos Shale unit top geometry was assumed to be the same as the bottom of the Gallup Sandstone. The bottom geometry for the lower Mancos Shale unit was estimated from the available structural data from Kernodle and others (1989) describing the top of the Dakota Sandstone, the stratigraphic unit located below the lower Mancos Shale.

Hydrologic Boundaries

The northeast boundary was simulated as a no-flow boundary because updated potentiometric surfaces of the Gallup Sandstone and additional references to previously published surfaces (Frenzel and Lyford, 1982; Kernodle, 1996) show that this boundary is roughly parallel to groundwater flow. Additionally, this boundary is the northeasternmost extent of the Gallup Sandstone where it pinches out, also justifying it as a no-flow boundary. The northwest boundary was simulated as a constant head boundary in the Gallup Sandstone layer, set to the head of the modern potentiometric surface. The southeast and southwest boundaries were simulated as flux-in boundaries in the Gallup Sandstone layer. Fluxes into the model were initially estimated from the hydraulic gradient, cell thickness, and initial hydraulic conductivity of the layer and were later calibrated (described in the “Calibration and Sensitivity Analysis” section). The model top and bottom were simulated as no flow boundaries because, as noted in previous work, flow into and out of the top and bottom of the Mancos Shale was assumed to be insignificant. There is generally no recharge to the Gallup Sandstone in this region from precipitation or surface water, indicated by the old age of the water (Linhoff and others, 2023) and the isolating thickness of the Mancos Shale (Stone and others, 1983).

Calibration and Sensitivity Analysis

Groundwater model calibration was conducted by matching simulated and observed heads. Calibration targets were created by extracting values from a GIS raster of the modern potentiometric surface derived from data collected for this study. These values were extracted at evenly spaced intervals throughout the model. A spacing of 20,000 ft was used for extraction points, which produced 69 evenly spaced points in the domain at a relatively high density for the size of the model area. For the calibration process, the hydraulic heads extracted at these points were used as the observed values, which were compared to the model simulated values at these points in the Gallup Sandstone (model layer 6).

Calibration and uncertainty analysis were carried out manually in GMV and automatically using the model-independent Parameter ESTimation (PEST) and uncertainty analysis program (Doherty, 2018a, b). PEST is a commonly used program for automatic environmental model calibration and uncertainty analysis. The program implements traditional parameter estimation approaches on the basis of a few parameters by running a model as many times as needed while adjusting parameters to reach a best fit between the simulated and observed data by minimizing weighted least squares errors. Parameter ranges for PEST input were set on the basis of the same literature values used in the manual calibration process (table 3). Adjustments to parameter values during calibration were based on the ability of model simulated heads, which were extracted from the model using the Head-Observation package (Hill and others, 2000), to reproduce the heads extracted from the modern potentiometric surface. Head observation weights within PEST were applied to give equal importance to each head observation. Model files, including those used for PEST calibration, can be found in the model archive associated with this report (Shephard and others, 2023).

Parameter inputs for the model and parameter calibration ranges were determined through values reported in previous works. A summary of available data is provided in tables 3 and 4. Primarily, the parameters altered during the calibration process were the vertical and horizontal hydraulic conductivity and the flux rates into the model, simulated as wells in MODFLOW. Specific yield and storativity were not calibrated because the model was run as steady state. Porosity was estimated from literature values for use in MODPATH, but it was not calibrated. However, this report does provide groundwater traveltime analysis for a range of porosity values. Parameters in the 10 Mancos Shale layers and single Gallup Sandstone layer were calibrated independently. The upper and lower Mancos Shale were treated as separate homogeneous and anisotropic zones, although they were calibrated at the same time with the same parameters. The flux-in boundary along the southwest and southeast parts of the model was divided evenly into three sections along the southwest and southeast boundaries, for a total of six reaches that simulated the boundary condition of groundwater flux into the model.

Sensitivity Analysis

A model sensitivity analysis of various model input parameters was conducted to identify the parameters that have the greatest influence on model performance. The sensitivity analysis presented in this report was conducted following model calibration; however, an initial sensitivity analysis was also conducted (not reported) to determine which model parameters should be calibrated. The sensitivity analysis was conducted using the auto-sensitivity function in GWV. Parameter values were varied by ±50 percent of the calibrated model parameter in 10 increments, resulting in 11 model runs that vary a single parameter 10 percent between each run, including the initial calibrated parameters. This process is repeated for each parameter of interest. For each of these model runs, the absolute residual mean (RM) between the observed and simulated data is calculated and plotted against the incremental parameter values. These RM plots are used to assess the sensitivity of each model parameter from a magnitude and rate of change perspective. The model parameters used in the sensitivity analysis were the horizontal and vertical conductivities of the Mancos Shale and Gallup Sandstone layers, and the flux-in boundary for the Gallup Sandstone model layer. For vertical and horizontal hydraulic conductivities, the parameters for all 10 Mancos Shale model layers were changed at the same time to assess Mancos Shale parameter sensitivity, and the Gallup Sandstone hydraulic conductivity was changed independently of the Mancos Shale, during different sensitivity runs. The sensitivity of the flux-in boundaries of each reach of the Gallup Sandstone was tested separately, resulting in six separate sensitivity analyses for the flux-in boundary. Additionally, the constant head boundary sensitivity was tested.

Model Application and Simulations

The calibrated model was used to simulate the groundwater flow paths of water that reach the Chaco well, and how proximity of potential HF operations to the Chaco well capture zone may influence the Chaco well. Modeling simulations primarily made use of forward and reverse particle tracking using MODPATH. The model and the scenarios only consider advective groundwater transport effects of potential conductivity changes caused by HF. Four simulations carried out following model calibration are described here.

Simulation 4. Potential Influence of Hydrologic Fracturing-Zones Forward Tracking, Two Pumping Rates, Two Hydraulic Conductivities

In simulation 3, hydraulic conductivity values were changed for entire layers representing Mancos Shale. In simulation 4, more localized changes to the Mancos Shale hydraulic conductivity were simulated. Potential future placement of HF wells is unknown, but this scenario provides a representative hypothetical example. Hydraulic conductivities in parameter zones within and adjacent to the capture zones determined in simulation 1 were increased to simulate hypothetical HF wells. These circular hypothetical HF zones were assigned a 10,000-ft radius to represent possible horizontally drilled distances and the zone of influence of an individual fracturing well (Meng, 2015). The zones were chosen specifically to be either within the two previously defined capture zones (Chaco well pumping rates of 1,200 and 10,000 ft³/d) or outside of the two capture zones. The model was run with Chaco well pumping rates of 1,200 and 10,000 ft³/d. For either scenario, the horizontal and vertical hydraulic conductivities were altered in the zones in layer 4 to 1 foot per day (ft/d) and 100 ft/d, respectively, resulting in a total of four separate model runs. These two values were used, given the uncertainty of the magnitude of hydraulic conductivity change caused by HF. Forward tracking was implemented by releasing particles in each of the cells occupying the zones and tracking their flow paths. A representative schematic of the 11 model layers and altered hydraulic conductivity values is shown in figure 3.

Model Calibration

Horizontal and vertical hydraulic conductivity values in the Mancos Shale and Gallup Sandstone and groundwater flux into the model were calibrated using modern head observations (table 5). The calibrated parameters are determined as being best fit by PEST, assuming the constant head boundary is equal to the head of the modern potentiometric surface in those model cells. Calibration ranges used as input for PEST were determined from well logs and literature sources (tables 2–4). Calibrated PEST hydraulic conductivity values in the Mancos Shale (horizontal hydraulic conductivity [HK] approximately [~] 0.0001 ft/d, vertical hydraulic conductivity [VK] ~0.0001 ft/d) were almost identical to parameters reported in Kernodle (1996) (HK = 0.0001 ft/d, VK = 0.0001 ft/d). HK was slightly lower than reported by Frenzel and Lyford (1982) (HK = 0.000864–0.00864 ft/d), but vertical hydraulic conductivity in this study was higher (Frenzel and Lyford [1982] reported VK = 0.0000000864 ft/d). Horizontal hydraulic conductivity in the Gallup Sandstone in this study (HK ~0.41 ft/d) was in-between the values of Kernodle (1996) (HK = 1 ft/d) and Frenzel and Lyford (1982) (HK = 0.1–0.14 ft/d), and vertical hydraulic conductivity was identical to that of Kernodle (1996) (VK = 0.002 ft/d).

The model was calibrated to the updated modern surface hydraulic heads (described in the “Calibration and Sensitivity Analysis” section), where calibration targets were extracted from an overlain grid at 20,000 ft intervals from the kriged modern potentiometric surface. Total model mass balance errors were very small (0.01 percent), and the model successfully converged with head change and residual convergence criteria of 0.001 ft and 1 ft³/d, respectively, using the Preconditioned Conjugate-Gradient solver. The mean residual (simulated minus observed) of the model was 2 ft, and the absolute residual mean was 33 ft, which is less than 5 percent of the total head drop across the active model area (~730 ft). The standard deviation of residuals was 42 ft. The minimum residual was −82 ft, and the maximum residual was 106 ft. Previously published models in the vicinity of this study area reported residuals of as much as 150 ft between model-derived potentiometric surface and observed head as “good fit,” given the total range of heads in the study area (Frenzel and Lyford, 1982), although some residuals were not reported (Kernodle, 1996). The R² between the simulated and observed head values was 0.93, indicating good agreement between simulated and observed values (fig. 9). However, comparison of the linear fit of the simulated against observed values, as well as the plot of the residuals versus observed values, indicate a slight tendency to overpredict at lower observed head values and a slight tendency to underpredict at higher observed head values (kriged head versus residuals R² = 0.14, slope = −0.10) (figs. 9 and 10).

Though RMs are generally low, there is some spatial heterogeneity among residual values (fig. 11). Residuals are generally most positive (>50 ft) along the northeast boundary of the active model area, near the extent of the Gallup Sandstone, and most negative just south of CCNHP (less than −50 ft). Large residuals could be due to kriging interpolation error. Residuals generally are closest to zero (−10 to 10 ft) along the northwest boundary, which is a constant head boundary, and scattered elsewhere between areas of higher residuals. Numerical groundwater models are simplified representations of natural systems, and these spatial biases are likely an artifact of the actual spatial heterogeneity of the vertical and hydraulic conductivity. The model treats the Mancos Shale and Gallup Sandstone layers as spatially homogeneous for the purpose of simplicity and because of the lack of data or other evidence to populate the model heterogeneously.

Several additional factors contribute to model uncertainty, including model error associated with the groundwater-level measurements and the standard error of the modern kriged surface (fig. 6). Parameter ranges for calibration were based largely on previously reported values (table 3), which were determined either physically or through parameters previously calibrated for earlier models, which also have uncertainty and various limitations. MODPATH, which was used to generate simulated groundwater flow paths, potentially yields simplified flow paths or uncertain traveltimes because it considers only advective groundwater flow and ignores preferential flow features and dispersive flow features of water-bearing formations (Baca, 1999). Additionally, subsurface features such as faults and anthropogenic preferential pathways such as abandoned well casings were not represented. Pumping wells other than the Chaco well were not included because of a lack of data but were assumed individually to have negligible influence because of their smaller pumping volumes. As noted earlier, there may be up to 10,000 acre-ft of groundwater used in McKinley County per year. The model simulations in this report are simplistic representations of how HF operation may influence subsurface processes and advective groundwater flow paths and thus may neglect other complex HF influences, as described in the “Introduction” section (Birdsell and others, 2015).

Sensitivity Analysis

Sensitivity analysis results are plotted as a percentage of parameter change from the calibrated parameter value versus the RM between simulated and observed heads (fig. 12A–C). The parameters included in the sensitivity analysis were horizontal and vertical hydraulic conductivities of the Mancos Shale and Gallup Sandstone model layers and the fluxes into the Gallup Sandstone. The sensitivity of particle traveltimes to porosity is discussed in the “Traveltimes” section of this report.

For hydraulic conductivity, only the HK RM in the Gallup Sandstone had high sensitivity (fig. 12B). The Gallup HK RM ranges from approximately 33 to 243 ft. The VK in the Gallup Sandstone as well as the HK and VK of the Mancos Shale are generally insensitive. Flux-in for reaches 3 and 4 were the most sensitive of the six flux-in parameters (fig. 12C). The reach 3 RM ranges from approximately 33 to 44 ft, and reach 4 ranges from 33 to 65 ft. The constant head boundary was found to be the most sensitive parameter in the entire sensitivity analysis, with the RM ranging from 33 to 3,033 ft (fig. 12A).

The Mancos Shale VK and HK have lower sensitivity, and as a result, altering the hydraulic conductivity through HF may have a lesser impact on model performance than initially anticipated. Reaches 3 and 4 are located along the steepest head gradient of the model, so the greater rate of change in RM for these reaches makes sense, because they contribute greater amounts of flux into the model than other reaches. The constant head boundary is the most sensitive parameter but also has the most quantitatively uncertain results through kriged standard error.

General Flow Directions and Gradients

Flow directions in the model are largely a function of head gradient derived from the modern kriged potentiometric surface, which may be subject to interpolation error. The flow directions and behavior in the model are also in part a function of the specified model boundaries. Water generally enters the model through the flux-in boundary along the southwest and southeast boundaries, flows parallel along the northeast boundary, and exits the model along the northwest boundary. Previously published potentiometric head contours (Frenzel and Lyford, 1982; Kernodle, 1996) are in relative agreement with this study in terms of gradient direction but are not digitized and are difficult to quantitatively compare with new results (Stewart, 2018). However, the hydraulic gradients in this study in the vicinity of the CCNHP (minimum and maximum values range ~730 ft) tend to be steeper (greater range) compared to previous contour maps of head (Kernodle, 1996), potentially from increased or continued groundwater extraction from the Gallup Sandstone aquifer, coupled with limited aquifer recharge.

For most of the active model area, flow vectors in cross section view indicate that groundwater tends to flow horizontally between model cells, with limited mixing between the Mancos Shale and Gallup Sandstone layers. However, model results indicate slight vertical mixing of water between the Mancos Shale and the Gallup Sandstone, primarily near the Chaco well where pumping is occurring (fig. 13A, B). Hydraulic head is lowered substantially in the model cells near the Chaco well, especially those cells extending vertically to the top and bottom of the active model area. This change in hydraulic gradient causes groundwater to flow from the Mancos Shale layers to the Gallup Sandstone layer, and the gradient increases when the pumping rate is increased from the reported 1,200 ft³/d to the theoretical 10,000 ft³/d. The model framework did not account for all pumping wells in the area and does not have the available data to account for potential abandoned or compromised wells that could act as permeable pathways between the Mancos Shale and Gallup Sandstone or other surrounding units (Birdsell and others, 2015). These gradients and pathways presumably exist elsewhere in the model domain, which could promote additional vertical mixing between the Gallup Sandstone and Mancos Shale layers.

Capture Zones

Capture zones were assessed by forward tracking particles placed along the flux-in boundary of the model domain until their capture by the Chaco well. Capture zones for the hypothetical Chaco well pumping rate of 10,000 ft³/d are shown in figure 14A, B in two-dimensional space and are depicted in three-dimensional space by the color of the pathways. These are representative figures of particles released in model layers 5 and 7 (Mancos Shale), representing average depths approximately 125 ft above and 175 ft below the Gallup Sandstone, respectively. The capture zone for either pumping rate tends to originate at the southeast flux-in boundary and take an arched northwest-trending path through the model until the particles are captured. Increasing the pumping rate substantially increases the capture zone width. The average width of the capture zone is approximately 1 mi for the 1,200 ft³/d pumping rate (not shown) and approximately 12 mi for the 10,000 ft³/d pumping rate.

A total of 2,442 particles were released in all 11 model layers in each cell along the southwest and southeast flux-in boundaries. Of these particles, 41 and 387 were captured by the well at pumping rates of 1,200 and 10,000 ft³/d, respectively. The remaining particles did not pass through or terminate in the cell containing the Chaco well. Particle traveltimes ranged between approximately 28,500 and 34.2 million years for captured particles at a pumping rate of 1,200 ft³/d and between approximately 24,300 and 76.4 million years at a pumping rate of 10,000 ft³/d. The longer traveltimes for the larger pumping rate are due to the resulting larger capture zone that is created. Carbon isotope dating of water from three Gallup Sandstone wells within the modeled domain indicated an age range from ~43,000 to >52,000 years before present, the latter being greater than the method can resolve (Linhoff and others, 2023). Although these sampled waters travel from an unknown recharge area, and the active model area only occupies a portion of the basin, the order of magnitude of ages and traveltimes suggest a level of consistency and accuracy between both methods. These traveltimes are theoretical and dependent on running the model in steady-state simulation until all particles reach their endpoint at the calibrated parameter set and initially determined assumed porosity values (table 6). Traveltimes are calculated for a range of porosity values suggested in the literature (examined in the “Traveltimes” section).

Particles released in layer 5 primarily remained in the upper Mancos Shale layers (layers 2–4) stratigraphically above the Gallup Sandstone and did not reach the Gallup Sandstone until they were within 1,500 ft of the Chaco well. Similarly, the particles released in layer 7 remained in the lower Mancos Shale layers (layers 7–10) stratigraphically below the Gallup Sandstone and did not reach the Gallup Sandstone until they were within 1,500 ft of the Chaco well. These particle tracking pathways indicate that mixing between the Gallup Sandstone and Mancos Shale likely occurs primarily near pumping wells or near very permeable pathways such as compromised well casings or fractured rock.

Traveltimes

Two additional model runs were conducted that used the calibrated model parameters and Chaco well pumping rate of 1,200 ft³/d with the minimum and maximum porosity values for the Mancos Shale and Gallup Sandstone from the literature (table 4; Martin, 2005) to assess potential minimum and maximum groundwater traveltimes. Similar to the process for capture zones, a total of 2,442 particles were released in all 11 model layers in each cell along the southwest and southeast flux-in boundaries and forward tracked in steady-state simulation. Traveltimes are generally slower with higher porosity, given that flow paths are less tortuous, and fluid velocity through a porous media is mathematically expressed as Darcy flux divided by the porosity. For minimum porosity values, traveltimes ranged from approximately 5,700 years to 30 million years. For maximum porosity values, traveltimes ranged between approximately 52,000 years and 36.3 million years. The minimum traveltimes reflect particles that originated in and primarily flowed through the Gallup Sandstone, whereas maximum values reflect particles that originated in and primarily flowed through the Mancos Shale.

These simulated values represent traveltimes from the boundary of the active model area, not total traveltimes from the time water enters the subsurface. Water typically enters the Gallup Sandstone from surface outcrops far to the south and southeast of the Chaco well and the active model area (Stone and others, 1983), meaning traveltimes to the Chaco well are longer than those determined by the model. However, computed model traveltimes for minimum porosity values and the assumed additional travel from outside of the active model area are consistent with computed carbon isotope ages of samples collected in the Gallup Sandstone within the modeled domain (43,000 to >52,000 years before present; Linhoff and others, 2023).

Hydraulic Fracturing Scenarios

Several scenarios representing different HF operations were investigated (table 6). Hydraulic conductivities of each Mancos Shale layer were altered individually by several orders of magnitude to simulate the maximum distance and depth that potential HF-induced changes in hydraulic conductivity could potentially influence the Chaco well. In table 6, each row represents a model layer (Mancos Shale = layers 1–5, 7–11; Gallup Sandstone = layer 6), and each column represents the hydraulic conductivity assigned to simulate potential HF operation. For each scenario model run, Mancos Shale layer hydraulic conductivity (horizontal and vertical) was altered separately by layer (one layer per model run), the model was executed, and then particles were tracked in reverse to record their path in space and time from model cells surrounding the Chaco well. If a particle passed through the altered Mancos Shale layer, this was an indication that HF operation in that particular layer could potentially influence the Chaco well. This simulation was conducted for Chaco well pumping rates of 1,200 and 10,000 ft³/d.

For the 1,200 ft³/d pumping rate, increased HK and VK in the Mancos Shale began to influence the Chaco well when hydraulic conductivity reached 1 ft/d, which is five orders of magnitude above the hydraulic conductivity determined through model calibration (HK and VK ~0.0001 ft/d). Hydraulic conductivity has been determined to change by three to five orders of magnitude from HF in laboratory settings (Gehne and Benson, 2019). Model results indicate that if HF alters hydraulic conductivity by less than five orders of magnitude, there will be no influence on the Chaco well. For a hydraulic conductivity of 1 ft/d, only the first Mancos Shale layers above and below the Gallup Sandstone (layers 5 and 7) influenced flow to the Chaco well, and for hydraulic conductivities above 10 ft/d, three layers above and below the Gallup Sandstone influenced flow to the Chaco well (layers 3–5 and 7–9). For hydraulic conductivities of 100 ft/d or more, four layers above and below the Gallup Sandstone influenced flow to the Chaco well (layers 2–5 and 7–10). However, the top and bottom layers of the model (layers 1 and 11) did not pass any particles that were captured by the Chaco well for any hydraulic conductivity values (as much as 1,000 ft/d) at the 1,200-ft³/d pumping rate.

For the 10,000-ft³/d pumping rate, flow to the Chaco well began to be influenced by high HK and VK in the Mancos Shale layers at a hydraulic conductivity of 0.1 ft/d, and more layers above and below the Gallup Sandstone began to influence flow to the Chaco well. At hydraulic conductivities of 1 ft/d or more, all of the layers above and below the Gallup Sandstone influenced water flowing to the Chaco well.

Changing the horizontal and vertical hydraulic conductivities of the entire layer may have a different effect on the behavior of subsurface flow than changing them in smaller, isolated areas. Changing the hydraulic conductivity of the Mancos Shale over an entire model layer to values above the Gallup Sandstone hydraulic conductivity may promote flow through the Mancos Shale layers, rather than promoting the movement of water from the Mancos Shale to the Gallup Sandstone. This behavior may not be indicative to how HF may influence flow pathways and traveltimes. Altering the hydraulic conductivity in small areas in the Mancos Shale layers, similar to how an HF well may alter hydraulic conductivity, effectively creates preferential pathways for groundwater flow between the Mancos Shale and the Gallup Sandstone.

For example, figure 15A, B shows the pathlines and starting locations in layer 4 for the particles captured by the Chaco well that originated from smaller areas (having a radius of about 10,000 ft) of altered hydraulic conductivity to emulate HF fracture zones. These particles passed vertically through the Mancos Shale, reached the Gallup Sandstone, and continued to pass through the Gallup Sandstone before reaching the Chaco well. Altering the hydraulic conductivity allowed groundwater to reach the Gallup Sandstone earlier than it otherwise would have and ultimately decreased the traveltimes to the Chaco well. This behavior occurred at modeled hydraulic conductivities of 1 and 100 ft/d, which suggests that a substantial change in hydraulic conductivity in a small area of the Mancos Shale may provide preferential pathways for groundwater travel to the Gallup Sandstone. These scenarios were run for both Chaco well pumping scenarios (1,200 and 10,000 ft³/d). Particles placed in smaller altered zones within the original capture zone (with original calibrated model parameters based on known, present-day, or observed conditions) in horizontal space (fig. 15A, B) were also captured by the Chaco well. Particles placed outside of the original capture zone (in horizontal space) were still not captured by the well, indicating that changes in hydraulic conductivity have limited influence on advective flow pathways outside of the HF zone.

These emulated HF zones decreased groundwater traveltimes from the Mancos Shale to the Gallup Sandstone by providing increased hydraulic conductivity and preferential pathways between the layers. Mean traveltimes for particles captured by the Chaco well that originated in altered HF zones were approximately 14,300 years for Chaco well pumping rates of 10,000 ft³/d compared to 18.5 million years for unaltered zones. The decrease in traveltime was primarily due to the HF zones changing the travel paths of the particles directly into the Gallup Sandstone.