Hydrogeology

The Cortland study area of 50.2 square miles (mi²) includes the Otter Creek and Factory Brook Basins along with small parts of West Branch and East Branch Tioughnioga River Basins and of Tioughnioga River Basin (fig. 3). Otter Creek and Factory Brook are tributaries to the West Branch Tioughnioga River, which in turn joins with the East Branch Tioughnioga River to form the Tioughnioga River. In the headwaters of Otter Creek Basin, groundwater and surface-water divides do not coincide (Miller and others, 1998). Groundwater in the valley-fill deposits from an adjoining basin flows northeastward toward Otter Creek Basin, whereas groundwater inflow to the valley from an adjoining upland hillslope flows away from Otter Creek Basin (fig. 3).

Glacial and postglacial deposits (fig. 3) overlie shale and siltstone bedrock in the study area. The altitude of the bedrock surface beneath these deposits is shown in figure 4; the bedrock surface rises 1,090 ft from the valleys to the highest areas in the uplands. A discontinuous layer of till deposited directly on the bedrock by glacial ice consists of an unsorted mixture of sediments ranging in size from clay to boulders.

In the valleys, an unconfined aquifer consisting of mostly sand and gravel outwash generally overlies a confining unit of lacustrine very fine sand, silt, and clay deposits, which in turn overlie a confined to semiconfined aquifer of sand and gravel basal deposits of mostly ice-contact kame deposits (Miller and Brooks, 1981; Miller and others, 1998). In some areas near the valley walls, surficial kame deposits directly overlie the lower sand and gravel unit. The vertical distribution and thickness of these surficial deposits along two generalized west–east geologic sections through Otter Creek Valley and Factory Brook Valley, drawn based on the surficial mapping and geologic sections in Miller and Brooks (1981) and Miller and others (1998), are shown in figure 5. In Otter Creek valley, the unconfined sand and gravel of mostly outwash ranges from 40 to 80 ft thick. The lacustrine very fine sand, silt, and clay deposits, where present, range from about 170 ft thick in the west part to 90 ft thick in the east part of Otter Creek valley. Buried kame deposits are generally 60 to 170 ft thick in Otter Creek valley but thinner (2 to 30 ft) north of Otter Creek valley in the West Branch Tioughnioga River valley. In Factory Brook valley, the unconfined aquifer is about 100 ft thick but there was limited lithologic information to determine the thickness of the lacustrine unit and confined aquifer. A detailed description of the hydrogeology of the study area, including its glacial depositional history, is available in Miller and Brooks (1981) and Miller and others (1998). The digital hydrogeologic maps developed as part of the present study that serve as the basis for the framework of the groundwater-flow model are available as a USGS data release (Woda and others, 2022). Data sources and the methods used to develop the hydrogeologic maps are described in Finkelstein and others (2022).

Numerous aquifer tests have been done in Otter Creek Basin to estimate horizontal hydraulic conductivity of the principal glacial deposits (Miller and others, 1998). Aquifer-pump tests, completed in the well-sorted sand and gravel outwash deposited by meltwater typically distal to the ice front, indicated hydraulic conductivity values ranging from 880 to 1,150 feet per day (ft/d) in the western and central parts of the valley and lower values of 220 to 380 ft/d in the eastern part of the valley and adjacent to Tioughnioga River. Analysis of slugtests in mostly outwash mixed with some alluvial sediments from upland hillslope runoff near the valley wall indicated hydraulic conductivity values ranging from 3 to 140 ft/d. Hydraulic conductivity values from aquifer-pump test results in kame deposits in the unconfined and confined aquifer ranged from 60 to 85 ft/d; these ice-contact deposits can range from poorly sorted to well-sorted sand and gravel. Hydraulic conductivity values of 2 ft/d from slugtest and grainsize analyses were estimated for lacustrine deposits of very fine sand and silt that were deposited in proglacial lakes.

For postglacial deposits (alluvium and streambed deposits) and shale and siltstone bedrock, measurements of hydraulic conductivity and transmissivity have been made in other areas of the Susquehanna River Basin. Randall (1978) calculated hydraulic conductivity of alluvium described as silty gravel at seven field sites that ranged from 13 to 134 ft/d and averaged 43 ft/d. Yager (1986) estimated hydraulic conductivity of coarse-grained bed sediments of the Susquehanna River as 1 and 6 ft/d. For the shale and siltstone bedrock, Randall and Mills (2020) calculated a median transmissivity of 42 feet squared per day (ft²/d) from short-term pump tests.

Recharge in upland till and bedrock areas such as in the Cortland study area is primarily from direct infiltration of precipitation but may also include leakage from streams. Recharge rates in upland settings are conceptually highly variable, ranging from near zero in low-permeability tills on steep topography, where the water table is near land surface, to values approaching the maximum water available for recharge (precipitation minus evaporation) in high-permeability tills on moderate slopes, where the water table is perennially below the land surface.

Sources of recharge to the valley-fill aquifers include direct infiltration from precipitation, runoff from adjacent upland hillslopes, and natural infiltration from streams that cross a valley from upland areas. In some cases, pumping by wells may also induce water from the streams. In Otter Creek Basin, streams that drain the uplands lose water as they flow across the transmissive valley-fill deposits, especially near the upland-valley contact (Miller and others, 1998).

Development of Groundwater-Flow Model

Groundwater levels and flows were simulated in this valley-fill and upland setting that encompasses 50.7 mi² in west-central Cortland County. The groundwater-flow model was calibrated based on 69 groundwater-level and streamflow (base flow) observations, most of which were measured May 28 to June 4, 1991, during average annual hydrologic conditions (Miller and others, 1998). It was assumed that long-term average annual conditions would be similar if more recent climatic conditions were also considered. Well withdrawals during this 1991 calibration period were simulated in the model and then replaced with 2010–15 average withdrawal rates to delineate the contributing areas to the production wells that supply drinking water. Most boundary conditions and hydraulic properties were represented by parameters (table 2) for calibration by nonlinear regression and for evaluating model-prediction uncertainty.

Hydrologic Boundaries

Hydrologic boundaries include the movement of water into and out of the groundwater-flow model. These hydrologic boundaries include recharge from precipitation, the interaction between streams and aquifers, and well withdrawals and associated return flows.

