Upper Raritan Aquifer

The upper Raritan aquifer was introduced by Stumm and others (2024) to represent a unit of Cretaceous sediments that are part of the Raritan Formation present above the clay and silt of the Raritan confining unit and below the basal gravel of the Magothy aquifer in parts of western Long Island. The upper Raritan aquifer was initially identified by Stumm and others (2024) through analysis of high resolution (5-ft interval) core descriptions and gamma logs from a dense network of about 50 vertical profile borings drilled into the Raritan confining unit to help define the comingled plumes in southeastern Nassau County (fig. 6). The analysis indicated a substantial deviation of the top surface of the Raritan confining unit as compared to that in previously published maps (Suter and others, 1949; Smolensky and others, 1990). Dense clay beds clearly were absent in the upper part of the Raritan confining unit, indicating a possible facies or depositional change from fine- to coarse-grained lithologies in this part of the Raritan Formation.

The newly mapped upper part of the Raritan Formation is generally characterized by fine to coarse sand with silt and minor clay lenses below the Magothy aquifer’s basal gravel, whereas the Raritan confining unit has been previously described as consisting of dense clay and silt with small sand interbeds (Suter and others, 1949; Smolensky and others, 1990). Further analysis of hundreds of core descriptions and gamma logs in Kings, Queens, and Nassau Counties by Stumm and others (2024) indicated the upper Raritan aquifer extends throughout the southernmost parts of these counties. The aquifer also appears to extend throughout the central part of Nassau County (fig. 6) and may extend eastward into Suffolk County.

Raritan Confining Unit

The unnamed clay member of the Cretaceous Raritan Formation forms the Raritan confining unit (fig. 7), which underlies most of Nassau County (including the entire inset area) and dips to the southeast (Suter and others, 1949). The unit consists of dense clay and silt that are multicolored, including gray, white, red, and tan (Suter and others, 1949; Stumm, 2001; Stumm and others, 2002, 2004). The top of the Raritan confining unit was defined in Stumm and others (2024) and in this study by the presence of a dense clay and silt that is at least 20 ft in thickness based on core samples and gamma logs.

Distribution of Hydraulic Properties

Nuclear magnetic resonance (NMR) logs measure the induced magnetic moment of hydrogen protons contained within the fluid-filled pore space of porous media. Unlike other geophysical logs that respond to the rock matrix and fluid properties and are strongly dependent on mineralogy, NMR logs respond to the presence of hydrogen protons in the formation fluid to determine water content (equivalent to porosity in the saturated zone) and pore-size distribution (derived from bound versus mobile water) (Dlubac and others, 2013). NMR logs can be used to estimate vertically continuous porosity and hydraulic conductivity in the adjacent formation using equations and compilations of regressions determined from published collocated hydraulic tests of unconsolidated aquifers (Walsh and others, 2013; Knight and others, 2016; Kendrick and others, 2021).

NMR logs were collected from five deep groundwater monitoring wells in Nassau County (fig. 8) and were analyzed to determine clay-bound, capillary-bound, and mobile-water content or fraction and to estimate hydraulic conductivity of aquifers and confining units. More details about the collection and analysis of NMR logs from these wells, including the estimation of hydraulic properties using published co-located hydraulic tests, are provided in Stumm and others (2024).

Upper Glacial Aquifer

The upper glacial aquifer south of the glacial moraines consists predominantly of well-sorted outwash and is more hydraulically conductive than sediments north of the moraine where the aquifer consists of variably silty ice-contact deposits (fig. 9). The geometric mean horizontal hydraulic conductivity (Kh) estimated from 121 pumping tests in the glacial outwash was 200 feet per day (ft/d), whereas the geometric mean horizontal hydraulic conductivity estimated from 22 pumping tests north of the moraine was 120 ft/d (Stumm and others, 2024). The total porosity outwash (south of the moraine) ranged from 30 to 40 percent based on laboratory analysis of several hundred samples from southern Long Island reported by Veatch and others (1906).

Magothy Aquifer

The geometric mean Kh of the Magothy aquifer estimated from 35 pumping tests was 120 ft/d (Stumm and others, 2024). The geometric mean Kh estimated from the NMR logs of monitoring wells N 12523, N 14421, DEC-HC05-VPB02 (Pollen well), and N RW8-MW01D3 (fig. 8) was 140, 81, 88, and 110 ft/d, respectively (U.S. Geological Survey, 2022a, 2022b, 2022c). The geometric mean Kh ranged from a low of 0.33 ft/d for finer-grained intervals to a high of 590 ft/d for coarse-grained intervals. The mean total porosities of the Magothy aquifer estimated from the NMR logs of the four wells were 0.33, 0.38, 0.38, and 0.38, respectively (Stumm and others, 2024). The estimated mean mobile-water fractions of the Magothy aquifer in the four wells were 0.26, 0.24, 0.25, and 0.28, respectively. The mean mobile and bound water fractions for coarser-grained intervals were 0.34 and 0.03, whereas the mean mobile and bound water fractions for finer-grained intervals were 0.02 and 0.38, respectively.

Texture Model of the Upper Glacial and Magothy Aquifers

We completed a detailed characterization of lithologic descriptions and the conversion of these qualitative data to more quantitative measures of important hydrogeologic characteristics for the upper glacial and Magothy aquifers. This detailed analysis was centered on the inset area and included the development of a local-scale model of changes in sediment characteristics referred to herein as a “texture model.” This analytical approach generally followed the one used by Walter and Finkelstein (2020) to develop a model of changes in the sediment distribution correlated to changes in hydraulic conductivity for the entire Long Island aquifer system.

Lithologic Data Analysis

We generally followed the analytical approach of Walter and Finkelstein (2020) with 22,523 lithologic descriptions from 256 boreholes in the inset area (fig. 10) to define standard lithologic codes (table 1 in Walter and Finkelstein, 2020) for each vertical interval in each borehole. These coded intervals were assigned values of Kh and Kv from previous investigations (Walter and Finkelstein, 2020) and were assigned a binary variable representing the presence or absence of silt and clay. Each borehole was divided into 10-ft lengths, and mean hydraulic conductivity values were calculated for each interval; values in each 10-ft layer represented a thickness-weighted mean horizontal and geometric mean vertical calculated from all lithologic intervals coincident with the layer. These borehole descriptions were compiled into a geographic information system database containing borehole locations and mean hydraulic conductivity values.

Texture Model Development

