Model Extent and Discretization

The MODFLOW model domain consists of a uniform grid of 680 rows and 619 columns of 50-m-square cells. The model has four layers that represent the glacially deposited sediments and bedrock, with layers 1, 2, and 3 representing surficial glacially-deposited sediments—including outwash, alluvium, kame, lacustrine sediments, and till—and layer 4 representing the underlying bedrock. Layers 1–3 varied in thickness across the model domain and were spatially discontinuous. The bedrock unit, represented by layer 4, was set to a uniform thickness of 30 m. Model top elevations were derived from lidar data (U.S. Geological Survey, 2015) by sampling the mean elevation within each cell and mapping that value to the model grid. Model layer elevations were developed using a variety of information, including well construction reports, boring logs, geologic cross-section data, Soil Survey Geographic Database (SSURGO) data, and lidar data. The methods used to develop model layer elevations and thicknesses are described in detail by Finkelstein and others (2022), and the data sets are made available in the USGS data release by Woda and others (2022).

The lithology of the valley-fill aquifer systems (represented in model layers 1–3) is complex. The individual units that comprise these layers have irregular geometries that interfinger and pinch out (fig. 2, for example). The lateral continuity of the units also is variable—valley-fill facies present in one valley reach may be different than those present in another. These complex aquifer geometries were simulated using the MODFLOW 6 idomain array to specify active (idomain=1), inactive (idomain=0), and vertical pass-through (idomain=−1) cells. Vertical pass-through cells were used to route water in the model between the uppermost active and the next highest active layer; for example, if the geologic units represented by model layers 2 and 3 pinched out along the edge of a valley, vertical passthrough cells applied to layers 2 and 3 allow water to be exchanged between layer 1 and layer 4 directly, without the need for the intermediate layers to be included. This improves the computational efficiency of the model and the accuracy of the model’s physical representation of the flow system.

The geographic extent of the model was established using the topographical divides defined by 5 adjacent 12-digit hydrologic unit code (HUC) drainage basins. These drainage basins include the Upper Neversink River (020401040305), Gumaer Brook-Basher Kill (020401040201), Basher Kill (020401040203), Beer Kill (020200070503), and Sandburg Creek (020200070504). The combined extent of these drainage basins defined the active model domain with two modifications: (1) the model boundary downgradient from the Neversink Reservoir was adjusted about 0.25 mile southeast—from the top of the Neversink dam structure to its base—to improve model stability; and (2) the model extent in Basher Kill River drainage basin was truncated south of Wurtsboro, N.Y., downgradient from any priority wells to improve the computational efficiency of the groundwater-flow model. The spatial extent of the model domain is shown in figure 1.

Boundary Conditions

Hydrologic boundary conditions represent flows of water into and out of the model domain. Boundary conditions in this model include areal recharge from net infiltration percolating below the root zone, groundwater flow across the model perimeter in the valley-fill areas, stream-aquifer interactions, and groundwater extraction at wells.

Groundwater Flow Across Model Perimeter

Groundwater flow across the model perimeter in valleys were simulated using the MODFLOW 6 specified head boundary package (CHD). CHD boundaries for this model were developed from two previously published regional groundwater-flow models (Zell and Sanford, 2020a; Zell and Sanford, 2020b). Zell and Sanford (2020a) produced 75 steady-state MODFLOW 6 groundwater-flow models that span the contiguous United States and simulate the long-term mean surficial groundwater system. The NS–RO model domain intersects parts of two Zell and Sanford (2020a) regional models. The discretization packages (DIS) and binary head output files from the two intersecting Zell and Sanford (2020a) models were combined to create a single synthetic “parent” model. NS–RO CHD values were extracted from the “parent” model head output using the modflow-setup Python package. The locations of CHD cells in the groundwater-flow model and their corresponding values are shown on figure 4.

Hydraulic Properties

Hydraulic properties in the model were defined by zones based on lithology data and surficial geology mapping. The initial hydraulic conductivity values assigned to these zones (and initial upper and lower calibration bounds) were estimated using published regional aquifer test data (Como, 2013; Smith and Williams, 2020). The lithologic zones are shown, by layer, in panel “B” of figures 5–12. Initial (prior) hydraulic conductivities, aggregated by geologic zones, are summarized in table 3. Information about the history matching process and the development of posterior hydraulic properties are presented in the “Parameter Estimation by Ensemble History Matching” section.