The spatial distribution of recharge rates was applied based on the results of a soil-water-balance model. Yager and others (2018) calculated mean recharge rates for 2000–13 for the glaciated conterminous United States at a resolution of 250 by 250 meters. The soil-water-balance method uses spatially distributed variables such as climate, soil type, land use, and soil water capacity (Westenbroek and others, 2010). Mean annual recharge rates for the 2000–13 period ranged spatially from 8.9 to 28.2 inches per year (in/yr) in the Cortland study area (fig. 6). Lower rates were calculated for upland deposits with low infiltration capacity and developed land cover, and higher rates were calculated for valley-fill deposits with high infiltration capacity and undeveloped land cover. The mean annual recharge rate averaged spatially over the modeled area was 13.4 in/yr. A dimensionless parameter R_Mult (table 2) with a value of 1 multiplies the spatially varying recharge rates. This parameter was not considered for calibration by nonlinear regression; instead, the uncertainty in these recharge rates was included in the uncertainty analysis of the simulated contributing areas to the wells.

Stream-aquifer interactions were simulated as a head-dependent flux boundary in layer 1 (fig. 7) by using the streamflow-routing package (Niswonger and Prudic, 2005) developed for MODFLOW. The streamflow-routing package accounts for gains and losses of water in each stream cell and routes streamflow from upstream cells to downstream cells. Stream altitudes were determined from light detection and ranging (lidar) imagery (Cortland County, New York, 2017). Each streamflow-routing cell requires a conductance term that incorporates the geometry and vertical hydraulic conductivity of the streambed. Water depths and bed thicknesses of 1 ft were used to determine the top and bottom bed altitudes from surface-water altitudes. Simulated streams were assigned widths ranging from 5 to 10 ft for small tributaries, 20 ft for some stream reaches of Factory Brook and Otter Creek, 50 ft for West Branch and East Branch Tioughnioga Rivers, and 100 ft for Tioughnioga River. A parameter, SB_K (table 2), defined the vertical hydraulic conductivity of the streambed materials. Base flows entering the model from areas not directly simulated (West Branch and East Branch Tioughnioga Rivers) were specified at the first boundary stream cell. Streamflows were specified based on base-flow statistics (base-flow index; average proportion of streamflow that is base flow) from the Tioughnioga River at Cortland, New York, streamgage (station number 01509000) (USGS, 2020a). The area-weighted average of base flow from this long-term streamgage was used to calculate base flow at the West Branch and East Branch Tioughnioga Rivers based on their respective drainage area sizes. The model contained a total of 4,286 streamflow-routing cells.

Well withdrawals during the 1991 calibration period in or near the Otter Creek Basin were from six high-capacity production wells (Miller and others, 1998). For well withdrawals in other areas of the modeled area, it was assumed that 2010–15 average well withdrawals from the Village of Homer (two production wells) and Village of Hamilton (two production wells) adequately represented pumping rates from these wells during 1991. A total of 10 wells totaling 4,752 gallons per minute (gal/min), were simulated in the model during the calibration period (table 1). Withdrawals from small-capacity domestic wells were not included in the model simulations. Most pumped water from domestic wells is returned to the aquifer through nearby onsite septic systems with little net change in flow.

Water withdrawn by the City of Cortland, Town of Cortlandville, and the Village of Homer was used in sewered areas and then treated at a wastewater facility before discharging to the Tioughnioga River (City of Cortland, New York, 2020). One of the high-capacity wells was an extraction well at a hazardous waste site; after onsite treatment, all water pumped from this extraction well was returned to the aquifer through recharge basins near the well (Miller and others, 1998). Infiltration of this pumped water was simulated in 10 model cells. For the remaining wells in the City of Cortland, it was assumed that pumped water was either treated at the wastewater facility or that by not including return flow it would have negligible effects on model calibration results.

Model Calibration

The steady-state groundwater-flow model was calibrated with the inverse modeling program UCODE–2014 (Poeter and others, 2014; Hill and Tiedeman, 2007) using nonlinear regression that minimizes the differences, or residuals, between field (observed) and simulated water levels and base flows to obtain an optimal set of parameter values. The quality of this calibration was determined by the reasonableness of the estimated parameter values and by analysis of the residuals. Some parameters might be insensitive to the available observations and therefore cannot be estimated by nonlinear regression. Values from the literature, primarily from Miller and others (1998), were used to specify parameter values that could not be estimated by nonlinear regression.

Parameter values in the groundwater-flow model were estimated using 67 groundwater-level observations and 2 base-flow observations (fig. 7). Observations were weighted to account for the difference in the type of observations and their relative effect in nonlinear regression. Observation weights are equal to the inverse of the variance (square of the standard deviation) of the measurement and model structure error. Groundwater-level observations were divided into two groups. Group A included 49 water levels measured between May 28 and June 4, 1991, at long-term average conditions (Miller and others, 1998) with surveyed measuring-point altitudes. These water levels were measured in valley-fill deposits mostly in outwash in layers 1 and 2 in and near Otter Creek Basin. Group B included 18 water levels measured between 1931 and 1991 at various hydrologic conditions (USGS, 2019) with measuring-point altitudes estimated from land-surface altitudes determined from USGS topographical contours. These water levels were also measured in valley-fill deposits mostly in outwash in layers 1 and 2 generally throughout the modeled area. The standard deviations used in model calibration ranged from 0.47 for Group A to 1 for Group B.

Base-flow measurements during the calibration period (May 28 to June 4, 1991) at four partial-record sites, two on Otter Creek and two on Dry Creek, were used for two base-flow observations. These two observations represented the net base-flow loss between the partial-record sites (1.1 ft³/s for Otter Creek and 0.8 ft³/s for Dry Creek). The accuracy of the base-flow measurement at each partial-record site was estimated to be 10 percent. The coefficient of variation (standard deviation divided by the mean) from which the observation weight was calculated was 11 percent for the base-flow observation on Otter Creek and 16 percent for the observation on Dry Creek.

Twenty-one model parameters were evaluated with parameter estimation: 10 for horizontal hydraulic conductivity (HK_Alluvial, HK_Outwash1, HK_Outwash2, HK_Outwash3, HK_Lac, HK_Basal, HK_Kame, HK_Till-Mor, HK_Till, and HK_Rock), 10 for their associated vertical anisotropy (VANI), and 1 for streambed vertical hydraulic conductivity (SB_K) (table 2). Because the recharge estimates used in the model were derived independently, they were not adjusted during calibration. Parameter sensitivities indicated the available groundwater-level and base-flow observations provided sufficient information to permit an estimate of three parameters (HK_Outwash1, HK_Outwash2, and SB_K). Most of the remaining 18 parameters were specified based on manual calibration and values in the literature (see section “Hydrogeology” for Cortland study area) (table 2). In the uplands, final values determined by manual calibration for horizontal hydraulic conductivity for upland till (HK_Till) and shale and siltstone bedrock (HK_Rock) were 3 and 0.1 ft/d, respectively, so as to simulate the water table approximately between the land surface and the bedrock surface. There were no literature values for vertical anisotropy in central New York but a ratio of 10:1 for all surficial deposits and 1:1 for bedrock were assumed to adequately represent these materials.