Following the methods of Walter and Finkelstein (2020), stacked grids of Kh and Kv were created, each representing 10 ft of aquifer thickness spanning the upper glacial and Magothy aquifers in the inset area (fig. 11). Grids representing the upper glacial aquifer were horizontal, whereas those representing the Magothy aquifer were sloping and were draped onto a surface representing the bottom of the unit (fig. 11). A single grid represented the hydraulic conductivity over the full thickness of the upper Raritan aquifer (fig. 11).

We queried a geographic information system database (Walter and Finkelstein, 2020) to extract those boreholes that extended to or beyond each successive grid. These queried subsets of data points were used to interpolate values of Kh and Kv by use of ordinary, spherical kriging. The vertical stack of interpolated grids constitutes a quasi-three-dimensional texture model of each hydrogeologic unit (fig. 12) that, when combined, form a detailed representation of the principal aquifer system of the inset area (fig. 13).

Distribution of Hydraulic Conductivity

We modeled Kh and Kv, both important controls on groundwater flow, following the methods of Walter and Finkelstein (2020). The spatial and vertical patterns in Kh within the inset area (figs. 11–13) generally are consistent with the depositional history of the regional aquifer system.

Values of hydraulic conductivity in the upper glacial aquifer generally are lower in northern parts of the inset area and are associated with glacial moraines (fig. 12A). Hydraulic conductivity is highest in outwash sediments south of the moraine (fig. 12A). Outwash sediments consist largely of well-sorted sand; moraine sediments generally consist of fine-grained and poorly sorted sediments. In addition, the glacial sediments generally become finer with depth (figs. 13 and 14).

The Magothy aquifer within the inset area is more heterogenous than the upper glacial aquifer, with fewer broad spatial patterns, except in its basal portion. Hydraulic conductivity in the basal portion of the Magothy aquifer (fig. 12C), where sediments likely were deposited in fluvial depositional environments, generally is higher than in overlying portions of the unit (figs. 12B and C). Hydraulic conductivity values generally are lowest in the middle part of the Magothy aquifer (fig. 12B) where fine-grained sands and silts with interbedded clay lenses likely were deposited in overbank lake and wetland environments.

Transition Probability Approach

The transition probability approach was developed to better describe the relation between observable features and model parameters. The transition probability measures the spatial variability (for example, given that a facies j is present at a location x, what is the probability that another [or the same] facies k occurs at location x+h). The approach can consider all cross-correlation information in three dimensions instead of just two dimensions (Carle, 1999).

Simulation Descriptions

T–PROGS was used to generate multiple equally probable models of aquifer heterogeneity all conditioned to the 256 boring logs. The goal was to create multiple material sets that support the stratigraphic framework of the groundwater model.

A three-dimensional grid was set up in GMS with 288 rows, 224 columns, and 26 layers to be coincident with the inset models (table 3, fig. 15). There are 256 boring logs with 22,838 individual intervals described that constitute the framework for the T–PROGS simulation. Each interval was assigned one of four categories (gravel, sand and gravel, sand, or fines; as described previously in the “Sediment Categories” section).

The most important step when running T–PROGS is defining the transition probability data for each material in the three primary directions: vertical, strike, and dip. The vertical transition is developed first based on the borehole data. The data in the strike and dip directions are then derived from the vertical data (Carle, 1999). The distribution of the boring logs within the grid shell, the vertical extent of the boring logs, and the sediment types within each boring log are shown in figure 16. Sand is the most prevalent category with a distribution percentage of 60 percent (table 4, fig. 16). Sand and gravel (fig. 16) is the second most dominant sediment type with a distribution percentage of 24 percent. Sand and sand and gravel combined account for 84 percent of the total sediment types within the boring logs. Gravel and fines make up a distribution percentage of 5 and 11 percent, respectively. These sediment type distribution percentages were similar for each potential geostatistical realization computed by T–PROGS.

Prior to creating a geostatistical realization, the GAMEAS utility was used to compute a set of transition probability curves as a function of lag distance for each material for a given sampling interval (Carle, 1999). Using a lag spacing of 0.3 ft, a maximum lag distance of 450 ft was computed, meaning the greatest distance a material correlates horizontally is 450 ft. In addition, the transition probability for each material with respect to each of the other materials is calculated and displayed graphically (fig. 17).

The probability is shown on each plot on the y-axis from 0 to 1, and the correlation distance on the x-axis (0 to 450; fig. 17). The top right graph indicates that there is very low probability to transition from gravel to fines throughout the 450-ft distance. The upper left plot indicates that there is a high probability of gravel adjacent to gravel at close distances and then the probability rapidly declines with distance. The goal of this portion of the analysis is to closely fit the Markov chain to the measured transition probability curve. The embedded transition probabilities derived from the calculated transition probability plots (fig. 17) are shown in table 5.

The next step in generating a realization is to define the Markov chains in the strike and dip directions. We applied Walther’s Law (Carle, 1999) to estimate the strike and dip Markov chains because borehole data in these directions are not sufficiently dense to calculate them in the same way as vertically. Walther’s Law states that vertical successions of deposited facies represent the lateral succession of environments of deposition. An assumption is made that the proportions are the same in all three directions. The lens length ratio is used to define the transition rate matrix with the diagonal transition rates defined from the lens lengths. These ratios calculate the lens length in the horizontal (strike and dip) directions based on the observed lens length in the vertical direction.

We used a trial-and-error method for determining the lens length ratios in the strike and dip directions to allow for appropriate connectivity between the borings (table 6). Different values were used for each sediment type in the lens length ratio and an example realization was created. This process was repeated until a realization was produced by T–PROGS that compared reasonably with the cross-sections interpreted by other methods, the borings, and knowledge of Long Island hydrogeology (Walter and Finkelstein, 2020). An orthogonal view of an example realization produced from the final lens length ratios used in the strike and dip directions is shown in figure 18.

Once the final lens length ratios were determined, T–PROGS was used to create 500 realizations based on the same input structure in a Monte Carlo type analysis. These 500 realizations were then used in the inset area groundwater-flow model to determine the potential range of, and most likely plume path based on, interpolated preferential flow paths.

Final Realizations

T–PROGS was run 500 times after a realization that compared reasonably to the data, cross-section, and knowledge of Long Island’s hydrogeology was created. These 500 realizations all held true to the boring logs at their locations and were interpolated in the vertical, strike, and dip between boring locations using the already determined transition probabilities and Markov chains. Six realization examples and the corresponding sediment type distributions for each example are shown in figures 19 and 20. Each realization varies slightly; however, the percent distributions remain consistent with the material distributions in the boring logs.