Observation Data

An observation dataset was developed to compare model outputs with field measurements of the system to facilitate ensemble-based history matching. Historical water-level measurements were obtained from the National Water Information System (NWIS) and NYSDEC water well completion reports (New York State Department of Environmental Conservation, 2020). Only a single water-level measurement was available in the NWIS database (at site 414525074360601; U.S. Geological Survey, 2020a). Where present, multiple hydraulic head measurements were averaged across time to derive a single, steady-state value for a given location. In instances where water-level measurements were reported as depth to water below land surface, water-level elevations were computed by subtracting depth to water values from the model top elevation. This subtraction is accompanied by substantial uncertainty given the steep terrain and the resolution of the elevation data informing the model’s top. In total, 448 hydraulic head targets were used for history matching. The single NWIS measurement includes more accurate elevation information than the NYSDEC measurements; however, the need to average over time causes enough uncertainty that the entire group of hydraulic head measurements was assumed to have the same uncertainty. This assumed hydraulic head measurement uncertainty was used to assign weights to hydraulic head observation for the history matching. During history matching, hydraulic head observations were compared against transmissivity-weighted hydraulic heads simulated by the model at the location of the observation well.

Steady-state base flows were computed from continuous records at two USGS streamgages 01436500 (Neversink River at Woodbourne, N.Y.), and 01366650 (Sandburg Creek at Ellenville, N.Y.). Long-term mean base-flow index (BFI; the average proportion of streamflow that is base flow) was computed using the USGS StreamStats web application (U.S. Geological Survey, 2020b). Mean annual flow rates were multiplied by the BFI values to compute an estimate of steady-state base flow (table 2). At streamgage 01366650, mean annual flow rates were computed using the entire period of record (October 1, 1957–September 29, 1977). At streamgage 01436500, only data collected after the construction of the Neversink Reservoir were used to compute mean annual flow rate (October 1, 1954–September 30, 1993). It was assumed that base flow computed over these periods is representative of long-term steady state conditions.

Land-surface elevations were also used as inequality targets on a grid every 50 model cells in each direction. We initially met a challenge common to this setting of hydraulic head solutions above the land-surface elevation (in other words, flooding) in early versions of the model, particularly before history matching. The locations of calibration targets are shown on figures 13–15.

We grouped the observations in three groups: flux, hydraulic head, and land surface. We assigned initial weights based on a general assumption of the observation quality (table 4), which corresponded to a coefficient of variation of 10 percent for flux observations, and standard deviation of 5 and 10 m for the hydraulic head and land-surface observations, respectively. These assumptions are qualitative and result in an imbalanced objective function (the squared sum of residuals between model outputs and observations, each multiplied by their weight, and referred to as the ‘PHI’ value; fig. 16A). Subjectively, we judged that a reasonable balance of the objective function is 20 percent for hydraulic heads, 30 percent for flux, and 50 percent for land surface. This results in the balanced objective function in figure 16B with adjusted weights shown in table 4. The adjusted weights imply a reversal of standard deviation between hydraulic head and land-surface observations (11.2 and 2.4 m, respectively). Furthermore, the implied standard deviation for the two flux observations decreased dramatically because of this rebalancing; however, without adjusting the weights, the history matching algorithm would largely ignore the flux observations, allowing great variability and not constraining the overall water balance of the model. While many hydraulic head observations were available, they are individually subject to substantial uncertainty as discussed previously.

Stepwise Estimation of Model Parameters

Model parameters that are estimated in the history matching process are those for which we assume are informed by the observations. In choosing which parameters to allow to be adjustable through history matching, we acknowledge that any parameters not allowed to be estimated through history matching are assumed perfectly known. In this section, we outline the strategy for selecting the parameters to adjust and estimate in the history matching process. We then describe the approach for evaluating the uncertainty of parameters and observations to set initial parameter values and then the algorithm for history matching.

Parameterizing the Model