Optimal values for two horizontal hydraulic conductivity parameters (HK_Outwash1, HK_Outwash2) and for streambed vertical hydraulic conductivity (SB_K) that were estimated by nonlinear regression are within a plausible range of values (table 2). Values for the two outwash parameters, HK_Outwash1 (555 ft/d) and HK_Outwash2 (543 ft/d), however, are nearly the same although aquifer tests (Miller and others, 1998) indicated that HK_Outwash1 might be much greater than HK_Outwash2. Although there were no direct measurements of streambed vertical conductivity in the study area, the SB_K value of 1.7 ft/d was considered reasonable based on the average value used by Miller and others (1998) of 3.5 ft/d determined from manual calibration, by field measurements of 1 and 6 ft/d by Yager (1986) in the Susquehanna River, and because of the number of terms required to define streambed geometry.

The uncertainties of the parameter estimates, indicated by the 95-percent linear confidence intervals (table 2), are within the ranges of reasonable values reported in the literature. A comparison of the relative precision of different parameter estimates can be made using the coefficient of variation (standard deviation of the estimated value divided by the optimal value; table 2); a smaller coefficient of variance indicates a more precisely estimated value for the parameter. The coefficient of variations ranged from 0.05 to 0.20. Parameter HK_Outwash1 was the most precisely estimated, whereas parameter SB_K was the least precisely estimated.

The quality of model calibration can be determined, both numerically and graphically, by comparison of the observations and simulated equivalents. Weighted residuals should be randomly distributed and close to zero. The average weighted residual was −0.78 ft for all groundwater-level and base-flow observations, and it ranged from a minimum of –11.58 ft to a maximum of 12.07 ft. The sum of squared weighted residuals was 188 for the calibrated model.

The relation between weighted observed values and weighted simulated values is shown in figure 8A, which symbolizes the weighted head residuals by the category of their groundwater-level group. Most values plot near a line with a 1:1 slope; the correlation between them was 0.99. Weighted water-level residuals are generally randomly distributed around zero for all unweighted simulated values (fig. 8B).

The spatial distribution of weighted base-flow and groundwater-level residuals is shown in figure 9. Unweighted residuals for groundwater levels would be one-half of the weighted residuals for Group A and the same for Group B observations (fig. 7). Weighted residuals are generally distributed randomly throughout most of the modeled area.

The two observations of base-flow loss between partial-record sites on Otter Creek and Dry Creek compare favorably with simulated values. The values of simulated net base-flow loss for Otter Creek (1.3 ft³/s) and for Dry Creek (1.1 ft³/s) were close to observed values of 1.1 and 0.8 ft³/s, respectively.

The altitude and configuration of the simulated water table for the steady-state calibrated model are shown in figure 9 at 100-ft contour intervals except for 50-ft intervals in the valleys and near the valley-upland contacts. These simulated water-table contours are consistent with the conceptual model of groundwater flow in the study area. Groundwater generally flows from topographically high areas and discharges to streams and surface-water bodies. In the uplands, the simulated water table approximately parallels the land surface, and simulated groundwater divides generally coincide with watershed divides. The water-table gradient is steepest in the till and bedrock uplands and in valley-fill deposits near the contact and then it flattens in the more transmissive valley-fill areas of stratified deposits.

The simulated groundwater budget for the calibrated model indicated that recharge from direct precipitation provides 51.2 cubic feet per second (ft³/s), which constitutes 74 percent of the inflow. Streamflow loss accounts for the remaining inflow, 18.2 ft³/s (26 percent). Most of this streamflow loss is from tributaries downstream from the upland-valley contact, which are near areas of abrupt changes in transmissivity and near large-capacity pumping wells. Of the total inflow, 58.9 ft³/s (85 percent) of the groundwater discharges to streams, and 10.6 ft³/s (15 percent) is withdrawn by the pumping wells during the calibration period. The percentage of streamflow loss as part of the total inflow to the aquifer is within the range of groundwater model simulations in the Otter Creek Valley (32 to 38 percent) (Miller and others, 1998) and in the Tioughnioga River Valley downstream from this study (19 to 23 percent) (Miller, 2000).

Simulation results of the calibrated groundwater model indicated that the model is acceptable for use in the study to delineate steady-state areas contributing recharge to drinking-water supply wells. Model-fit statistics indicated that simulated values for groundwater levels and for base flows are generally close to observed values. Optimal parameter values are realistic, and their confidence intervals include reasonable values. The simulated water table and water budget are consistent with the conceptual understanding of valley-fill systems bordered by till uplands.

Areas Contributing Recharge and Prediction Uncertainty Analysis

Calibration of the groundwater-flow model by inverse modeling by using nonlinear regression provided an optimal set of parameter values. This optimal parameter set was estimated by minimizing the weighted residuals between the observation dataset (67 groundwater levels and 2 base flows) and simulated values. A predicted area contributing recharge to the production wells based on this optimal parameter set provides a single, most likely contributing area (deterministic contributing area) based on the current understanding of the geologic and hydrologic features of the modeled area. The parameter values, however, were estimated with different levels of uncertainty; this uncertainty in the optimal values was based on the information that the observation dataset provided on the parameters. Parameter uncertainty and its associated effects on model predictions (spatial variability of the simulated contributing area to a well) can be evaluated by a stochastic Monte Carlo analysis. The parameter variance-covariance matrix from nonlinear regression can be used to create plausible parameter sets for the Monte Carlo analysis (Starn and others, 2010). The parameter variance-covariance matrix incorporates the uncertainty of the parameter estimates and the correlation among parameters from the calibrated model. The Monte Carlo analysis was done by replacing the parameter set in the calibrated model by a plausible parameter set multiple times. The probability of a location being in the contributing area to a production well was calculated from these multiple model simulations.

Areas contributing recharge were determined for the 10 production wells that supply drinking water based on their 2010–15 average withdrawal rates (table 1). The two Village of Hamilton wells had very low average withdrawal rates (3 and 5 gal/min). At this low withdrawal rate only a small part of the water that flows into the model cell with the well screen is removed by the production well (weak sink). Because it cannot be determined whether water that enters this cell is withdrawn or continues through the groundwater system, simulations of particles for these two wells were stopped at these two weak sinks (instead of allowed to pass through weak sinks). Thus, the contributing areas to these two wells are larger, or more conservative, than if it was not a weak sink cell. In addition, the Village of Hamilton wells were not incorporated into the uncertainty analysis of the contributing areas because of their low pumping rates and the fact the simulated contributing areas to these wells were considered conservative.