One of the 500 potential realizations shows channels of gravel and sand-gravel that stretch throughout the grid as would be expected in a deltaic deposit highlighting the preferential pathways that the plume at the NWIRP and NGBF sites could potentially follow (fig. 21). The areas between isosurfaces are areas of lower hydraulic conductivity and therefore are not considered preferential flow paths.

The T–PROGS realization is constrained by the location of the boring log information. Beyond the grid (fig. 16), the results are no longer as constrained because there are no borings for T–PROGS to use to guide the realization. The realization beyond the boundary of the existing borings is extrapolated from the boring information using the Markov chains and distribution percentages. Overall, the final compiled realizations match the boring logs closely and represent the known geology of the area.

Transient Recharge

We estimated recharge for each month from January 2005 through December 2019. Recharge components included natural and redirected recharge, return flow, and leakage from water-supply infrastructure. These recharge components were combined to estimate the total monthly recharge for the 2005–19 period.

Natural and Redirected Recharge

Finkelstein and others (2022) used the SWB model (Westenbroek and others, 2010) to estimate natural and redirected groundwater recharge for the Long Island aquifer system. Monthly recharge generally follows a seasonal cycle with most of the recharge occurring during the nongrowing season (fig. 22). Anomalies include Tropical Storm Tammy, which added significant recharge during October 2005; the wet summer of 2006; and a March 2010 Nor’easter (Finkelstein and others, 2022). Tropical Storm Tammy resulted in about a sixfold increase in the monthly mean recharge rate, and the March 2010 Nor’easter resulted in about a fourfold increase in this rate.

In the more urbanized areas of Long Island, impervious surfaces intercept precipitation, resulting in less aquifer recharge as compared to predevelopment (prior to 1900 in Nassau County) conditions. The largest recharge deficits attributed to interception from impervious surfaces occur in the urbanized areas of New York City and southern Nassau County; however, much of the recharge potentially lost to impervious surfaces in Nassau County can recharge the aquifer through sumps, dry wells, and a network of more than 6,000 recharge basins (Walter and others, 2020).

Monthly estimates of redirected recharge from impervious surfaces were computed by subtracting SWB computations of the steady-state predevelopment condition from the 2005 to 2019 condition, resulting in a total redirected recharge of about 3 percent of natural recharge under present conditions (Finkelstein and others, 2022). The temporal pattern of redirected recharge (fig. 23) is like that of natural recharge (fig. 22) because of the correlation of predevelopment and 2005–19 conditions recharge with seasonal patterns of precipitation. In addition to Tropical Storm Tammy and the March 2010 Nor’easter, the Halloween Nor’easter of October 2011 also generated notable redirected recharge volume.

Wastewater Return Flow

We estimated wastewater return flow using population density derived from the 2010 census in conjunction with annual water use (Walter and others, 2020). About 85 percent of public-supply water use on Long Island is returned annually as wastewater return flow to the aquifer system by way of onsite domestic septic systems and constitutes a spatially and temporally distributed supplemental recharge. In areas served by public sewering, we assumed that 1 percent of wastewater recharged the aquifer from leaky sanitary sewer lines (Walter and others, 2020). Wastewater return flow was about 5 percent of total recharge to the water table during the 2005–19 period. The monthly wastewater return flow component of total recharge increased steadily from 130 to 140 cubic feet per second (ft³/s) during the 2005–19 period (fig. 24).

Water-Supply Infrastructure Leakage

An additional source of artificial recharge to the aquifer system is leakage from public-supply infrastructure. Most of Long Island is served by public supply; private water supplies do exist but are limited to small areas of eastern Suffolk County and the northern shore of Nassau County. The distribution of public-supply lines was estimated from the extent of roads within individual model cells (Walter and others, 2020). The potential for recharge from leaky infrastructure (as indicated by road length) generally reflects patterns of population and development and is largest in New York City and parts of Nassau County. In Kings and Queens Counties, about 700 Mgal/d of water is imported from a reservoir system in upstate New York, and assuming a 10-percent loss of water from the water distribution system, may result in about 70 Mgal/d of additional recharge on western Long Island (Misut and Monti, 1999).

Water-supply infrastructure leakage accounted for about 2 percent of natural recharge during the 2005 to 2019 period. The large seasonal range in this leakage term (fig. 25) is correlated with the seasonal variation in water use (fig. 26).

Groundwater Withdrawals

An average of about 424 Mgal/d of groundwater was withdrawn annually from 2005 to 2015 from the Long Island aquifer system for multiple uses, including public supply, agriculture, and industry. The public supply of drinking water accounted for nearly all (95 percent) of the total annual groundwater withdrawal on Long Island. Agricultural withdrawals for crop irrigation generally were limited to eastern Suffolk County, with limited reporting (Walter and others, 2020). However, potential agricultural pumping was represented using the SWB model (Finkelstein and others, 2022) by estimating water demand by crop type and available precipitation, resulting in an estimated demand of about 2 Mgal/d. Seasonal cycles in pumping (fig. 26) relate to increased water demand during the summer caused by increased inhabitation, irrigation, and lawn irrigation demand.

During the 2005 to 2019 period, a total of about 8 Mgal/d of groundwater on average was withdrawn annually from the Long Island aquifer system for contaminant remediation (Corson-Dosch and Fienen, 2023). This remedial pumping generally does not vary seasonally and has increased with time; all water generally is returned to the aquifer system after treatment either by reinjection wells or discharge to recharge basins. The vast majority (7.9 Mgal/d) of the total was withdrawn as part of Navy Grumman groundwater plume remediation in southern Nassau County.

Hydraulic Properties

The transient LIRM model requires specification of unconfined and confined storage. We assigned a uniform value of 0.25 to represent specific yield (Sy; unconfined storage) and a uniform value of 0.00001 to represent specific storage (Ss; confined storage), which are reported for Long Island sediments (Smolensky and others, 1990). All other hydraulic properties are unchanged from the values reported in the steady-state model of the LIRM (Walter and others, 2020).

Inset Model Domain