We followed the technique proposed by White and others (2020) to use multiple scales of variability in the parameterization. We typically used multipliers to allow adjustment at multiple scales—for example, single multipliers applied at the model layer scale combined with zone-wise multipliers applied to local areas (zones) within the model, in the case of parameters representing hydrogeologic properties (for example, hydraulic conductivity). Within zones, multipliers on pilot points (Doherty, 2003) were used to allow for finer-scale variability in hydraulic conductivity and recharge. Table 5 shows the categories of parameters and the number of adjustable parameters for each type. For horizontal and vertical hydraulic conductivity, we assigned zone-based homogeneous multipliers based on the mapped lithologic zones. Additionally, pilot points were placed every 20 cells in layers 1 and 4 and every 5 cells in layers 2 and 3. The higher frequency in layers 2 and 3 was to accommodate the irregular geometry of the geological zonation in those layers. The hydraulic conductivity zones and pilot point networks are displayed in panel “C” of figures 5–12. In a similar multiscale approach, we implemented a single recharge multiplier for the entire recharge array to allow for relatively large-scale adjustment, augmented by pilot points to allow for finer scale adjustment. The recharge pilot points were also placed every 20 model cells. The recharge pilot points are collocated with the layer 4 hydraulic conductivity pilot points as displayed in figure 17. The final parameters we implemented are the pumping rate in each well (as a multiplier) and a multiplier for the SFR streambed conductance by reach.

Table 5.    

Summaries of initial parameter values and bounds during the prior Monte Carlo, expanded Monte Carlo, and iterative ensemble smoother history matching, aggregated by parameter type.

[iES, iterative ensemble smoother; SFR; MODFLOW 6 Streamflow Routing package]

Table 5.    Summaries of initial parameter values and bounds during the prior Monte Carlo, expanded Monte Carlo, and iterative ensemble smoother history matching, aggregated by parameter type.

Parameter		Prior Monte Carlo			Expanded Monte Carlo			iES start
Type	Count	Initial value	Lower bound	Upper bound	Initial value	Lower bound	Upper bound	Initial value	Lower bound	Upper bound
Pilot point multipliers
Horizontal hydraulic conductivity: Layer 1	202	1	0.442759	1.55724	1	0.442759	1.55724	0.942248 to 1.05538	0.75596 to 0.866352	1.13672 to 1.31561
Horizontal hydraulic conductivity: Layer 2	560	1	0.442759	1.55724	1	0.442759	1.55724	0.934627 to 1.06573	0.741177 to 0.874849	1.1411 to 1.33064
Horizontal hydraulic conductivity: Layer 3	473	1	0.442759	1.55724	1	0.442759	1.55724	0.935518 to 1.07049	0.744879 to 0.870171	1.13199 to 1.33375
Horizontal hydraulic conductivity: Layer 4	470	1	0.442759	1.55724	1	0.442759	1.55724	0.928528 to 1.0651	0.754228 to 0.865706	1.15635 to 1.33452
Vertical hydraulic conductivity: Layer 1	202	1	0.442759	1.55724	1	0.442759	1.55724	0.935312 to 1.06418	0.759673 to 0.879614	1.15972 to 1.3271
Vertical hydraulic conductivity: Layer 2	560	1	0.442759	1.55724	1	0.442759	1.55724	0.933356 to 1.06665	0.746704 to 0.880208	1.13367 to 1.32458
Vertical hydraulic conductivity: Layer 3	473	1	0.442759	1.55724	1	0.442759	1.55724	0.899499 to 1.06398	0.750198 to 0.876427	1.13878 to 1.32987
Vertical hydraulic conductivity: Layer 4	470	1	0.442759	1.55724	1	0.442759	1.55724	0.934159 to 1.08109	0.749258 to 0.884273	1.14991 to 1.36895
Recharge: Highest active layer	470	1	0.8	1.2	1	0.8	1.2	0.978882 to 1.01934	0.912194 to 0.953566	1.04398 to 1.09877
Multipliers by other attributes
Recharge for entire array	1	1	0.9	1.1	1	0.9	1.1	0.992711	0.953912	1.02459
SFR streambed conductance by reach	915	1	0.1	10	1	0.1	10	0.799918 to 1.29901	0.349099 to 0.637529	1.6714 to 2.77499
Pumping rate by well	31	1	0.8	1.2	1	0.8	1.2	0.987671 to 1.01061	0.91623 to 0.949433	1.03776 to 1.09147
Horizontal hydraulic conductivity by zone	42	1	0.181818 to 0.75	1.25 to 1.81818	0.266667 to 50	0.00266667 to 0.5	26.6667 to 5000	0.195096 to 65.7042	0.0380862 to 12.6817	1.03622 to 314.01
Vertical hydraulic conductivity by zone	42	1	0.181818 to 0.75	1.25 to 1.81818	0.266667 to 50	0.00266667 to 0.5	26.6667 to 5000	0.252877 to 67.3208	0.0611225 to 14.1789	1.48613 to 309.894