Areas Contributing Recharge

Simulated deterministic areas contributing recharge and simulated groundwater traveltimes to the production wells were determined on the basis of the calibrated steady-state model for simulated pumping conditions and tracking of pathlines with the MODPATH particle-tracking program. Simulated areas contributing recharge for the 2010–15 average withdrawal rates extend upgradient from the wells to groundwater divides and topographical divides, some of which serve as model boundaries (fig. 10). Most areas contributing recharge extend beneath and beyond streams and include upland till deposits. Areas contributing recharge include isolated areas remote from the wells including small areas the size of a model cell. In addition, most wells are not overlain by their contributing area, including Cortlandville well 2 and well 7, which are overlain by the contributing area to the downgradient Cortland well field. The wells in the study area are screened in the deep part of the outwash deposits (layer 2), and recharge near the wells travels along flow paths above and around the screened interval of the well.

The size of the simulated areas contributing recharge to the Cortland and Cortlandville wells in Otter Creek Basin covered a total area of 5.86 mi². For the large total pumping rates for these wells, the contributing areas covered most of the valley southwest of the Cortland well field where there is substantial amount of urban and agriculture land uses (National Land Cover Database, 2016). Land cover in the contributing area to the Cortland well field and to Cortlandville Terrace Road well have urban land use considered high density. In the uplands, the land cover in the contributing areas to these wells is primarily agriculture and forested. The size of the contributing area to Cortlandville Terrace Road well is relatively small for its pumping rate; this well, which is situated downgradient of a long losing reach of Otter Creek that parallels the valley-upland contact, captures a substantial amount of this streamflow loss. The contributing areas to Cortlandville well 2 and well 7 are in the least developed part of the valley but include some agriculture and urban land uses. A nearby sand and gravel mining operation southeast of Cortlandville well 2 and well 7 is within the simulated contributing area to the Cortland well field.

In Factory Brook Basin, the size of the simulated areas contributing recharge to the Homer well field covered a total area of 0.99 mi². This well field, which is screened in the deep part of the outwash deposits, is on the west edge of the Village of Homer but the land cover in its contributing area, which is distal to the well field, is primarily agriculture and forested (National Land Cover Database, 2016). However, small parts of a sand and gravel mining operation, shown in World Imagery (2018) along the southwest side of Factory Brook valley, were in simulated parts of the contributing area to the well field (fig. 10). Recharging water between the well field and its contributing area travels along shallow and intermediate depth flow paths before discharging to streams. The size of the simulated contributing area to the Hamilton wells covered 0.08 mi². Land cover in this contributing area is agriculture and forested.

Simulated traveltime estimates from recharging locations to the production wells for the 2010–15 pumping rates are shown in figure 11. Porosity values (table 2) were specified in the model for MODPATH, but they affect only groundwater velocity and do not change the contributing areas to the wells in a steady-state simulation. Glacial deposits and postglacial deposits were assigned a porosity of 0.35 based on values determined for similar deposits in Rhode Island (Allen and others, 1963), Massachusetts (Garabedian and others, 1991), and southern New England (Melvin and others, 1992). For bedrock, a porosity of 0.05 was assigned based on values for shale bedrock (Freeze and Cherry, 1979).

Traveltimes generally depend on where recharge enters the aquifer in relation to the production wells. Water that recharges the aquifer near the wells has the shortest traveltimes and youngest water, whereas water originating in the till uplands has the longest traveltimes and the oldest water. Traveltimes ranged from less than 1 year to more than 300 years; 50 percent of the traveltimes were 10 years or less. The median traveltime was 10.2 years for the wells in Otter Creek Basin and 8.7 years for the Homer well field. For the Hamilton well field, which has very low pumpage and weak sinks, the median traveltime was 20.5 years. This percentage for traveltimes 10 years or less and the relatively short median traveltimes indicate that the wells are vulnerable to contamination from activities on the land surface.

Probabilistic Areas Contributing Recharge

The effects of parameter uncertainty on model predictions (the predicted contributing area) were quantitatively measured by a Monte Carlo analysis. A Monte Carlo analysis was used to obtain the probability of a recharge location being in the contributing area of the well fields. The probability distribution is related to the information that the observation dataset provided on the estimated parameters, to prior information on specified parameters, and to the sensitivity of the simulated contributing area to the parameters. Hundreds of parameter sets generated from summary statistics of the calibrated steady-state model were used to run hundreds of model simulations in the Monte Carlo analysis. Because combinations of reasonable parameter values might result in unrealistic groundwater levels and streamflows, the parameter sets were evaluated using the pumping rates and associated observation dataset from the calibrated model. Those parameter sets that simulated realistic results were then used in Monte Carlo analysis for the 2010–15 average pumping rate scenario.

Parameter values for the Monte Carlo analysis were created from the parameter values in the model and its parameter variance-covariance matrix (Starn and others, 2010; Starn and Bagtzoglou, 2012). Parameter values that could not be estimated by nonlinear regression and thus were not included in the parameter variance-covariance matrix of the calibrated model might still be important for model predictions (for this study, the size, shape, and location of the area contributing recharge to the production wells). Horizontal hydraulic conductivity parameters of the principal glacial deposits (HK_Outwash3, HK_Lac, HK_Basal, HK_Kame, and HK_Till) were incorporated into the parameter variance-covariance matrix using their specified parameter values. Recharge rates, although calculated independently from the model, were also included in the parameter variance-covariance matrix because of the importance of this boundary condition on the simulated contributing areas to the wells. Vertical anisotropy of the surficial sediments and bedrock was not included in this analysis because they have little to no effect on the simulated contributing area to a well in a valley-fill setting for a plausible range of values (Friesz, 2010). Incorporating these five horizontal hydraulic parameters and one recharge parameter into the parameter variance-covariance matrix, however, caused large unrealistic uncertainties around the specified parameter values because the information the observations provided was insufficient. Prior information on these parameters was used to constrain this uncertainty to include most plausible values (table 2). Parameter uncertainties are from the observation dataset, but also from prior information on parameters that the modeler provided.