The inset model domain encompasses an area of about of 36 mi² in southeastern Nassau County surrounding the plume. Surface elevations are highest in the northeastern part of the inset domain near the hamlet of Plainview and Bethpage State Park and generally decrease from north to south towards South Oyster Bay and the Atlantic Ocean. Streams drain the inset domain and typically flow parallel to the north-to-south slope. Larger streams in the inset domain include Bellmore Creek and Massapequa Creek, smaller streams include Seamans Creek, Seaford Creek, and Carman Creek. Although stream channels (expressed as a drainage pattern in the surface elevation model) extend into the northern part of the inset domain, stream reaches with perennial flow are only present in the southern part of the inset domain (fig. 32). There are wetlands adjacent to anastomosing sections of Bellmore Creek and Massapequa Creek near the coast in the southern part of the inset domain, which have shallow depth to groundwater.

Hydrogeologic units in the inset model domain include unconsolidated glacial and coastal-plain deposits and are described in detail in the “Hydrogeologic Framework” section. Permeable sand and gravel glacial deposits make up the surficial (upper glacial) aquifer. The glacial deposits are underlain by a thin, discontinuous marine clay layer in the southern-most portion of the domain (Gardiners Clay) that may act as a confining unit. Both units are underlain by a thick, regionally extensive coastal-plain deposit (Magothy aquifer) consisting of mostly sand and silt with a basal layer of gravel. This basal gravel is underlain by another unit of sand and silt (upper Raritan aquifer), which in turn is underlain by a fine-grained, clay-rich deposit (Raritan confining unit) that acts as a confining unit to the Lloyd aquifer below. Impermeable crystalline bedrock underlies the Lloyd aquifer and is considered to be the base of the aquifer system.

The area encompassed by the inset model domain is highly developed with extensive residential, commercial, and industrial land use. Aerial recharge is limited by impervious surfaces that generate excess stormwater runoff. Many large, human-constructed basins have been installed in developed areas throughout the model domain for stormwater retention and infiltration (Aronson and Seaborn, 1974). Numerous sole-source water-supply wells are active in the inset model domain, managed by distinct water utilities. An existing network of remedial pumping and reinjection wells are also active as part of a groundwater extraction, treatment, and return remedial system (New York State Department of Environmental Conservation, 2024a, 2024b).

Inset Model Discretization

The inset model extent adopted the bounds of the southeastern Nassau County inset area shown in figure 2. The layering adopted the same lithologic definitions as in the transient LIRM with the exception that layer 24 of the transient LIRM was subdivided into two layers in the inset model. The purpose of the two layers was to divide the Raritan confining unit into a clayey unit (inset model 25, the “lower Raritan”) and a sandy unit (inset model 24, the “upper Raritan”). Inset model layer 26 represents the Lloyd aquifer, which is consistent with layer 25 in the transient LIRM (table 9). All layers shallower than layer 24 are lithologically consistent with the transient LIRM (Walter and others, 2020).

We refined the horizontal discretization from the square 500 ft by 500 ft cells in the transient LIRM to square 125 ft by 125 ft cells. This is a 16-fold refinement of the transient LIRM geometry and is intended to facilitate higher resolution representation of the relation among remedial wells, the plume, surface waters, and boundary conditions.

The transient temporal discretization adopted the same 180 monthly stress periods as used in the transient LIRM, representing each month from January 1, 2005, through December 31, 2019. We implemented an initial steady-state stress period at the beginning of the simulation, which was constructed using the average stresses and boundary conditions of the entire 2005–19 period. More detailed information on the transient stress periods that were represented in the transient LIRM is presented in the “Conversion from Steady State to Transient Simulations” section.

Boundary Conditions

Boundary conditions that do not change from the transient LIRM model are mapped by MODFLOW-setup to the overlapping model cells in the inset model. These boundary conditions are presented on a map of the inset model area in figure 32.

Recharge

Recharge was inherited from the transient LIRM at the 500 ft by 500 ft resolution represented by the transient LIRM. As a result, each 16-cell block of inset model cells contained within a single transient LIRM model cell used the same recharge value. Average annual recharge rates for the simulation period (2005–19) are shown in figure 33.

Focused Recharge

Several recharge basins (fig. 34) were simulated in the inset model to represent locations where treated water extracted in the remedial system was returned to the subsurface. We represented these in the inset model using the MODFLOW 6 Well (WEL) package with a return-flow (or reinjection) well implemented in each cell in the shallowest layer that overlaps with the outline of a recharge basin. We set hydraulic conductivity in these basins to an initial value that was high enough to accommodate the known return flows of treated water without simulating excessive mounding. This hydraulic conductivity parameter was allowed to vary in the history-matching process.

Hydraulic Properties

The hydraulic properties applied to the inset model include hydraulic conductivity, storage, and porosity. The initial values of these hydraulic properties were adjusted using history matching, as described in the “Parameter Estimation by Ensemble History Matching” section. Hydraulic conductivity was represented using two methods: a single realization generated using geostatistical interpolation referred to as the “texture model,” and an ensemble of 500 realizations generated using the transition probability geostatistical method (T–PROGS; Carle, 1999). Both conceptualizations are described in detail in the “Hydrogeologic Framework” section.

We specified specific yield as a constant value of 0.25 throughout the model and specific storage also was specified as a constant value of 1x10⁻⁵. Both uniform values were inherited from the transient LIRM. A single value for porosity of 0.35 was selected for the analysis in this report. Travel times were not reported in the particle-tracking simulations, and those model outputs are the only outputs that would be affected by changes in porosity.

Observation Data

Observations of groundwater heads and base flows were used for history matching for the 2005–19 period. The observation dataset consisted of three categories: (1) base flow in streams, (2) hydraulic heads in wells, and (3) temporal hydraulic head differences. Base-flow observations were derived from streamflow records at the streamgage on Massapequa Creek (01309500) and two streamgages on Bellmore Creek (01309950 and 01309990) shown in figure 27 and table 8 (U.S. Geological Survey, 2023).

Hydraulic head observations were collected from multiple sources, as outlined in Corson-Dosch and Fienen (2023). These observations were supplied at various temporal intervals ranging from occasional instantaneous values to some collected every 15 minutes for many months. To map these values to the monthly temporal discretization of the inset model, we selected the final (in time) measurement of groundwater head in each month as representative of that month in the model. Temporal differences were finally derived by subtracting sequential groundwater head observations in each well from each other where more than one monthly measurement was available. These temporal differences can be valuable in amplifying the effects of storage on the timing of aquifer dynamics (Doherty and Hunt, 2010; Anderson and others, 2015).