Using Uncertainty Analysis to Set Initial Values

An initial prior Monte Carlo analysis was completed on the parameterized model. A parameter ensemble was generated assuming a normal distribution of values and that the parameter bounds represent a spread of 4 standard deviations (in other words, 95 percent confidence) about the initial parameter values, as shown in table 5 in the section “Prior Monte Carlo.” The bounds were intentionally limited to keep parameters within initial ranges based on literature values. This was particularly problematic for the bedrock hydraulic conductivity values that assume deep, generally unweathered, and sparsely fractured bedrock when, in reality, the bedrock units include a shallow zone of greater weathering and fracturing that must be represented by the same unit.

The initial parameter bounds were evaluated by comparing model outputs from the ensemble of parameter samples in the prior Monte Carlo against observation data. We used an ensemble approach to sample the “noise” of the observations such that each realization uses a set of observation values sampled from a normal distribution with the mean being the reported observation value and the standard deviation being a user-defined variance representative of the uncertainty of the observations (table 4).

The left column in figure 18 shows histograms for several representative observations. The blue histograms represent model outputs from the ensemble of parameter samples in prior Monte Carlo. The orange histograms represent the ensemble of observation values, with the vertical solid line representing the reported observation value, and the vertical dashed lines representing ±1 standard deviation of the observation ensemble. While observations A and B show that flux is overestimated in one location and underestimated in another, observations C and D show that with the initial parameter assumptions, hydraulic heads are often too high; in the case of observation D, hydraulic head is 25 to 75 m higher than land surface. As a result, we expanded the bounds for hydraulic conductivity parameters by zone and reran the Monte Carlo analysis resulting in generally lower heads (the middle panels of fig. 18) although with relatively minor changes to the flux residuals. Based on this improvement, we used the parameter ensemble from the expanded Monte Carlo as both starting parameter values and upper and lower bounds for the history matching. The 25th, 50th, and 75th quantiles from the expanded Monte Carlo parameter ensemble were assigned as the lower bound, initial parameter value, and upper bound for history matching, respectively. These values are summarized in table 5 in the section “iES start.”

History Matching with the Iterative Ensemble Smoother (iES)

We selected the iES (White, 2018) as implemented in PEST++ (White and others, 2021) for history matching. This technique builds on the Monte Carlo approach by evaluating an ensemble of parameter values, yielding a residual vector and calculating a weighted objective function for each member of the parameter ensemble. From this information, the algorithm calculates an updated ensemble of parameters to evaluate, and iterations are done until a posterior conditioned ensemble is obtained. This results in an ensemble of objective function values from iteration to iteration (fig. 19). In figure 19, the red curve represents a base ensemble member that is taken as representative when a single realization is desired. The iES approach does not depend on local linearization of the inverse problem so inequality observations are possible. In this case, we defined water levels at the land surface on a regular network as inequality observations for the land-surface observation group. If the inequality is not invoked (in other words, if the model output hydraulic head elevation is lower than the collocated observation of land-surface elevation) then the weight assigned to this observation is 0.0 and thus the residual does not play a role in the objective function calculation. However, if the simulated hydraulic head exceeded these elevations, a quadratic penalty in the objective function was incurred at their assigned weights. Furthermore, the observation noise distributions allow flexibility in the algorithm to consider the uncertainty of the observations which, as discussed previously, is considered significant in this case where we have many observation values but none of them may be perfectly representative (for example, uncertainty in their location, timing versus the steady state model, and other aspects).