For the Monte Carlo analysis, the model was first run with 500 parameter sets and with the 1991 pumping rates used in calibrating the model. The nine hydraulic and recharge parameter values in each dataset replaced the corresponding parameter values in the calibrated model. Three criteria for accepting a given parameter set were used: (1) the model converged, (2) the model mass balance was 1 percent or less, and (3) a model-fit statistic (calculated error variance, which is the sum of squared weighted residuals divided by the difference between the number of observations and the number of parameters estimated by nonlinear regression) was less than a specified value. The third acceptance criterion was used so that model-prediction uncertainty would not be overestimated by using a parameter set that produced unrealistic groundwater levels or streamflows compared with that for the calibrated model. The value used for this criterion, however, can be model dependent and subjective. For this model application of the Monte Carlo analysis, a calculated error variance of 7 was selected for the third criterion, or about 2.5 times the calculated error variance of the calibrated model. Of the 500 parameter sets run with MODFLOW, 470 sets (94 percent) fit these three criteria.

A Monte Carlo analysis was then done using the parameter sets that fit the acceptance criteria for the 1991 average withdrawal rates but using the 2010–15 average withdrawal rates. The criteria for the Monte Carlo analysis that used the 470 parameter sets and these pumping rates were that the water budget error be 1 percent or less for models that converged. For the average withdrawal rates from 2010 to 2015, 467 parameter sets (99 percent) fit these criteria and thus were run with the particle-tracking program. The probability that a recharge location would be in the area contributing recharge to the production wells was determined by dividing the number of times a particle at a given location was captured by a well by the total number of accepted particle-tracking simulations; this probability was expressed as a percentage.

The probabilistic areas contributing recharge to the production wells in Otter Creek Basin and for the Homer well field for the average pumping rates from 2010 to 2015 are shown in figure 12. For this analysis, the probabilistic distribution for the closely spaced Cortland and Cortlandville well fields in Otter Creek Basin was done for the wells together. Probabilistic contributing areas to individual well fields or wells might overlap, even though the deterministic contributing areas do not overlap under a steady-state simulation. The total size of the probabilistic contributing area for the Otter Creek Basin wells and for the Homer wells was larger than the deterministic contributing areas because of the effects of parameter uncertainty. The probabilistic contributing area analysis indicated that some areas not in the deterministic contributing area, including additional areas of urban and agricultural land use, might actually be in the contributing area.

In Otter Creek Basin, in most cases, areas with high probabilities (greater than 50 percent) generally coincide with the deterministic contributing areas for the Cortland and Cortlandville wells. The probabilistic contributing areas to the wells indicate additional areas in the valley and uplands might be the contributing area. Dry Creek Basin, an adjoining basin, which was not in the deterministic contributing area, might have areas that are in the contributing area to the wells. Low probabilities are generally in areas distant from the pumping wells and in areas where simulated streams intercepted precipitation recharge in the deterministic model. Additional areas that might be in the contributing area include recharge originating near the pumping wells that have relatively short groundwater-flow paths and traveltimes. In some cases, small isolated areas remote from the pumping wells that are in the deterministic contributing area are associated with low probabilities.

For the Homer well field, areas with probabilities of greater than 25 percent generally coincide with the deterministic contributing area. The probabilistic contributing area to this well field has a large spread in low probabilities, in part because in this area of the model most parameters could not be estimated by nonlinear regression. Instead these parameters were added to the uncertainty analysis with prior information, which included most plausible values. Streamflow measurements and additional groundwater levels in Factory Brook Basin might help reduce the uncertainty in the simulated contributing area to the well field by increasing the number of parameters that could be estimated by nonlinear regression. The probabilistic contributing areas indicated that precipitation recharge that discharges to Factory Brook and its tributaries in the deterministic model might instead go directly to the well field. These areas include probabilities of greater than 25 to 50 percent between the wells and the deterministic contributing area with relatively short groundwater-flow paths and traveltimes. That part of the sand and gravel mining operation along the southwest side of Factory Brook valley that is not in the deterministic contributing area is in the probabilistic contributing area associated with low probabilities of 10 percent or less (fig. 12).

Hydrogeology

The Chenango River drains a northeast–southwest trending valley-fill system before it eventually joins the Susquehanna River. The surface of the shale and siltstone bedrock rises steeply from the valley by as much as 980 ft to the topographical divide in the uplands that defines the study area extent in the north (fig. 14). Glacial till mantles the shale and siltstone bedrock surface in most places except in central areas of the valleys but can range from several feet to tens of feet in the uplands (Hetcher-Aquila and Miller, 2005).

In the valleys, glacial deposits (outwash and ice contact of mostly kame deposits) and postglacial deposits (mostly alluvium) form an unconfined sand and gravel aquifer (Hetcher-Aquila and Miller, 2005) (fig. 13). Beneath this unconfined aquifer, along most of the northwest side of the main Chenango River valley and in central areas parallel to the main valley axis, lacustrine deposits of very fine sand, silt, and clay overlie a confined aquifer of ice-contact deposits that are poorly to well-sorted sand and gravel. The three production wells that supply potable water are screened in this confined aquifer in the center of the valley. The vertical distribution and thickness of the unconsolidated deposits, along a generalized west–east geologic section through Chenango River valley near the production wells and drawn based on the surficial mapping and geologic sections in Hetcher-Aquila and Miller (2005), are shown in figure 15. Along the southeast side of the main valley, the unconfined aquifer consists of ice-contact kame deposits extending vertically to the bedrock surface and these deposits are connected to the confined aquifer. The unconfined aquifer is generally less than 20 ft thick along the northwest and central areas of the main valley but thicker along the southeast side of the valley where kame deposits overlie bedrock. The lacustrine confining unit is as much as 160 ft thick. The confined aquifer is generally 10 to 30 ft thick in the central parts of the valley but along the northwest side of the valley buried ice-contact deposits are less than 10 ft thick or not present. A detailed description of the hydrogeology of the study area, including its glacial depositional history, is available in Hetcher-Aquila and Miller (2005). The digital hydrogeologic maps developed as part of the present study that serve as the basis for the framework of the groundwater-flow model are available from Woda and others (2022). Data sources and the methods used to develop the hydrogeologic maps are described in Finkelstein and others (2022).

Sources of water to the confined or semiconfined aquifer, and thus the production wells, are an important factor in the location of their contributing areas. The fine-grained lacustrine deposits impede vertical flow to the confined aquifer from surface sources such as direct infiltration of precipitation and streamflow loss. Under pumping conditions, the increased vertical hydraulic gradient between the unconfined and confined aquifers may increase leakage through the confining unit, especially in areas of thin lacustrine deposits. Recharge to the unconfined kame deposits east of the wells likely is an important source of water because these deposits are connected to the confined aquifer. Westward and northwestward of the wells, thin ice-contact deposits may extend up along the valley wall to the land surface where recharge can also enter the confined aquifer. Another source of water to the confined aquifer may be groundwater inflow from adjacent uplands through till and shallow bedrock.

Development of Groundwater-Flow Model