Observations are assigned weights to indicate how closely model outputs are expected to correspond to the observed values. We initially assigned these weights based on assumed standard deviations of 10-percent coefficient of variation for the base flow observations, 1.5 ft for groundwater head observations, and 0.5 ft for temporal groundwater head difference observations. These initial weights were calculated as the reciprocals of the assumed standard deviations and resulted in an imbalanced objective function. Information about observations and weights are summarized by group in table 10. Subjectively, weights were adjusted to rebalance the objective function for the texture model and T–PROGS models (fig. 35) so that base flow, hydraulic heads, and temporal hydraulic head differences make up 40, 40, and 20 percent of the initial objective function, respectively.

The ensemble approach of iES also allows for realizations of variability around the reported observation values, referred to below as “observation noise.” By default, the reciprocal of the weights is used as the standard deviation which, combined with the observation value as the mean, can be used by the PEST++ software to formulate a normal distribution about which a sample including observation noise for each observation is drawn for each ensemble realization.

Rebalancing the objective function and thus changing weights can result in standard deviation specifications that are inappropriate for the assumed variability around the observation values. Further, the default behavior in PEST++ does not include autocorrelation of the observation noise. To account for this, we calculated the covariance of the time series for all transient observations assuming an exponential variogram with unit variance and a correlation length of 140 days, and then scaled by the variance (the standard deviation values ascribed to each observation in table 10, squared). Examples of the observation noise realizations calculated in this way are shown in figure 36. In figure 36A, which includes temporal autocorrelation, the realizations are smoother, and more consistent with error on the time series that should vary smoothly, whereas in figure 36B the uncorrelated observation noise can suggest errors jumping above and below the reported values.

Estimating Model Parameters

We adopted a multiscale approach to parameterizing inputs of properties in the inset model based on the methods described in White and others (2020) and Fienen and others (2022b). In this approach, most parameter values were defined as multipliers that are multiplied together against the initial reference property values. In the specific case of hydraulic conductivity, this allowed us to apply parameter adjustment at two scales: homogeneous zones and pilot points (Doherty, 2003). Different scales of parameters respond differently to information contained in observations so there is an advantage to designing the history-matching setup to accommodate this. The summary of initial parameter values, bounds, and numbers of parameters for the texture model and T–PROGS inset models are shown in tables 11 and 12, respectively.

We developed separate inset models with the hydraulic conductivity distribution from the texture model and T–PROGS. The texture model was created using more and thinner layers than the layering in the inset model. As a result, the texture model was mapped to the inset model layers using a thickness-weighted arithmetic mean for Kh and a thickness-weighted geometric mean for Kv. The T–PROGS model was generated with the same layer configuration as the inset model, so the mapping was one-to-one in this case except for the Lloyd and Raritan units, which were omitted from T–PROGS calculations. We parameterized these two layers from a homogeneous starting point.

In both cases, the hydraulic conductivity values were adjusted such that their mean value was the same as the mean of the corresponding layer in the transient LIRM. This adjustment was done for all layers except for layer 4 (the Gardiners Clay, which was represented as a single value) and layer 24 (the representing the upper Raritan, which was not represented in the transient LIRM). As a result, each of the 500 initial realizations for the T–PROGS model has the same bulk mean hydraulic conductivity per layer, and because the relative proportions of each of the four zones in each layer are consistent among the realizations, the initial hydraulic conductivity values per zone per layer were nearly identical at the start.

For the texture model, we assigned a zonal multiplier to each layer and a network of pilot points in each layer for Kh and Kv. A subsequent analysis confirmed that Kv was less than Kh in all cells. For the T–PROGS model, each layer has four discrete values, and as a result, unique values were assigned to each of the four zones in each layer (for layers 1–23 except for layer 4, which was simulated as homogeneous) in addition to the same network of pilot points and layer-wise constant multipliers as for the texture model.

Each of the 500 realizations used for the texture model represented multipliers that expressed variability around an initial hydraulic conductivity distribution. An additional parameter was required in the T–PROGS inset model to indicate which of the initial 500 realizations of hydraulic conductivity was the starting point for multipliers to express variability around them. This parameter was set at an integer value and not allowed to change through the history-matching process.

The remaining parameterization was the same for the texture and T–PROGS models. We implemented a single multiplier, each, for 181 transient parameters for the boundary CHD, WEL, and recharge (RCH; Langevin and others, 2017) MODFLOW packages for each stress period. These multiplier parameters were restricted to the range of 1.0 to 1.2 for CHD, 0.9 to 1.1 for WEL, and 0.6 to 1.1 for RCH, respectively (table 11). Direct values (not multipliers) were assigned for hydraulic conductivity in the three recharge basins where the return of treated water to the aquifer was simulated, resulting in three Kh and three Kv parameters used in the history-matching process.

The same pilot-point network was used for all Ss, Sy, Kh, and Kv parameters. This pilot-point network included one pilot point assigned for every 15th row and column, resulting in a total of 6,650 parameters for each of the four properties spaced every 1,875 ft apart. Finally, vertical multipliers were assigned by model layer for Ss and Sy. The total number of adjustable parameters was 27,317 for the texture model and 27,501 for the T–PROGS model. The additional parameters in the T–PROGS formulation account for individual homogeneous Kh values we assigned to the four discrete zones per layer, which we allowed to vary independently from the pilot-point multipliers.

History Matching with the Iterative Ensemble Smoother

We selected the iES (White, 2018) as implemented in PEST++ (White and others, 2021) for history matching. This technique builds on the prior Monte Carlo analysis by evaluating an ensemble of parameter values, yielding a residual matrix, 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 repeated until a posterior conditioned ensemble is obtained. This process results in an ensemble of objective function values from iteration to iteration.

The “Base” line on figure 37 represents a base realization ensemble member (discussed further in the “Base Realization” section) that is considered representative of the entire ensemble when a single realization is desired. This base ensemble member is calculated by iES by applying the approximate upgrade information to the initial parameter values and using observations unperturbed by noise to approximate a single minimum-error-variance solution. A base realization can be estimated for the texture model inset because it started with a single initial realization of hydraulic conductivity. However, a base realization of hydraulic conductivity cannot be estimated for the T–PROGS model inset because a single initial realization does not exist.

History-Matching Ensemble Results