We selected the results of the second iES iteration as optimal based on a subjective judgement that the objective function had lowered to an acceptable level of fit while still retaining variability in the ensemble thus still capturing uncertainty in the system. As an additional step, we applied a rejection sampling process in which objective function values that are unacceptably high are trimmed from the ensemble. As shown in figure 20, a PHI cutoff was selected that removed anomalously high PHI results and results in an approximately normal distribution of PHI results. Figure 21 shows the correspondence between model outputs and observations, broken down by group. These results are shown to have little bias and low variance relative to the uncertainty of the observations. In figure 21, the blue symbols show fit for the representative base model while the bars indicate the 95-percent credible interval as calculated from the trimmed ensemble. Figure 18 shows the posterior ensemble of model outputs compared with the ensemble of observation values, the distribution of outputs and observations are consistent with their uncertainty. Generally, for the examples shown, the distributions of model outputs and observation uncertainty overlap with similar mean values.

Figures 13–15 show the spatial patterns of the residuals for the hydraulic head, flux, and land-surface observation groups, respectively. In each figure (13-15), the color of the symbol indicates the residual value (defined as observed value minus simulated value) for the base ensemble member, and the size of the symbol indicates the standard deviation of the residuals for all the ensemble members. In figure 14, one residual value is positive and one is negative. In figure 13, the residuals are positive and negative as well, highlighting that while there are some clusters of simulated values exceeding observed and vice versa, there is not a consistent spatial pattern that might indicate a structural issue with the model or the estimated parameters (for example, if all valleys were simulated as flooded). This mixed pattern of residuals is the desired outcome from history matching and indicates that the results are generally unbiased. The pattern is different for the land-surface residuals in figure 15. In this case, there are very few negative residuals and nearly all are biased such that the simulated results are lower value than the observed. Because these are inequality constraints and were intended to prevent parameters from driving the solution to simulate hydraulic heads higher than land surface, this bias is the desired outcome.

Base Realization

The final values of the “base” realization represent the central tendency of the posterior parameter distribution (White and others, 2021). This realization differs from the other realizations that are obtained by making a stochastic sample using the bounds of the initial parameter values. Instead, the base realization is initially made up of the specific starting values provided to PEST++. As a result, this realization starts at values that are considered “most likely” by the users. This realization is then subjected to adjustment during each iteration of the iES algorithm, but the central tendency typically remains, and this single realization typically is less variable than the other realizations. In this section, we present the representative base realization model properties, including hydraulic conductivities and boundary conditions that are representative of the optimal iES ensemble. We also discuss model results that were simulated using base realization parameters, including the simulated water balance, water table, and streamflow. For the analysis of areas supplying recharge, however, the entire ensemble is used to be sure that rather than focusing on central tendency, we are exploring the range of outcomes consistent with prior knowledge of the parameters and the information from the observation data.

Hydraulic Conductivity

Base realization hydraulic conductivity arrays, zones, and pilot points are shown in figures 5–12. Horizonal hydraulic conductivity arrays from model layers 1–4 are shown in figures 5–8, and vertical hydraulic conductivity arrays from model layers 1–4 are show in figures 9–12. Base realization hydraulic conductivities, aggregated by geologic zones, are summarized in table 3. Base realization horizonal hydraulic conductivity values were generally larger than vertical values, but in some areas of model layers 1 and 3, adjusted vertical hydraulic conductivities (Kv) were larger than horizonal conductivities (Kh). Prior Kv values were set equal to 10 percent of Kh values. This relation was not enforced during history matching, and Kv and Kh were adjusted independently. These areas were limited to relatively small areas in layer 1 (in the highly permeable [sand, outwash] and [or] generally unstructured [kame, till] deposits) where it is not unreasonable for Kv to be similar or slightly larger than Kh. An analysis of model anisotropy—the ratio of Kh to Kv—aggregated by lithologic zone is shown in figure 22. A simple sensitivity analysis was done to evaluate the relative importance of cells with anisotropy greater than 1 to overall model performance. A version of the base realization model was created where Kh was increased to equal Kv in all cells where anisotropy was greater than 1.0. This adjustment did not result in an appreciable change to model performance (less than 0.5 percent change in the objective function PHI), so the few cells with anisotropy greater than 1.0 were allowed to remain for the rest of the analysis.