Groundwater flow in the unconsolidated deposits and the underlying shale and siltstone bedrock was simulated by a five-layer model with a uniformly spaced grid. The groundwater-flow model was manually calibrated to 11 groundwater-level altitudes and to a generally shallow water table that approximates the land-surface configuration in the till and bedrock uplands. The 11 groundwater levels were measured from 1959 to 1992 at various hydrologic conditions (Hetcher-Aquila and Miller, 2005). Five of these measurements were from wells screened in the unconfined aquifer and the remaining six measurements from wells screened in the deep part of the valley in the confined aquifer. The altitudes of the measuring points were based on the altitude of the land surface determined from lidar imagery (Terrapoint USA, 2008). Because of the limited quantity of available groundwater levels and no streamflow measurements, model calibration was considered a rough estimate of hydraulic properties. Model characteristics, including hydraulic properties and recharge rates, are summarized in table 4.

The groundwater-flow model extended to natural hydrologic boundaries beyond the likely areas contributing recharge to the production wells. This extent also minimizes the effects of boundaries on simulated heads near groundwater divides separating Chenango River from Wheeler Brook east of the production wells. The lateral model boundaries included Wheeler Brook on the east; elsewhere, the edge of the model extended to presumed groundwater divides in the low permeable till uplands in most areas. The active model grid represented an area of 32.8 mi², consisted of 505 rows and 192 columns, and included a total of 292,340 cells each with 125 ft on a side.

Vertical discretization of the five-layered model was based on lithology (fig. 15) from hydrogeologic mapping by Hetcher-Aquila and Miller (2005) and updated by Woda and others (2022). The unconsolidated sediments were subdivided vertically into four model layers. In the valley-fill deposits the top two layers (layers 1 and 2) represented the unconfined aquifer, which consisted mostly of alluvium, outwash, and kame deposits. Layer 3 in the valleys represented lacustrine and kame deposits. Layer 4 in the valleys represented the confined aquifer, which consists of kame deposits. The three production wells were simulated in layer 4. In the uplands, till deposits were represented in all four layers; shallow bedrock areas less than 13 ft from the land surface were simulated as till deposits. The bottom layer (layer 5) represented bedrock with a constant thickness of 100 ft.

The model of the Greene study area specified the same types of boundary conditions that the Cortland study area did to represent sources of recharge and areas of discharge. Recharge rates were applied based on the results of a soil-water-balance model (Yager and others, 2018). Mean annual recharge rates calculated by the soil-water-balance model for 2000–13 ranged spatially from 8.3 to 28.0 in/yr in the Greene study area (fig. 16). Lower recharge rates were calculated for upland till deposits with low infiltration capacity and higher rates were calculated for valley-fill deposits with high infiltration capacity. The mean annual recharge rate, averaged spatially over the modeled area, was 11.7 in/yr.

Interactions between surface water and groundwater were simulated as a head-dependent flux boundary in layer 1 (fig. 17) by using the streamflow-routing package (Niswonger and Prudic, 2005). Stream altitudes were determined from lidar imagery. Water depths and bed thicknesses of 1 ft were used to determine the top and bottom bed altitudes from surface-water altitudes. Simulated streams were assigned widths ranging from 10 ft for small tributaries to as much as 30 ft for tributaries to Chenango River where they cross the valley, and 150 ft for the Chenango River itself. A vertical hydraulic conductivity of 3 ft/d was used to represent bed sediments. Base flow entering the model from areas not directly simulated was specified at the first boundary stream cell. Streamflows were specified based on base-flow statistics (base-flow index) at the Chenango River at Greene, New York, streamgage (station number 01507000) (USGS, 2020b). The model contained a total of 3,470 stream cells.

Average annual withdrawal rates for 2010–15 from the three Village of Greene wells (table 3), a total of 169 gal/min, were used to represent long-term average rates. The three production wells are screened in kame deposits of the confined aquifer. Withdrawals from small-capacity domestic wells were not included in the model simulation because most of the pumped water is returned to the aquifer through nearby onsite septic systems. Water withdrawn by the production wells is used in sewered areas and then treated at a wastewater facility before discharging to the Chenango River (Town and Village of Greene, New York, 2020).

Horizontal hydraulic conductivity values were assigned on the basis of lithology (table 4). For the valley-fill deposits, final values ranged from 2 ft/d for lacustrine deposits to 200 ft/d for the outwash deposits. Till deposits ranged from 0.5 to 4 ft/d. The ratio of horizontal to vertical hydraulic conductivity was 10 for all surficial deposits and 1 for bedrock.

The altitude and configuration of the simulated water table for long-term, steady-state conditions are shown in figure 17 along with the water level residuals from the limited calibration dataset. The simulated water-table contours and flow directions are consistent with the conceptual model of groundwater flow in the study area. Groundwater flows from topographically high areas and discharges to streams. In the uplands, the water table generally parallels the land surface, and simulated groundwater divides generally coincide with watershed divides. The hydraulic gradient is steepest beneath the upland till and valley-fill deposits near the contact and flattens in the more transmissive valley-fill areas. Simulated groundwater flow in the valley-fill deposits is nearly perpendicular to the main valley stream, Chenango River, and nearly parallel to the small tributary streams that drain the uplands and become minimally gaining or losing as they flow over more transmissive valley-fill deposits. The average residual was −0.3 ft and it ranged from a minimum of −5.2 ft to a maximum of 9.1 ft.

The simulated long-term, steady-state groundwater budget for the modeled area indicated that precipitation recharge accounts for most of the total inflow (28.3 ft³/s or 83.7 percent). Streamflow loss accounts for the remaining inflow (5.5 ft³/s or 16.3 percent). Of the total inflow, 33.4 ft³/s (98.8 percent) of the groundwater discharges to streams and the remainder, 0.4 ft³/s (1.2 percent), is withdrawn by the production wells.

Areas Contributing Recharge

Areas contributing recharge and groundwater traveltimes to the Village of Greene production wells were determined on the basis of the steady-state model for simulated pumping conditions and tracking of pathlines with the MODPATH particle-tracking program. The area contributing recharge to the production wells for the 2010–15 average pumping rates (total of 169 gal/min; table 3) covers 0.35 mi² and includes recharge originating on both sides of the Chenango River (fig. 18). The contributing area extends southeastward from the wells to the simulated groundwater divide in the uplands. The contributing area also includes isolated areas remote from the well from recharge originating on the opposite side of the river from the wells in the till uplands and near the valley-upland contact. Land cover in the contributing area to the wells is mostly agriculture and forested with some urban land uses near the valley-upland contact on both sides of the river (National Land Cover Database, 2016). On the east side of the river this urban land use is considered low density whereas on the west side of the river, this urban land use is low to high density.