The results of the first iES iteration were selected as optimal based on a subjective judgement that assumes the objective function had achieved an acceptable level of fit while still retaining enough variability in the ensemble to characterize the uncertainty in the system. Individual model runs may fail to converge during the history-matching process because of the randomized generation of parameter sets. When a model run fails to converge, iES drops that run from the ensemble. Most texture model inset realizations converged through the first iteration of the history-matching process, but many T–PROGS model inset realizations (209 out of 500) did not. These nonconverging runs likely were due to the large variability in the T–PROGS model inset hydraulic conductivity (Kh) zonations that produced unrealistic model properties when combined with randomized parameter sets.

Ensemble members with unacceptably high phi values were removed through the process of rejection sampling as a final step towards obtaining a posterior parameter ensemble. We selected a separate phi cutoff for each of the texture and T–PROGS model insets that removed anomalously high phi results and yielded an approximately normal distribution (fig. 38). This reduced the number of realizations retained in the posterior ensembles from the initial size of 500 realizations to 459 for the texture model inset and 268 for T–PROGS model inset.

Anomalously high phi values are common in the prior Monte Carlo analysis, and because high phi realizations are either eliminated or move toward the mean of the ensemble, more normally distributed phi values become apparent as the algorithm progresses through iterations. The T–PROGS model inset realizations are less constrained than the texture model inset, resulting in more realizations with phi values being dropped and fewer realizations remaining for the posterior ensemble than for the texture model inset.

The history matching discussion includes the full range of results from the entire posterior parameter ensemble and the results from the single representative base realization. The entire posterior parameter ensemble is used in the discussion of head and base-flow observations and for boundary condition parameters, where the visualization of the full range is possible. The base realization is preferrable for visualizing properties and water budget results. Results from both models are discussed in this section; however, only texture model inset results are shown for a base realization because the T–PROGS model inset had no base realization and, therefore, no single representative realization.

Heads and Flows Residuals

Heads and flows residuals are calculated from the entire ensemble rather than focusing on a single realization; therefore, these residuals represent a range of outcomes consistent with prior knowledge of the parameters and the information from the observation data. The correspondence between model outputs and observations are grouped separately for the texture model and T–PROGS model insets. The base realization and the spread from all the realizations are shown on each plot as a dot and whiskers, respectively (figs. 39 and 40). The fit for the texture model inset base realization is shown in figure 39, whereas the fit for the T–PROGS model inset realizations closest to the median of the ensemble, based on weighted residuals, are shown in figure 40. There is more variance and bias in modeled outcomes with the T–PROGS model inset than with the texture model inset (figs. 39 and 40); however, these results are shown to have little bias and low variance relative to the uncertainty of the observations.

The correspondence between modeled flow in SFR cells at streamgage locations and observed base flow indicate that the ensemble of model outputs generally overlaps with the ensemble of observations (fig. 41). There is less scatter for all cases in the modeled results than in the observed data (for example, the two streamgages for the texture and T–PROGS model insets); overall, there is little bias, and the fit is considered reasonable.

Transient hydraulic head and temporal hydraulic head-difference results (figs. 42 and 43) for selected locations shown in figures 44 and 45 were generated for the texture and T–PROGS model insets. Although there is some bias relative to the measured observation values, the modeled results are generally completely enveloped by the range of observation noise as shown by the modeled ensemble limits region being completely inside the observation ensemble limits region (fig. 42C, D, G, H).

A similar pattern is observed in the temporal hydraulic head-difference targets as seen with the transient hydraulic head targets; however, the temporal hydraulic head-difference results are less biased than those of the transient hydraulic head results (figs. 42A, B, E, F and 43A, B, E, F). Temporal hydraulic head-difference observations provide insight into the dynamics of the aquifer system over time, so even when there is some overall bias in transient hydraulic head observations, the dynamics can still be represented in meaningful ways to improve model prediction capabilities.

The average residuals (modeled minus observed) of hydraulic heads and flux observations were computed over time and over the entire ensemble to display these results spatially. Maps of the spatial residuals grouped by hydrogeologic unit are shown on figures 44 and 45 for the texture and T–PROGS model insets, respectively. The residuals generally are low (close to zero) and randomly distributed for the texture and T–PROGS model insets, indicating low bias, which is a goal of the history-matching process.

Hydraulic Properties

The post-history matching Kh and Kv from the base realization are shown in figures 46 and 47 for the texture model inset. Texture model inset Kh values agree favorably with literature values (Walter and Finkelstein, 2020) and values used in this area in the transient LIRM. The highest Kh values were typically in the upper glacial aquifer in represented in model layer 1, with a maximum value of 375 ft/d. The lowest values were in the layers representing fine-grained confining units including the Gardner’s Clay (layer 4, average Kh 0.009 ft/d) and the lower Raritan (layer 25, average Kh 0.01 ft/d). Kv is typically one to two orders of magnitude lower than Kh in the texture model inset. The average vertical anisotropy (the ratio of average Kh to average Kv) across all model layers was about 62 to 1 and ranged from 120 to 1 in model layer 1 to about 9 to 1 in model layer 9. These ranges are considered appropriate for the layered glacial- and marine-deposited sediments that make up the aquifers in the inset model domain. Kv was adjusted below recharge basins to reduce water-table flooding in these locations and to better represent measured infiltration rates by Seaburn and Aronson (1974). Kv was adjusted to 200 feet per day in layer 1 cells that contained recharge basins (in texture model and T–PROGS model insets). These adjusted cells exceed the upper bound of the color bar on this figure, which was intentionally limited to 40 feet per day to better illustrate Kv variably within and between other model layers.

The Ss and Sy for the optimal base realization of the texture model are shown in figures 48 and 49, respectively. Literature values and the ranges used in the transient LIRM were used to set history-matching bounds for Ss and Sy because of the limited data available to constrain these storage properties (Walter and others, 2024). In MODFLOW 6, both storage properties (Ss and Sy) must be defined for each model cell. The Sy value then is used to calculate water exchange with storage when the water table is below the top of the cell; otherwise, Ss is used for this calculation. The water table is almost exclusively present within the top three layers in this model, so Sy values are summarized and visualized only in layers 1–3, whereas Ss values are summarized and visualized in layers 4–26. The average Ss was 9.7e-5 ft⁻¹ across layers 4–26 and reached the lower and upper bounds of 1.0e-6 ft⁻¹ and 1.0e-3 ft⁻¹ in at least some places (fig. 49). Generally, Ss values are consistent with changes in lithology; however, close alignment of estimated Ss with lithology is not always achieved during the history-matching process. The average Sy was 0.276 in layers 1–3 and was highest in the top two model layers representing the upper portion of the upper glacial aquifer (fig. 48). Some Sy values reached the lower and upper bounds of 0.238 to 0.290, respectively; however, the Ss and Sy values resulting from the history-matching process generally are in good agreement with values reported for similar sediments in the literature (Batu, 1998).