Boundary Conditions

The base realization recharge array and parameterization is shown in figure 17. The base realization recharge resulted in slightly greater variability from Yager and others (2019) and had a slightly lower mean recharge rate of 14.1 in/yr across the study area. Well pumping rates were adjusted during calibration, but base realization pumping rates compare closely to reported values, with a mean difference of just 1.4 percent from initial pumping rates among all wells. Figure 23 shows the simulated base flow in the streams represented by the SFR package throughout the model. Negative values indicate flow into the streams from the groundwater-flow system (gaining streams), and positive values indicate flow from the streams into the groundwater-flow system from the streams (losing streams). Most of the model domain is characterized by gaining streams, although near the outlet of the Neversink Reservoir in the northwest corner of the model, in some of the steepest valley wall slopes, and at valley margins, losing streams are simulated. Some dry SFR reaches were simulated in upland areas and likely represent intermittent, low-order streams.

Water Balance, Water Table, and Streamflow

The main components of the water balance, as simulated by the MODFLOW 6 model, are presented in figure 24 for the base realization. As expected from the conceptual model of the system, the main inflow component is recharge and the main outflow component is discharge to streams represented by SFR. Smaller inflow components are recharge from streams represented by SFR and inflow at CHD boundaries. The smallest outflow component is the water removed through well pumping, followed by flow out through CHD boundaries. As required in the formulation of MODFLOW 6, the fluxes in and out of the model balance to zero.

The simulated water table, where hydrostatic pressure is equal to atmospheric pressure, is shown in figure 25 as contours that depict a muted representation of the land-surface elevation, stylized in shaded gray. Figure 25 also shows areas where the model is “flooded”—locations where the simulated water table is higher than the top of the model. Such simulated “flooding” is a common issue in valley-fill regions with substantial areas of thinly covered bedrock. In reality, flooded areas may correspond to wetlands where groundwater is discharging at the land surface. In early versions of the model, literature-based parameter values, particularly for the bedrock, resulted in extensive simulated flooding. Furthermore, some areas with inadequate SFR representation resulted in the model’s inability to convey water from low-lying areas into streams. With the refined stream network and parameters informed through history matching, the extent of flooding is limited to areas deep in the valleys and some highlands that are generally coincident with wetlands and ponds or lakes. Routing of water simulated above land surface could be implemented through the MODFLOW 6 unsaturated zone flow (UZF) package (Feinstein and others, 2020); however, we judged the extent of flooding to be sufficiently limited that the additional package, accompanied by more computational complexity, was not necessary.

Providing an Ensemble for Monte Carlo Analysis of Areas Contributing Recharge

The end result of the history matching effort was an ensemble of parameter values for the groundwater-flow model that spans the range of values within prior uncertainty bounds and that also is informed by the historical observation data, within a reasonable range of uncertainty on those observations. We carried this ensemble through to the simulations to delineate the areas contributing recharge by running the groundwater-flow and particle tracking models for each ensemble member, as we discuss in subsequent sections.

Probabilistic Areas Contributing Recharge

The areas contributing recharge to the priority production wells were delineated using the groundwater-flow model and a particle-tracking simulation. Parameter sets produced by the history matching process were used to run MODFLOW. The unique flow field generated by each parameter realization and MODFLOW were then used to build MODPATH7 (Pollock, 2016) particle tracking input files. For the MODPATH7 simulations, one particle was placed in each cell in the uppermost active layer (in the same way that recharge is applied to the model) at the beginning of the simulation. Particle tracking was then run in the forward direction. The movement of the particles simulate the movement of recharge through the aquifer system, moving downgradient toward sinks (such as wells) and other model boundaries. MODPATH7 tracks the movement of the particles during the simulation and records where they terminate at model sinks. In this study, we instructed particles to stop when they met weak sink cells, which was necessary because some of the priority wells in the model are weak sinks (see Abrams and others, 2013 for additional information about weak sinks). At the conclusion of the particle-tracking simulation, the particles that terminated in priority well sinks were traced back to their starting locations and mapped to the model grid. This produced a single, deterministic map of the model cells that do (and do not) contribute recharge to priority wells.