Simulated traveltime estimates from recharge locations to the production wells, based on a porosity of 0.35 for the surficial deposits and 0.05 for bedrock, are shown in figure 19. Traveltimes ranged from 3 years to more than 300 years; 9 percent of the traveltimes were 10 years or less. The median traveltime to the wells screened in this relatively deep, confined aquifer was 17.9 years. In addition to the distance from recharge location to the wells, other factors that might affect traveltimes include pumping rates, recharge rates, hydraulic conductivity, and porosity. For example, east of and near the Village of Greene well 1 and well 2, there is a traveltime interval of 20 to 50 years. Groundwater-flow paths from recharge in this interval travel to the wells partly through the low hydraulic conductivity lacustrine deposits. On the opposite side of the river from the wells and near the valley-upland contact, traveltime intervals of 3 to 10 years are adjacent to an interval of greater than 50 years that includes traveltimes of 150 to 170 years. Particle tracks show that recharge originating in the lower traveltime interval mostly travels through the confined kame deposits, whereas recharge originating in the higher traveltime interval travels through the lacustrine and the buried kame deposits.

Hydrogeology

The Cincinnatus study area is characterized by till-covered bedrock uplands that border a north–south trending valley drained by the Otselic River. The altitude of the surface of the shale and siltstone bedrock underlying the surficial materials is shown in figure 21; bedrock-surface altitudes range from about 890 ft in the valley to 1,970 ft in the uplands.

There was minimal information available in the Cincinnatus study area to define the stratigraphy. Thus, the lateral extent and thickness of the lithologic units are uncertain in many areas. Surficial geology (U.S. Department of Agriculture, 2014), lithologic logs, configuration of the land-surface topography, and a regional understanding of how glacial deposits were laid down by the northward retreating ice were used to define the extent and thickness of glacial and postglacial deposits.

The valley consists of unconfined and confined sand and gravel aquifers separated by a laterally extensive fine-grained lacustrine unit; a generalized west-east geologic section through the Otselic River valley between the production wells illustrates these units (fig. 22). The unconfined aquifer consists of mostly outwash deposits, but it also includes kame deposits in several places along the valley walls. In addition, alluvial fans are present where major upland tributaries empty into the valley near the production wells (fig. 20). This unconfined aquifer is generally 30 ft thick in center areas parallel to the valley axis. Lacustrine deposits consisting of very fine sand, silt, and clay can be as much as 40 ft thick but generally are 20 ft thick in center areas of the valley. The confined aquifer, consisting of buried kame deposits, is generally thickest (20 to 40 ft) in center areas of the valley where the production wells are screened. The unconfined kame deposits and postglacial deposits along the valley walls are connected to these buried ice-contact deposits in several areas. The digital hydrogeologic maps developed as part of the present study that serve as the basis for the framework of the groundwater-flow model are available from Woda and others (2022). Data sources and the methods used to develop the hydrogeologic maps are described in Finkelstein and others (2022).

Conceptually, the hydraulic connection between the confined aquifer and the overlying unconfined aquifer, and areas where the buried kame deposits connect to deposits that extend to the surface along the valley walls, are important factors for the location of the contributing areas to the production wells. The intervening fine-grained lacustrine deposits may impede vertical flow from sources of recharge to the unconfined aquifer. For the Town of Cincinnatus wells (well 1 and well 4), recharge to the unconfined kame deposits northeastward and kame deposits and alluvial fan deposits southwestward of the production wells may be an important source of water to the wells. For the Town of Cincinnatus well, well 3, recharge to the confined aquifer from kame and alluvial fan deposits along the valley wall southeastward of the well may be an important source of water. Finally, groundwater inflow from adjacent upland till and shallow bedrock may contribute to the water withdrawn by the production wells.

Development of Groundwater-Flow Model

The groundwater-flow model for the Cincinnatus study area was designed to simulate long-term, steady-state groundwater levels, flow paths, and traveltimes. A basic, manual calibration was done to six groundwater-level altitudes measured in the confined aquifer from 1945 to 2004 (USGS, 2019) at various hydrologic conditions and to a generally shallow water table that approximates the land-surface configuration in the till uplands. The altitudes of the groundwater-level measuring points were based on the altitude of the land surface determined from lidar imagery (Axis GeoSpatial, LLC, 2015; Cortland County, New York, 2017). Because of the limited number and types of observations, model calibration was considered a rough estimate of hydraulic properties. Model characteristics, including hydraulic properties and recharge rates, are summarized in table 4.

The groundwater-flow model subdivided the study area of 62.8 mi² into 472 rows and 560 columns and included a total of 560,620 active cells. Each cell was 125 by 125 ft. In most areas the boundary of the model coincided with watershed divides in the uplands. The large extent of the model minimized the effects of these boundaries on model predictions.

Vertical discretization of the five-layered model was based on lithology (fig. 22). The unconsolidated sediments were subdivided vertically into four model layers. In the valley-fill deposits the top two layers (layers 1 and 2) represented the unconfined aquifer, which consisted mostly of outwash and kame deposits. Layer 3 in the valley represented the semiconfining unit of mostly lacustrine and some kame deposits. Layer 4 in the valleys represented the confined aquifer, which consists of kame deposits. The three production wells were simulated in layer 4. In the uplands, till deposits were represented in all four layers; shallow bedrock areas less than 13 ft from the land surface were simulated as till deposits. Layer 5, the bottom layer, represented only bedrock; the bottom altitude of layer 5 is 100 ft below the bedrock surface.

The model of the Cincinnatus study area specified the same types of boundary conditions that the Cortland and Greene study areas did to represent sources of recharge and areas of discharge. Recharge rates were applied based on the results of a soil-water-balance model (Yager and others, 2018). Mean annual recharge rates calculated by the soil-water-balance model for 2000–13 ranged spatially from 8.9 to 30.2 in/yr in the Cincinnatus study area (fig. 23). In general, lower recharge rates were computed for upland till deposits with low infiltration capacity and higher rates were computed for valley-fill deposits with high infiltration capacity. Mean annual recharge rates averaged spatially over the modeled area was 12.7 in/yr.