Boundary Conditions

We parameterized the principal model boundary conditions including specified flux (well pumping and recharge) and specified head (constant head), using multipliers that were constant for all stresses in a given stress period but varied with time between stress periods. These parameters were not constrained with autocorrelation because errors were not interpreted and uncertainty in their values would be associated from one time to another. For example, recharge might be adjusted in a stress period with a storm that occurs near the end of a stress period such that the effect of that storm is more appropriately applied to the model in the next stress period. In such a case, one stress period would be adjusted to decrease, whereas the next would be adjusted to increase recharge. The values of these multipliers across stress periods and across ensemble realizations for the texture and T–PROGS model insets are shown in figure 50. The limited range (20 percent) of flexibility is evident on these charts and, with a few exceptions, the multipliers did not reach their bounds.

Hydrologic Budget

The main components of the hydrologic budget for the texture model inset base realization are shown by model stress period in figure 51 and by calendar year in figure 52. The primary sources of groundwater inflow to the model include recharge (“RCHA In,” figs. 51–52), groundwater inflow across lateral boundary faces (“CHD In,” figs. 51–52), and through injection wells at recharge basins (“WEL In,” figs. 51–52). The primary outflows from the groundwater system were groundwater outflow across layer boundary faces (“CHD Out,” figs. 51–52), groundwater pumping at wells (“WEL Out,” figs. 51–52), and through discharge to gaining streams (“SRF Out,” figs. 51–52).

Storage terms are consistent with MODFLOW conventions where storage out is water entering storage and thus leaving the flowing groundwater system (STO–SS out and STO–SY out) and storage in (STO–SS In and STO–SY In) is water leaving storage and thus rejoining the flowing groundwater system (figs. 51 and 52; Langevin and others, 2017). Recharge was the largest source of groundwater inflow on an annual basis (fig. 52), though the volume of recharge inflow varied substantially by month (stress period), reflecting seasonal and event-driven dynamics in precipitation. Groundwater outflow across lateral boundaries was typically the largest annual groundwater sink. Well pumping was also a major outflow of groundwater and was often the largest outflow component during summer months (fig. 51).

Converting Model Stresses to Steady State

We calculated long-term average model stress and boundary values to create the steady-state inset model. Recharge and pumping rates, and constant head boundary heads that link the inset model to the transient LIRM, were averaged during 2005–19 to represent long-term average conditions. Water-supply pumping rates were adjusted to represent expected long-term average rates. Additionally, long-term average remedial well pumping rates (referred to by NYSDEC and hereafter in this report as the “commitments” pumping rates) represented the starting values for simulating the remedial system. With these changes in pumping, this steady-state inset model differs from the steady-state LIRM (Walter and others, 2020). Pumping rates for the commitments wells (the remediation wells for which commitments pumping rates have been assigned) were adjusted from these starting values to explore the tradeoffs of multiple objectives while maintaining constant pumping at the water-supply wells. The term “commitments” is used interchangeably with “baseline” in subsequent discussion to represent the initial assumed pumping regime from which optimization analysis starts.

Running Monte Carlo on Nontransient Estimated Parameters

The history-matching process described in the “Parameter Estimation by Ensemble History Matching” section resulted in two ensembles of model parameters—one each for the texture model inset and the T–PROGS model inset. These include constant parameter values for hydraulic conductivity and streambed conductance, and time-varying parameter values for recharge and constant head boundaries. The transient history-matching model runs include a steady-state initial stress period, so transient parameter multipliers for the constant head boundaries and recharge were adopted from that initial stress period. The storage parameters are not applicable to the steady-state model and thus are disregarded.

The performance of the ensemble as described above in the history-matching effort is based on correspondence between transient model outputs and the observation dataset. However, for evaluating remedial system performance, a different set of metrics are needed to explore the uncertainty of model outputs and to determine a risk with respect to optimization. The metrics of interest to the NYSDEC are metrics of capture by remedial wells and water-supply wells.

We added a particle-tracking algorithm, MODPATH 7 (Pollock, 2016) to the steady-state MODFLOW 6 model to evaluate plume capture. Forward trajectories of particles were simulated as starting in all model cells that intersect the plume shell (fig. 2). At the end of the simulation, the final destination (that is, endpoint) of the 113,328 particles was used to calculate the fraction of particles ending in various locations including the pumping wells or other model boundaries. The entire posterior Monte Carlo generated parameter ensemble was used for this particle-tracking analysis to provide an evaluation of the particle-capture uncertainty.

Selection of a Risk Level for Optimization

The posterior Monte Carlo analysis resulted in a range of values for metrics of capture by various classes of wells. The range of particles captured by commitment remedial wells (as a percentage of all particles released from the footprint of the plume shell) are shown on figure 53. The median of percent particles captured represents the risk-neutral stance. Only one realization from the ensemble of parameter sets is consistent with that median behavior. All realizations indicating percent capture less than the median can be considered risk-averse in the sense that the uncertainty of the results indicates a potential underestimate of capture. Conversely, realizations indicating percent capture greater than the median can be considered risk-tolerant in that they represent realizations that may overestimate capture. Evaluation of the uncertainty in this way is a form of correcting the model outputs considering uncertainty through “chance constraints” (Wagner and Gorelick, 1987). The risk-neutral stance selected for this analysis required the use of the realization that resulted in particle capture by the remedial wells closest to the median for each of the conceptual models (texture model and T–PROGS model insets).

Parameter Sweep

The initial optimization setup included adjusting remedial well pumping rates to maximize total plume capture. We initially started simulations with no remedial pumping and then adjusted pumping rates upward in increments of 5 percent of the total baseline or commitments remedial pumping, applied to all remedial wells equally, as much as 200 percent of the baseline remedial pumping. The increase in total capture is only from 80 percent to nearly 100 percent for the texture model inset and 95 percent to nearly 100 percent for T–PROGS model inset (fig. 54). However, as remedial pumping decreases from the baseline (100 percent of commitments) pumping, water-supply well capture increases owing to the location of several water-supply wells within the plume shell boundary as the remedial well capture decreases. In other words, adjusting remedial pumping substantially affects redistribution of capture between water-supply and remedial wells but only slightly affects total capture. Therefore, using total capture as an objective for optimization would likely result in an unacceptable level of water-supply well capture of the plume, which highlights the importance of evaluating system behavior prior to running through a sophisticated analysis.