Prediction Uncertainty Analysis

Results of the history-matching process produced an ensemble of 500 groundwater model parameter realizations that account for both model parameter and calibration data uncertainty. Of these 500 parameter realizations, 142 were trimmed from the ensemble during rejection sampling analysis to remove anomalously high PHI values (fig. 20) and thereby parameter sets that result in unrealistic groundwater levels and streamflows. The remaining 358 parameter realizations were used to complete a Monte Carlo analysis to quantify the uncertainty of the contributing recharge area to priority water-supply wells. The MODPATH particle tracking process, described in the “Probabilistic Areas Contributing Recharge” section, was repeated for each of the 358 realizations in the history-matched parameter ensemble. Deterministic contributing recharge areas produced from these repeated runs were used to calculate the probability of a location being in the area contributing recharge to a priority well.

In the history matching part of the project, best estimates of mean water use were used for the well package in MODFLOW 6; however, these estimates are typically not at the highest range of capacity for the wells, in particular for the priority wells (table 1). As a result, variations in pumping rates related to population growth and other future decisions were added as an additional consideration beyond model parameters that contribute to the overall uncertainty of the contributing recharge area. The uncertainty in pumping rates was represented using six pumping scenarios:

The Monte Carlo analysis (358 simulations) was repeated for each of the 6 pumping scenarios (for a total of 2,148 simulations) to incorporate the uncertainty related to the effects of pumping on the areas contributing recharge to priority wells. The results of these analyses are shown in figures 26–31. The general shapes of the contributing recharge areas were consistent among pumping scenarios, with increasing pumping rates leading to larger contributing recharge areas and larger areas of high probability that a location contributes recharge to priority wells; for example, the footprints simulated during the “Maxratio” scenario (fig. 27) cover larger regions and have larger areas of high probability than those simulated during the “Baseratio” scenario (fig. 26).

Panel B of figures 26–31 shows the probabilistic contributing recharge area to the 4 priority wells near Ellenville, N.Y. The Ellenville contributing recharge area is the largest in size for all scenarios, owing the number (4) and density of wells and their combined maximum withdrawal rates (table 1). The Ellenville footprint extends to the south, to near Spring Glen, N.Y., and the upper Sandburg Creek drainage (fig. 1). The footprint also extends laterally to the east and west up the valley walls into the uplands, particularly to the east of Ellenville between the tributary streams draining Shawangunk Ridge. The highest probability areas of the Ellenville contributing recharge area were typically simulated along the eastern part of the valley floor and the eastern uplands, with some high probability areas also simulated in the western uplands.

Panel C of figures 26–31 shows the probabilistic contributing recharge areas to the 3 priority wells near Woodridge (2) and Mountain Dale (1), N.Y. Contributing recharge areas near the two western wells near Woodridge cover parts of the valley floor near the wells and, to a larger extent, areas in the uplands to the north, east, and west of the wells. The contributing recharge area for the eastern well, pumping from a confined aquifer unit near Mountain Dale, includes an area south of the well in the uplands and a small part of the valley floor near the well. During intermediate and higher pumping scenarios, additional contributing areas are simulated in the uplands to the north and south of the well (fig. 27, for example); however, these areas have relatively low probability of contributing recharge.

Panel D of figures 26–31 shows the probabilistic contributing recharge area to the priority well near Wurtsboro. N.Y. The contributing recharge area for this well includes a part of the valley floor near the well and the uplands to the west of Wurtsboro. The extent of this well’s contributing recharge footprint does not change greatly among pumping scenarios, but results indicate that increased pumping at this well leads to increased probability of recharge contribution.