Interactions between surface water and groundwater were simulated as a head-dependent flux boundary in layer 1 (fig. 24) by using the streamflow-routing package (Niswonger and Prudic, 2005). Stream altitudes were determined from lidar imagery. Water depth and bed thickness were set equal to 1 ft so that the top and bottom bed altitude could be determined from surface-water altitudes. Simulated streams were assigned widths ranging from 3 to 10 ft for small tributaries, as much as 20 ft for tributaries to Otselic River where they flow across the valley, and 70 ft for the Otselic River. A vertical hydraulic conductivity of 3 ft/d was used to represent bed sediments. Base flow entering the model from the Otselic River area not directly simulated was specified at the first boundary stream cell. Streamflows were specified based on base-flow statistics at the Otselic River at Cincinnatus, New York, streamgage (station number 01510000) (USGS, 2020c). The model contained 5,618 stream cells.

The average annual withdrawal rate for 2010–15 from the three Town of Cincinnatus wells (table 3), which totaled only 34 gal/min, was assumed to represent long-term average rates. The three production wells are screened in kame deposits of the confined aquifer. Withdrawals from small-capacity domestic wells were not included in the model simulations because most of the pumped water is returned to the aquifer through nearby onsite septic systems.

Hydraulic conductivity values were assigned on the basis of lithology. Final calibrated hydraulic conductivity values were similar to the Greene study area (table 4). For the valley-fill deposits, final values ranged from 2 ft/d for lacustrine deposits to 200 ft/d for the outwash deposits. Till deposits were simulated as 5 ft/d. The ratio of horizontal to vertical hydraulic conductivity was 10 for all surficial deposits and 1 for bedrock.

The altitude and configuration of the simulated water table for long-term, steady-state conditions are shown in figure 24 along with the groundwater-level residuals from the limited calibration dataset. The simulated water-table contours and flow directions are consistent with the conceptual model of groundwater flow in the study area. Groundwater flows from topographically high areas and discharges to streams. In the uplands, the water table approximately parallels the land surface, and simulated groundwater divides generally coincide with watershed divides. The water-table gradient is steepest in the upland till and valley-fill deposits near the contact and flattens in the more transmissive valley-fill areas. Simulated groundwater flow in the valley-fill deposits is nearly perpendicular to the main valley stream, Otselic River, and nearly parallel to the small tributary streams that drain the uplands and become minimally gaining or losing as they flow over more transmissive valley-fill deposits. The average residual was −0.1 ft and it ranged from −5.1 to 6.5 ft.

The simulated long-term steady-state groundwater budget for the study area indicated that precipitation recharge provides 58.9 ft³/s or 82.1 percent of the total inflow. Streamflow loss accounts for the remaining inflow, 12.8 ft³/s or 17.9 percent of the total inflow. Of the total inflow, 71.6 ft³/s or 99.9 percent of the groundwater discharges to streams and the remainder, 0.07 ft³/s or 0.1 percent, is withdrawn by the production wells.

Areas Contributing Recharge

Areas contributing recharge and groundwater traveltimes to the Town of Cincinnatus production wells were determined on the basis of the steady-state model, for simulated pumping conditions, and for tracking of pathlines with the MODPATH particle-tracking program. The areas contributing recharge to the production wells for the 2010–15 average pumping rates (total of 34 gal/min; table 3) cover 0.06 mi² (fig. 25). Although the pumping rates for the three production wells are small, the groundwater flux through the confined aquifer is also relatively small in comparison to the unconfined aquifer, thus the wells do not simulate as weak sinks. At these low pumping rates, the areas contributing recharge to each well field are only on the same side of the Otselic River as the wells. Contributing areas to the wells are from recharge originating in the upland till and near the valley-upland contact. The contributing area to the well field on the west side of the river includes small isolated areas the size of a model cell. Analysis of particle tracks indicate that recharge originating in the uplands travels through the buried ice-contact kame deposits to the wells or through the unconfined aquifer, then vertically through the relatively thin lacustrine deposits (3 to 20 ft thick) before travelling to the wells through the confined aquifer. Land cover in the contributing areas to the wells is agriculture and forested (National Land Cover Database, 2016).

Simulated traveltime estimates from recharge locations to the production wells, based on a porosity of 0.35 for the surficial deposits and 0.05 for bedrock, are shown in figure 26. Estimated traveltimes ranged from about 3 years to 180 years; 44 percent of the traveltimes were 10 years or less. The median traveltime to the wells was 11.7 years.

Sensitivity Analysis

The semiconfining unit of fine-grained lacustrine deposits impedes leakage between surface sources of recharge and the confined aquifer. However, there was minimal lithologic data to determine the lateral extent and thickness of the lacustrine deposits and there were no measurements of its vertical hydraulic conductivity, which controls the connection between the unconfined and confined aquifers. A sensitivity analysis of the effects of the degree of hydraulic connection between the unconfined and confined aquifers on the simulated areas contributing recharge to the wells was done by (1) decreasing the vertical hydraulic conductivity of the lacustrine deposits and by (2) a less laterally extensive lacustrine unit near the production wells. These two scenarios were simulated with the 2010–15 total average pumping rate of 34 gal/min but also with a higher total rate of 150 gal/min (50 gal/min for each well) to increase vertical hydraulic gradients between the two aquifers. The alternate simulations are intended to provide insights into the importance of the confining lacustrine unit on delineating the contributing areas to the wells in a narrow, valley-fill setting.

The hydraulic connection between the surface sources of recharge and the confined aquifer was decreased by reducing the vertical hydraulic conductivity of the lacustrine deposits from 0.2 to 0.02 ft/d (ratio of horizontal to vertical hydraulic conductivity from 10:1 to 100:1, respectively), which is a plausible alternate value for these fine-grained, layered deposits. At the low pumping rate of 34 gal/min, the reduced hydraulic connection increased the size of the contributing area to the wells slightly (0.056 to 0.064 mi²) but the wells west of the Otselic River intercepted recharge originating in the uplands more to the northwest than west of the wells in comparison to the higher vertical hydraulic conductivity (fig. 27A). At the higher pumping rate of 150 gal/min in comparison to the lower pumping rate, the contributing area expands to capture enough water to balance the increased pumping rate, including in the till uplands on the opposite side of the Otselic River from well 1 and well 4 (fig. 27). The reduced hydraulic connection increased the size of the contributing areas to the wells for the higher withdrawal rate from 0.21 to 0.29 mi² (fig. 27B).

A plausible, less laterally extensive lacustrine unit near the well fields was simulated with both pumping rates. Lacustrine deposits beneath the alluvial fans west of well 1 and well 4 and east of well 3 (figs. 20 and 22) were instead simulated as kame deposits. This less restrictive connection to recharge sources caused the contributing areas to the wells to decrease from 0.06 to 0.04 mi² at the lower pumping rate and from 0.21 to 0.17 mi² at the higher pumping rate (fig. 28). Where the contributing areas to the wells were present, they generally coincided with the contributing area for the more laterally extensive lacustrine unit.