Formal Nonlinear Multiobjective Management Optimization

The tradeoffs between remedial and water-supply well pumping and plume capture necessitate a more complicated analysis than traditional single-objective optimization where a single objective is to be minimized or maximized, such as water-supply well and remedial well plume capture. Fienen and others (2024) provide a comparison between single and multiple objective optimization for this case. Instead, exploring the tradeoffs between these potentially competing objectives using a multiobjective approach may be more appropriate. In this analysis, the Non-Dominated Sorting Genetic Algorithm II (NSGA–II) (Deb and others, 2002) as implemented in PEST++ through the Multiple Objective Optimization under Uncertainty (MOU) version (White and others, 2022) was used.

The NSGA–II is a genetic algorithm requiring multiple generations of computations to converge on an optimal set of solutions that traverse the alternatives that are not dominated by any other alternative (a pareto frontier; Goodarzi and others, 2014). Solutions on this pareto frontier are defined ideally as points that are optimal for all objectives. However, exact pareto optimality with respect to all objectives is not achievable, so the frontier solutions are actually those for which only one objective can be improved.

Defining Decision Variables

“Decision variables” in management optimization are defined as the model inputs that are adjusted to explore optimal model solutions with respect to the multiple objectives and are subject to hard constraints. The only decision variables in this analysis were pumping extraction rates in remedial wells.

Simulating the Remedial System

The remedial system is defined by remedial pumping wells, treatment systems, injection wells, and infiltration basins at which treated water is returned to the subsurface. Extraction rates (decision variables) are selected by the algorithm to simulate each candidate solution. The amount of water extracted by each well then must be simulated as conveyed through the appropriate treatment system. The pumped water in each treatment system is then returned to the subsurface through the appropriate infiltration basin or injection well. The conveyance of extracted water to specific treatment systems is outlined in table 13 and shown schematically in figure 55.

Enforcing Constraints

Hard constraints are a combination of functions of model inputs and outputs and can restrict the search space for decision variables. For this analysis, the constraints enforced on model inputs included the following:

These constraints can be operational (for example, limitations on treatment system capacity) or more subjective to set limits on the feasible search space of the algorithm (such as setting an upper limit on particle capture by supply wells, which was determined using early versions of the optimization).

Defining Objectives

Five potentially competing objectives that the algorithm was intended to optimize were defined. The objectives were the following:

Analysis of Optimal Solutions for Plume Remediation

The MOU algorithm was started with 800 initial randomly selected sets of decision variable values (that is, individuals in a population). This population was allowed to progress for 70 generations in which various attributes of the individuals are combined and their fitness is evaluated. Fitness is defined as metrics in line with the objectives (for example, higher total particle capture and lower remedial pumping are both considered to have higher fitness).

The solutions that are chosen for management are “optimal” in the sense of having the greatest “fitness” as described above. However, particularly with multiple competing objectives, the determination of “optimal” is subjective. For the discussion here, we chose representative solutions as optimal for the purpose of discussion, but other analysts may judge the tradeoffs differently and choose different optimal solutions. The final set of optimal solutions, chosen for the purpose of discussion, for the texture model and T–PROGS model insets, are shown in figure 56A and B, respectively. The final set of optimal solutions for all the model simulations that were evaluated throughout the process are shown in figures 56C and D for the texture model and T–PROGS model insets, respectively. There is a nearly linear tradeoff between water supply and remedial capture, with the maximum capture by remedial wells corresponding with the lowest water-supply capture (fig. 54) in figure 56A and B; however, this selection does not consider the cost of remedial pumping or the costs of treatment that may be required at water-supply downgradient wells.

A key advantage to using multiple objectives is that optimal solutions can be refined by incorporating additional information. For example, figure 57 shows the same information as figure 56, but the y-axes are changed from remediation capture to total remediation pumping. Remedial pumping is expressed as a percent of the baseline pumping rates where values less than 100 indicate a decrease in overall pumping and values greater than 100 indicate an increase in overall pumping. In contrast with the relation between water-supply and remedial capture, which is largely linear, the tradeoff is nonlinear when comparing water-supply capture with total remedial pumping.

The shape of the curves in figure 57 is less linear than in figure 56 and is characteristic of a pareto frontier. Near the origin of each graph, improving one objective degrades another objective—this is referred to as a pareto optimal point. For example, the decrease in water-supply plume capture from 21 to 18 percent results in an increase in remedial pumping of nearly 30 percent. From another perspective, increasing remedial pumping from 90 to 120 percent of the baseline only represents a 2-percent change in water-supply capture.

However, the straighter edges of the curve distant from the origin indicate locations where improvement can be made in one objective without incurring cost in the other. For example, there is a potential outcome where the water-supply capture can be decreased from 27 to 21 percent with only a small increase in total pumping (fig. 57A) because the MOU algorithm is able to redistribute pumping among the wells to optimize this outcome without significantly increasing the total pumping.

The previous discussion presented the pareto optimal point example when only two competing objectives were considered. A third objective—downgradient capture—is also shown in figures 56 and 57. All three objectives (water-supply capture, remedial pumping, and downgradient capture) are observed at low values near the origin where a trough of downgradient capture is present when all the model runs are evaluated along the path to optimality (fig. 57C and D).

The pumping rates for two selected solutions from the optimization are shown in figures 58 and 59. The spatial distribution of pumping for the solutions were chosen (figs. 58B–C and 59B–C, respectively) based on the criterion of maximizing remedial capture without regard to total remedial pumping or downgradient capture (fig. 58), and the tradeoff point is identified in figure 56. A comparison of the optimal solutions with the baseline indicates slight increases in remedial well pumping along the outer edge of the plume. The subtlety of these changes largely is due to the strict constraints on pumping that allowed for limited redistribution of pumping among a subset of the remedial wells.

When total remedial pumping cost and downgradient plume capture were considered (fig. 59), the optimal solution was identified as closest to the pareto optimal point indicated in figure 57. The subtle changes from the baseline needed to achieve optimal pumping are similar to changes shown in figure 56 but with pumping generally redistributed to focus more on the center of the plume, which allows for more water-supply capture to the east in favor of preventing further downgradient migration. Baseline and optimal pumping rates for the two optimal solutions for each model are listed in table 14.