Each pumping scenario represents an outcome that, when used for planning purposes, corresponds to a different level of risk tolerance; for example, pumping rates in the “Maxratio” scenario are scaled to maximum capacities, creating the largest contributing recharge footprint and large areas of high probability near wells. Management decisions made using results from this scenario (fig. 27) would be relatively risk averse. Alternatively, pumping rates in the “Baseratio” scenario are adjusted only by history matching and remain near 2010–2015 reported mean levels, creating contributing recharge footprints that are relatively small and smaller high probability areas. Management decisions made using “Baseratio” results (fig. 26) would be relatively risk tolerant. Management decisions made using and the intermediate scaled and ranged pumping scenarios (“0.25ratio,” “0.5ratio,” “0.75ratio,” and “Maxrange”) therefore represent intermate levels of risk tolerance. Taken together, these results offer a view of the overall uncertainty of the areas contributing recharge to priority wells in the NS–RO study area and provide a tool for risk-based decision making related to the source water protection of these wells.

Model Discretization

The model discretizes the groundwater-flow system into a simplified model grid, consisting of cells that are 50 m on a side and to more than 100 m in thickness. Values in these model cells therefore represent averages of the respective cell volumes. Similarly, simulated outputs from each cell represent average conditions within that cell.

Steady State Assumption

For the purposes of simulating contributing recharge areas to priority wells it was assumed that the steady state simulations were appropriate. In this model, recharge estimates from 2000 to 2013 overlap with water use estimates developed from the 2010–2015 pumping data. Perfect alignment of stresses and observations with the steady-state period is difficult to obtain which can affect the results. Furthermore, steady state conditions represent a time-averaged delineation of source water areas and can smooth over temporal variability (Rayne and others, 2014). However, the relatively high hydraulic conductivity in the glacial sediments—where the priority wells are screened—likely negates much of the transience that can cause elasticity of the source water areas.

Data Gaps

Observation data used for history matching were few and often of poor or unknown quality. This was especially true for hydraulic head observations, where all but one observation was from NYSDEC well completion reports. Flux observations were developed from BFI estimates, which are also uncertain. We adapted to this by including realizations of observation values consistent with their uncertainty for the history matching effort.

Only two flux calibration targets were available. Furthermore, the flux and hydraulic head datasets did not always overlap in time with the steady-state period. The groundwater-flow model simulated base flow only, so the targets are expected to be representative of the groundwater contribution to flow only. However, the presence of the Neversink Reservoir just upstream from streamgage 01436500 made it difficult to assess the groundwater contribution only. The uncertainty of the reservoir outlet importance at streamgage 01436500 was, in part, incorporated into the weight for that observation. Nonetheless, more streamflow data on this stream and others in the model would better constrain the water balance.

Particle Tracking

At base pumping levels, some of the priority wells represented in the groundwater-flow model are strong sinks and others are weak sinks. Strong sinks extract water from the entire cell depth, whereas weak sink well cells extract water from some of the cell’s depth but also let some water pass through the cell without being extracted (Abrams and others, 2013). Allowing particles to pass through weak sinks in this model could therefore lead to an underprediction of the areas contributing recharge to wells. To avoid underpredicting the size of the contributing recharge areas, the particle tracking analysis in this study stopped all particles when they reached weak sinks. This approach offers a conservative estimate of the areas that contribute recharge to the priority wells but may slightly overpredict the size of and the probability that these areas contribute recharge to the wells.

Each particle tracking Monte Carlo analysis consisted of 384 particle tracking simulations—one for each realization in the posterior parameter ensemble. In this study, it was assumed that the 384 realizations were sufficient to reach statistical stationarity in the probabilistic contributing recharge areas. An evaluation of the effect of ensemble size on contributing recharge area stationarity was completed by Fienen and others (2021). This evaluation indicated that higher probability areas tend to stabilize with fewer realizations than the lower probability areas and that the threshold for contributing recharge area stationarity depends on the user’s adopted level of risk tolerance. The probabilistic contributing recharge areas maps produced by this study (figs. 26–31) are considered generally representative but may slightly underpredict the extent of low probability contributing recharge areas. Use of a larger ensemble for the particle tracking Monte Carlo analysis would likely have resulted in similar probability maps—particularly in high probability areas—but with some additional low-probability areas. Further discussion on the effect of ensemble size to map stationarity is presented in Fienen and others (2021).