Spatial and Temporal Discretization

The finite-difference grid for the numerical model consists of 145 rows, 245 columns, and 21 layers. The model grid size is variable, with a grid spacing of 0.25 mi² (the smallest) in the northern and southwestern parts of the New Jersey Coastal Plain, 0.31 mi² in the southeastern part of the Coastal Plain, and as much as 3.16 mi² in the offshore areas. The model area is approximately 16,597 mi² and extends from the Fall Line (fig. 1), which is the transition from the consolidated rocks of the upland Piedmont Physiographic Province to the unconsolidated sediments of the Coastal Plain Physiographic Province, to the approximate downdip and (or) freshwater-saltwater interface in each model layer. (The location and extent of the Fall Line in southeastern New York is not shown in figures.) The grid is oriented approximately parallel to the Fall Line and to the strike of the New Jersey Coastal Plain hydrogeologic units (Zapecza, 1989).

A schematic diagram of the model layers used to represent the hydrogeologic units in the New Jersey Coastal Plain and the corresponding hydrogeologic units in Delaware and eastern Maryland is shown in figure 4. Confining units are modeled as individual layers by assigning them appropriate hydrologic properties.

The model area extends vertically from the water table to consolidated bedrock and includes 11 distinct major aquifers and 10 confining units of the New Jersey Coastal Plain. The hydrogeologic framework described by Zapecza (1989) was revised for this study to improve the simulation of the groundwater-flow system. Extents of each model layer in Delaware and eastern Maryland are based on hydrogeologic data from Meisler (1989), Vroblesky and Fleck (1991) and Andreasen and others (2013). The hydrogeologic framework is described in more detail in the data release for this report (Carleton and others, 2023). This data release contains an ArcGIS shapefile containing the top, bottom, thickness, and extent of each layer in the model used to define the aquifers and confining units incorporated into the New Jersey Coastal Plain model (Carleton and others, 2023).

If a hydrogeologic layer does not extend across a given model layer in the MODFLOW-2005 finite-difference solution method but instead the hydrogeologic layer pinches out, but the model layer above and (or) below is an active area in the model, an area was created in the model layer with the limited extent to make that layer continuous with the above and (or) below layer, if needed, and was given a small thickness (1 foot or more). The hydraulic properties of the overlying model layer were assigned so that the layer would be continuous, thereby achieving a hydraulic connection between the layers above and below. For example, model layer 11 which represents the Vincentown aquifer, an aquifer with a small area, but the areas of the overlying and underlying layers extend beyond the area of the Vincentown aquifer, therefore an area is needed to simulate the flow from overlying aquifers to underlying aquifers to more accurately represent the vertical head distribution in the model layers (fig. 1). Additionally, in the updip area of the Piney Point aquifer (model layer 9) where the aquifer pinches out and does not outcrop, an area is needed to simulate the flow from overlying units to underlying units to more accurately represent the vertical head distribution in the model layers (fig. 1). The same is true for the updip area in model layer 5, which in the downdip area (fig. 1) represents the confined Rio Grande water-bearing zone, and the overlying and underlying confining units (model layers 4 and 6, respectively). In the updip parts of those three model layers, an area is needed to simulate the flow from the overlying model layer 3 to the underlying model layer 7 to more accurately represent the vertical connection in the model layers.

A total of 36 stress periods were simulated—1 steady-state period and 35 transient periods. Stress period 1 is steady-state and stress period 2 is transient and 3 years long. Stress periods 1 and 2 are transitional stress periods that generally correspond to 1970–76 and 1977–79, respectively. Stress periods 3 to 36 are 1 year long and generally represent the period 1980–2013.

Lateral and Lower Boundary Conditions

The model is bounded to the north and northwest and below by the contact between the unconsolidated Coastal Plain deposits and consolidated bedrock that underlies the Coastal Plain sediments. This boundary was represented numerically as a no-flow boundary condition. Lateral boundaries for each layer in the groundwater-flow model are shown in a generalized schematic diagram in figure 4. Although water may flow into the model area along the Fall Line between the subsurface fractured rock and the Coastal Plain sediments, this amount is assumed to be negligible and, therefore, is not considered in this report.

In the two previous New Jersey RASA groundwater-flow models (Martin, 1998; Voronin, 2004), the lateral boundaries in the northeast and southwest are specified-flux boundaries that originally were derived from a model of the NACP (Leahy and Martin, 1993). Some modifications were made to the boundary fluxes near Delaware Bay in the rediscretized RASA model (Voronin, 2004). The Leahy and Martin (1993) model area extended from Long Island, New York, to North Carolina and includes the New Jersey Coastal Plain. In the model described in this report, the lateral model boundaries in the northeast and southwest are specified-flux boundaries derived from the NACP aquifer system model (Masterson and others, 2016a), which also extends from Long Island, New York, to North Carolina, by using the MODFLOW Flow and Head Boundary (FHB) package (Leake and Lilly, 1997). Specified fluxes are applied in areas where the aquifers and confining units intersect the flow boundary. Flows from the NACP model (Masterson and others, 2016a) are incorporated into the FHB package for stress periods 1 through 32. Because the end of the NACP model simulation coincides with stress period 32 of the current New Jersey Coastal Plain model, fluxes for stress period 32 are repeated for stress periods 33 through 36. The FHB package is used to specify lateral boundary fluxes in each active cell along column 1 of the model grid, which traverses eastern Maryland and a small part of southwestern Delaware (fig. 1). Some active cells along column 1 where the flux was zero were not included in the FBH model file. Specified fluxes are also input in New Jersey using the FHB package along the northeastern edge of the active model area for some layers in column 245 for the offshore extent of layers 7 and 8 and in layers 16 to 19 for the offshore part of the upper and middle Potomac-Raritan-Magothy aquifers (layers 17 and 19), the Merchantville-Woodbury confining unit (layer 16), and the confining unit between the upper and middle Potomac-Raritan-Magothy aquifers (layer 18). The remaining layers have a no-flow or a constant-head boundary along northeastern limit of the model area (fig. 4).

The southeastern (downdip) boundary of the aquifers and confining-unit model layers in New Jersey, Delaware, and Maryland are simulated as no-flow boundaries except as described in the “General Head Boundary” section of this report. The downdip extents of the aquifers and confining-unit model layers in Delaware and Maryland are from Meisler (1989), Vroblesky and Fleck (1991), and Andreasen and others (2013). The southeastern downdip boundary of the aquifers and confining-unit layers in New Jersey are generally located at the downdip limit of freshwater in the aquifer (10,000 mg/L isochlor) determined by Pope and Gordon (1999). Pope and Gordon (1999) used the SHARP computer model (Essaid, 1990) to simulate the groundwater-flow system, including the location and movement of the freshwater-saltwater interface in nine aquifers and eight intervening confining units in the Coastal Plain. Although the interface position is not static in the model, the actual amount of movement cannot be quantified because of the model’s coarse cell size. Because the SHARP model is quasi-three-dimensional, the downdip limit of the 10,000-mg/L isochlor in those aquifers was assumed to be the same for the overlying confining unit. However, the Englishtown aquifer system and Vincentown aquifer are not continuous throughout the New Jersey Coastal Plain and contain freshwater throughout their confined extent in New Jersey. However, chloride concentrations exceeded 15,000 mg/L at an observation well in the Englishtown aquifer system in Sandy Hook in northeastern Monmouth County (DePaul and Rosman, 2015).

General Head Boundary

The General Head Boundary (GHB) package (McDonald and Harbaugh, 1988) was used to represent flow in the updip direction from limited parts of the downdip lateral boundaries of layers 7, 9, 11, 13, 15, 17, 19, and 21 (fig. 4). However, in model layer 11, the GHB boundary was used in New Jersey only. The GHB head is constant for all stress periods and is set to 0.0 ft (NAVD 88); no attempt was made to estimate the equivalent freshwater head where groundwater density is greater than that of freshwater. A single value of bed conductance of 1 square foot per second (ft²/s) was applied to all GHB cells.

The GHB boundaries were included to prevent the occurrence of unrealistically low heads downdip from withdrawals; because the boundaries are generally far from the withdrawals, they have only a minor effect on the overall flow system. In 2013, flow into the model from GHB cells is about 0.6 percent of the total flow budget.

Upper Boundary Conditions

The upper boundaries in the model represent the streams and recharge in onshore areas and the Delaware and Raritan Bays and Atlantic Ocean (fig. 1) in offshore areas. In offshore areas, the upper boundary is a constant freshwater equivalent water level specified by using the Time-Variant Specified-Head (CHD) package (Harbaugh, 2005) and is applied across the top active layer of model cells used to represent the Atlantic Ocean and the Delaware and Raritan Bays. CHD cells represent the outcrop of layers 2 through 5 in the Delaware Bay, the outcrop in layers 14 through 19 (excluding layer 18) in Raritan Bay, and the outcrop of layers 2, 3, 7, and 10 through 12 in the Atlantic Ocean (fig. 4).

Recharge

The Recharge (RCH) package of the USGS finite-difference groundwater model MODFLOW–2005 was used to apply recharge to the topmost active layer of each onshore model cell. This upper boundary is spatially variable recharge from precipitation estimated by using the Soil-Water Balance (SWB) model (Westenbroek and others, 2010). The SWB model is based on a modified Thornthwaite-Mather method (Thornthwaite and Mather, 1957) that incorporates spatially distributed land cover, soil properties, and daily climate data to produce a detailed spatial and temporal variability of recharge. SWB model-calculated recharge rates from 1980 to 2013 were used as the initial recharge values for New Jersey in the groundwater-flow model. The SWB model is described in appendix 1. Recharge was adjusted as part of the model parameter-estimation process to allow for more variability than that calculated by the SWB model and adjusted for comparison with base-flow estimates at various surface-water drainage basins in the New Jersey Coastal Plain. Average annual precipitation in southern New Jersey (Office of the New Jersey State Climatologist, 2021) during 1970–76 and 1977–79 were 1.23 and 1.31 times the 1980 annual average, respectively. Therefore, the recharge rates applied in model stress periods 1 (associated with 1970–76) and 2 (1977–79) are 1.23 and 1.31 times those applied in stress period 3 (1980). Stress periods 1 and 2 are used in the model to set antecedent conditions for the other years.

Because the surficial aquifer in Delaware and Maryland is not simulated, cell-by-cell vertical flows into or out of active cells in model layer 1 in Delaware and Maryland were calculated from the NACP model (Masterson and others, 2016a) and applied with the recharge.

Drain Package

The MODFLOW Drain package (DRN) is used to represent the rivers and streams in the model where groundwater flows (drains) out of the model area (fig. 5). The DRN package was used for continually gaining stream reaches in the New Jersey Coastal Plain which are not a source of groundwater recharge except on a localized scale (for example, upstream from a dam). Use of the DRN package allows for the representation of headwater reaches, where the streambed elevation may be above the water table in places, preventing simulated small streams from being a source of potentially unlimited volumes of water, as would be the case with the MODFLOW River package (RIV) (discussed below). The drain elevation is the average elevation of a stream or stream tributary in the topmost active layer in each onshore model cell. Streams in New Jersey are those used in the previous model version from Voronin (2004). Model cells containing streams in New Jersey were identified by intersecting the model-grid polygons with streams digitized from USGS 1:24,000-scale topographic maps. Streams in Maryland and Delaware were identified by intersecting the grid with NHDPlus (U.S. Geological Survey, 2005), a comprehensive geospatial dataset developed by the U.S. Environmental Protection Agency and the USGS.

The elevations of drain boundaries are kept constant for all stress periods in the transient model. The drain elevations are long-term averages and do not include annual or seasonal variations. Seasonal variations cannot be considered when annual stress periods are used; neglecting the variation of average elevation from year to year is minor in the context of surface-water elevations averaged over model cells of 0.25 mi² or larger.

The DRN package also requires specification of streambed conductance, an aggregate parameter that represents the product of the streambed hydraulic conductivity, the length of the stream reach, and the stream width divided by the streambed thickness. The same streambed conductance is applied to all drains and was adjusted during calibration.

River Package

The MODFLOW River package (RIV) package was used to represent the Delaware River (fig. 5). The RIV package (unlike the DRN package) allows the Delaware River to be a source of recharge to the aquifer near pumping centers. This is appropriate because the Delaware River carries substantial flow, is tidal in the simulated area, and can provide unlimited flow to the aquifer. The river stage is constant for all stress periods and increases from 0.0 ft (NAVD 88) at the southwestern extent of model river cells in Delaware Bay to 1.0 ft near Trenton by 0.1 ft in segments of approximately the same length. A single value of riverbed conductance was applied to all river cells, representing the product of the riverbed hydraulic conductivity, length, and width of the river reach in each model cell divided by riverbed thickness. Riverbed conductance was modified during model calibration.

Hydraulic Properties

The hydraulic properties hydraulic conductivity and storage coefficient are described in the “Hydraulic Properties of Aquifers and Confining Units” section near the beginning of this report. As previously mentioned, data for the hydraulic properties of confining units are limited. The hydraulic properties are specified by parameter zones within each layer that are assumed to have similar properties. The initial values are comparable with those from Martin (1998), Voronin (2004) and Masterson and others (2016a). During model calibration, the range of hydraulic properties shown in tables 2 and 3 were used to provide supplemental information for values in areas of New Jersey. The hydraulic conductivities and storage coefficients in table 4 were used to provide supplemental information for values in Delaware and Maryland. Hydraulic-conductivity and storage-coefficient parameters were adjusted manually and by using the parameter-estimation software UCODE-2014 (Poeter and others, 2014) during model calibration as needed to improve the model fit between observed data and simulated results.

Groundwater Withdrawals

Groundwater-withdrawal data for the New Jersey Coastal Plain were tabulated and mapped to assess the volume of water pumped from wells in each of the aquifers. These wells include those used for public supply, large-scale agricultural (irrigation), and commercial or industrial purposes. Withdrawals from small-capacity wells, such as those used for domestic supply, are not incorporated into the model. This is a limitation of the model in the outcrop areas of confined aquifers or the Kirkwood-Cohansey aquifer system in areas where a substantial number of domestic-supply wells are located, but domestic self-supply withdrawals do not make up a substantial part of the groundwater withdrawals. Statewide production from domestic wells accounts for about 3 percent of total withdrawals in New Jersey (New Jersey Department of Environmental Protection, 2017b).

Withdrawal data for 1980–1998 were obtained from Voronin (2004) and were reviewed for updates that may have been made since 2003, such as an update to an aquifer code or surveyed well location. Withdrawal data for 1999–2003 were obtained from data reported to the NJDEP. These data were quality reviewed and incorporated into a water-use database available in the data release accompanying this report (Carleton and others, 2023). This withdrawal information includes permitted data only—that is, only data from wells in which daily withdrawals meet or exceed 100,000 gallons for a period of more than 30 days in a consecutive 365-day period. Withdrawal data for 2004–13 were obtained from the NJWATr database developed by the USGS and maintained by the NJDEP to track water withdrawals, use, treatment, and discharge in New Jersey (New Jersey Department of Environmental Protection, 2017a). Total withdrawals from New Jersey Coastal Plain aquifers input to the model for 2004 to 2013 are shown in figure 6. The largest withdrawals occurred in the unconfined Kirkwood-Cohansey aquifer system, the Atlantic City 800-foot sand, and the upper, middle, and lower Potomac-Raritan-Magothy aquifers (model layers 3, 7, 17, 19, and 21).

Annual withdrawal data for 1980–2013 for Delaware and Maryland were obtained from a dataset of pumping rates for 1980–2010 from the groundwater-flow model developed for a groundwater-availability study of the NACP by Masterson and others (2016a; Jack Monti, Jr., U.S. Geological Survey, written commun., 2014). Agricultural withdrawal data for Delaware and Maryland may be less complete and less well documented (Masterson and others, 2016a) than New Jersey agricultural withdrawals. Masterson and others (2016a) provide a detailed description of the methods used to tabulate and estimate groundwater withdrawals in their study. Withdrawals from the unconfined surficial aquifer in Delaware and Maryland were not included in the New Jersey Coastal Plain groundwater-flow model. The unconfined surficial aquifer in Delaware and eastern Maryland was not included in the model hydrogeologic framework because model layer 1 represents the Holly Beach water-bearing zone in Cape May County, New Jersey, which is not continuous into Delaware.

Stress periods 1 (steady-state) and 2 (transient) are transitional stress periods that generally correspond with 1970–76 and 1977–79, respectively. Stress periods 3 and 4 generally represent 1980 and 1981, respectively. Withdrawals for stress periods 1, 2, and 3 are 73, 86, and 93 percent, respectively, of those in stress period 4, and account for the estimated increase in withdrawals from 1970 to 1980. Withdrawals for stress period 4 generally are an average of 1978–81 withdrawals in Voronin (2004) and represent 1981 withdrawal rates in the groundwater-flow model. Average annual withdrawals were used in the model for stress periods 5 to 36, which are 1 year long, and the average annual withdrawals are equivalent to total annual withdrawals during each of these stress periods. Stress periods 5 through 36 represent each year from 1982 to 2013, respectively.

Water-Level Observations

The model was calibrated by using water-level observations that include observations derived from selected water-level measurements described below. The water-level observations in New Jersey are a subset of measurements collected during the seven quinquennial New Jersey Coastal Plain water-level synoptic studies (1983, 1988, 1993, 1998, 2003, 2008, and 2013) and were selected to provide a spatial distribution in 10 of the simulated aquifers. The water-level synoptic studies did not include the unconfined Holly Beach water-bearing zone and the unconfined part of the Kirkwood-Cohansey aquifer system. Water levels in correlated aquifers in Delaware were obtained from the Delaware Geological Survey (Delaware Geological Survey, 2017). Most of the water levels were measured during October–December after large summer withdrawals had decreased but before water levels had recovered to those measured during spring, thereby representing an approximation of the annual average water level. Water levels were measured at least once during 1983–2013 in 392 wells in New Jersey and 48 wells in Delaware (fig. 7), yielding a total of 2,533 water-level observations in 440 wells. Hydraulic water levels at each well were calculated by subtracting the water level, in feet below land surface, from the land-surface altitude, in feet above the North American Vertical Datum of 1988 (NAVD 88); the result of this calculation is shown in this report as the water-level altitude.

Average annual water levels were calculated for 29 wells in New Jersey (fig. 7) where continuous water levels were measured: 19 wells with observations for every year during 1980–2013 and 10 wells with data for most of that period (average of 27 years). A total of 920 observations from the 29 wells were included in the model calibration; these data are shown in time-series hydrographs shown in this report in the “Calibration of Water Levels” section.

In addition to the water-level observations used in the model calibration, derived observations also were included by calculating the vertical gradient at nested observation wells in New Jersey. Vertical gradients were calculated for 33 pairs of nested wells, for a total of 210 observations. (Nested well pairs are wells that are near each other but are open to different aquifers.) For each sequential synoptic water-level measurement period, vertical gradients calculated by subtracting water levels in colocated wells screened in different aquifers provided information for calibrating the hydraulic conductivity of the confining units. By attempting to match the gradient across the confining unit, the parameter-estimation function improved the accuracy of the vertical-conductance value, for example, in a situation where simulated water levels in both aquifers were somewhat high but the gradient was exact. Also, the change in water level over time was calculated by subtracting water levels in selected wells in which measurements had been made in sequential water-level synoptic studies; the difference was used in calibrating storage in aquifers and confining units. Changes in water levels over time were calculated for 138 wells, 134 wells in New Jersey and four wells in Delaware, where water levels had varied substantially (approximately 10 ft) over the 30-year span of synoptic water-level measurements, for a total of 767 observations.

Water-level observations and vertical-gradient and change-over-time observations calculated from water-level measurements were all assigned the same weight, except for annual observations from the 29 wells with a continuous water-level recorder. Most water-level hydrographs for selected wells had 34 observations, whereas many synoptically measured wells had at most 7 observations; therefore, synoptic water-level observations and calculated observations were assigned a weight of 5, and hydrograph observations were assigned a weight of 1.1. Hill and Tiedeman (2007) describe methods for assigning weights on the basis of estimated error ranges so that observations with substantial potential errors exert a smaller effect on the parameter estimation function than the more accurate measurements do. This approach was not practical for the model in this study because water-level observations were collected under a wide range of conditions, often unknown. For example, some water-level measurements were made in wells with varying degrees of recent and nearby pumping, whereas others were made in unused wells known to be thousands of feet from current or recent withdrawals; still other water-level measurements were made at production wells that had been pumped as recently as 1 hour earlier. Attempting to quantify potential errors based on proximity in time or space to known and unknown withdrawals for synoptic water-level measurements conducted over a 30-year period was infeasible. Most land-surface-elevation estimates were determined from light detection and ranging (lidar) data considered to be accurate to within 2 ft, a modest error compared to the uncertainties related to proximal withdrawals at the time of measurement.

Base-Flow Observations

The New Jersey Coastal Plain model was calibrated by using a total of 1,485 estimates of base flow from 47 stream basins in New Jersey. The locations of the 47 basins are shown in figure 8. The surface-water basins were selected because of location and availability of data during the period 1980—2013. Annual base flow was calculated by using the hydrograph-separation technique PART in Barlow and others (2015) for each of the 21 continuous-record gaging stations for each year of record. For 26 partial-record gaging stations in the New Jersey Coastal Plain, an annual base flow was calculated by correlating manual low-flow measurements made at the partial-record stations with same-day flows at index stations by using the MOVE.1 program (Colarullo and others, 2018). Annual base-flow estimates for each of the 47 basins were averaged to obtain a period-of-record average base flow for each basin.

Base-flow observations determined by using the PART hydrograph-separation technique in Barlow and others (2015) of continuous-record-station data were assigned a weight of 1.0. Observations determined from partial-record-station data correlated with measurements from continuous-record index stations were assigned a weight less than 1.0 based on the percent standard error of estimate (%SEE), which is the percent of unexplained variance relative to the value of the estimate. The %SEE is a better measure of uncertainty than the standard error of estimate (SEE) because it accounts for the proportional effect encountered in natural streamflow processes, where variance increases with flow (Amy McHugh, U.S. Geological Survey, written commun., 2020). Colarullo and others (2018) defined weights as the reciprocal of %SEE to lessen the weight of predictions made using MOVE.1 regressions that show a high degree of unexplained variability in flows and increase the weight of predictions made using regressions with a small amount of unexplained flow variability. The weights of measurements made at partial-record stations were set to a number less than 1.0 by dividing 4.5 by the squared average %SEE for each partial record station; resulting weights ranged from 0.95 to 0.22.

Model Parameters

A parameter is defined in this report as a single value assigned to a variable used in the finite-difference groundwater-flow equation at one or more model cells. The 422 parameters used in the model fall into seven categories: horizontal hydraulic conductivity (184 parameters), vertical hydraulic conductivity (184 parameters), storage coefficient (42 parameters), riverbed conductance (1 parameter), streambed conductance (1 parameter), general-head-boundary conductance (8 parameters), and recharge (2 parameter). The initial parameter values for several aquifers and confining units in New Jersey, Delaware, and eastern Maryland were changed during calibration to improve the match with observations and improve model results.

The final calibrated values for horizontal and vertical hydraulic conductivity and storage coefficients for aquifers and confining units are listed in tables 6 through 9. Table 6 lists the horizontal hydraulic conductivity (K_h) for the aquifer layers, table 7 lists the horizontal hydraulic conductivity for the confining-unit layers, and table 8 lists the vertical hydraulic conductivity (K_v) for the model layers. Table 9 lists the storage coefficients for each model layer. Parameter zones in each model layer were created using these values to allow for spatial variability of K_h and K_v. Several parameter zones were added or modified during model calibration. The conversion to a fully three-dimensional model framework and varied recharge resulted in the need for additional K_h parameter zones in the outcrop areas of all layers with updip outcrops following the strike of the hydrogeologic layers, corresponding to model layers 3 through 8 and 10 through 21. During the model calibration, some zones acquired a parameter value equal to or similar to that of an adjacent zone but remained a separate zone and was not incorporated into the adjacent zone. The ratio of horizontal to vertical hydraulic conductivity (K_h:K_v) was initially set to 10:1 for aquifers and 1:1 for confining units. K_h and K_v were adjusted independently during calibration; the sensitivity of K_v is low in most aquifer parameter zones and the sensitivity of K_h is low in most confining-unit parameter zones. The zones are designated by the parameters listed in tables 6 through 9. The location and extent of the K_h zones for the aquifers and K_v zones for the confining units are shown in figures 9 through 19. The parameter names for horizontal hydraulic conductivity for the confining units are zones with prefix “CU”. The location and extent of the horizontal hydraulic conductivity parameter zones for the confining units (table 7) are not shown in figures. The location and extent of the horizontal hydraulic conductivity parameter zones for the confining units coincide with the location and extent of the vertical hydraulic conductivity parameter zones (zones with prefix “VK”) for the confining units shown in figures 9B through 18B for those model layers.

The final parameters for horizontal and vertical hydraulic conductivity were compared to the values shown in tables 2 and 3 of this report. Many of the parameter values of horizontal hydraulic conductivity, such as those for the confined aquifers and the Holly Beach water-bearing zone in New Jersey, fall within the ranges of the value for horizontal hydraulic conductivity shown in the tables 2 and 3 and also table 1 of Masterson and others (2016a). Some horizontal hydraulic conductivity assigned to the outcrop areas of certain model layers for aquifers were higher (greater than 500 ft/d, table 6) because those areas are thin and the higher value increased model stability in these areas but did not affect the flow farther downdip. Higher values of horizontal hydraulic conductivity were also assigned to the outcrop areas of some confining-unit model layers (table 7) for the same reason. Lower values of horizontal hydraulic conductivity (less than 1 ft/d) were assigned in downdip areas of some aquifers—for example model layer 11 which represents the Vincentown aquifer—because the model layer is not present downdip and the small horizontal hydraulic conductivity value allowed for the lateral connection between the updip and downdip parts of the model layer and (or) a vertical connection between the model layer with the overlying and underlying model layers.

The parameters for riverbed conductance for the Delaware River cells and for streambed conductance for drains cells each have a single value. Although hydrogeologic and stream-width data could be used to justify multiple parameters for these properties, these data were considered secondary to the calibration of the regional-scale Coastal Plain model. During model calibration, simulation results indicated the model was insensitive to the riverbed-conductance parameter. Additionally, the comparison of simulated to estimated base flows was reasonable when a single streambed-conductance parameter for drain cells was used; therefore, neither parameter was subdivided to represent smaller or local areas for the river or drain cells. This process is described further in the “Parameter Sensitivity and Correlation” and “Parameter Estimation and Residual Analysis” sections below.

Two parameters were used to vary the recharge rates in the model calibration. One parameter was used for New Jersey and another parameter was used for recharge in Delaware and Maryland. The recharge parameter for New Jersey is a multiplier applied equally to the recharge value estimated for each model cell by using the Soil-Water-Balance model (described in the “Recharge” section earlier in this report). The recharge parameter was increased from an initial value of 1.0 to a final value of 1.2 to provide a reasonable match between simulated and estimated base flow in the 47 New Jersey surface-water basins (fig. 8).

Parameter Sensitivity and Correlation

A sensitivity analysis of the calibrated model was conducted to evaluate the relative effects of the various parameters on the match between simulated and observed water levels and simulated and estimated base flows. The sensitivity analysis is used to identify which parameters in the model, when modified, improve model results. Parameters which are least sensitive, when modified, have little effect on model results. The sensitivity of model parameters was calculated by using the parameter-estimation software UCODE–2014 (Poeter and others, 2014). Parameter sensitivities were calculated two ways: by using base-flow observations only and by using water-level observations only. This section discusses composite scaled sensitivities (CSS) and parameter correlation coefficients. The CSS values show how sensitive water-level and base-flow observations in the model are to changes in each of the parameters.

According to composite-scaled sensitivity analysis, the parameters most sensitive to base-flow observations (excluding recharge) are primarily streambed conductance, horizontal hydraulic conductivity, and storage in zones in the outcrop areas of several model layers (fig. 20). Recharge is the dominant parameter with respect to base-flow observations, and the recharge multiplier parameter was adjusted manually such that average simulated base flow equaled the estimated base flow of the New Jersey Coastal Plain of about 16 inches per year (in/yr), a value determined from surficial-aquifer studies completed in surface-water basins in the New Jersey Coastal Plain during 1992–2003 (Watt and Johnson, 1992; Watt and others, 1994; Lacombe and Rosman, 1995; Johnson and Watt, 1996; Johnson and Charles, 1997; Charles and others, 2001; Watt and others, 2003; Gordon, 200418). The streambed-conductance parameter for the drains (DRN_Par1) is ranked the 4th most sensitive parameter to base-flow observations (fig. 20). As indicated above, there is only one parameter for streambed conductance and any change to that parameter is applied to all streams uniformly rather than to individual streams. When the streambed conductance was decreased, groundwater levels tended to increase but the location or timing of discharge to base flow did not necessarily change. The 21 most sensitive parameters (fig. 20) other than streambed conductance (DRN_Par1) are in, or laterally or vertically adjacent to, an outcrop area.

Results of the CSS analysis indicate that the parameters most sensitive to water-level observations are recharge, and horizontal and vertical hydraulic conductivity (fig. 21); riverbed and general-head-boundary conductances are insensitive. Of the 27 most sensitive parameters (including the two recharge parameters (RCH_SWB and RchDeFlux), 19 (70 percent) are horizontal hydraulic conductivity, and 6 (22 percent) are vertical hydraulic conductivity.

Six pairs of the 130 most sensitive horizontal and vertical hydraulic conductivity parameters are correlated when recharge is excluded and only head observations are included in the parameter estimation (fig. 22). Most of the correlations are between vertically adjacent zones (for example, layers 20 and 21 in Delaware (figs. 18B and 19)) or horizontal and vertical hydraulic conductivity of the same zone (for example, layer 11 in northwestern Ocean County and northeastern Burlington County (figs. 1424.)). Four of the six pairs are parameter zones in Delaware and three of the six pairs are parameter zones in outcrop areas. Although the parameters were correlated, each correlated parameter was allowed to be estimated separately. There is not sufficient water-level data in the parameter zones to resolve the negative correlation between the horizontal hydraulic conductivity and (or) vertical hydraulic conductivity of each of the correlated parameter pairs.

Parameter Estimation and Residual Analysis

The parameters used in the New Jersey Coastal Plain model were tested and adjusted in an iterative calibration process. Initial parameter-estimation runs with UCODE–2014 (Poeter and others, 2014) were conducted by using all observations (water-level and base-flow). Results indicated that matching water levels in the aquifers would be less likely if base-flow observations were included. The challenge of adjusting the weights of base-flow observations to avoid their having too little or too great an effect on water levels was judged to be not effective because of the single streambed parameter for drains and the relative coarseness of the representation of the stream network in the DRN package. Therefore, after some adjustments were made to sensitive parameters to obtain an initial optimal calibration with all observations, parameter-estimation runs were done by using only base-flow observations. The recharge multiplier had the greatest effect on simulated base flows; once the recharge multiplier was held constant, the model response was generally more sensitive to parameters in the outcrop areas of aquifers and confining units. Subsequent parameter-estimation runs were done by using only water-level observations. For these runs, the parameters that were sensitive to drain (base-flow) observations, but not to water-level observations, typically remained at the value determined in the previous parameter-estimation runs in which only drain observations were used.

The number of simulated and water-level observations from wells in New Jersey used in the NJCP groundwater-flow model each total 4,243. The total includes 2,361 water-level observations from 392 wells, 920 historical water-level trends for 29 wells during 1980–2013, 210 observations of vertical gradients in 33 pairs of nested wells, and 752 observations of changes in water levels over time calculated for 134 wells. When plotted against each other (fig. 23), these values are clustered around the 1:1 correlation line where the observed water levels have a range of 411 ft and more than 96 percent of the simulated water levels fall within +/− 10 percent of the total range during 1980–2013. The residuals (simulated minus observed water level) are within + or – 10 feet for 58 percent of the water-level observations for New Jersey, and within + or – 20 feet for 84 percent. The average water-level residual is 1.5 ft; the minimum; first, second, and third quartile; and maximum residuals are −127.4, −5.5, 1.6, 10.4, and 179.4 ft, respectively.

Calibration of Water Levels in the Holly Beach Water-Bearing Zone

Simulated water-level contours in the Holly Beach water-bearing zone (model layer 1, Cape May County) are shown in figure 27. The Holly Beach water-bearing zone is a minor aquifer in New Jersey; the observed water levels shown in Lacombe and Carleton (2002) were used to visually compare to the simulated water-level contours.

Calibration of Water Levels in the Confined Cohansey Aquifer and Unconfined Kirkwood-Cohansey Aquifer System

Simulated water-level contours in the unconfined Kirkwood-Cohansey aquifer shown in the updip area of model layer 3 closely match the observed water-table contours in general direction and slope (fig. 28). The observed water-table contours are an aggregate of data from multiple surficial-aquifer studies in the following surface-water basins of the New Jersey Coastal Plain: Great Egg Harbor River Basin (Watt and Johnson, 1992); Toms River, Metedeconk River, and Kettle Creek Basins (Watt and others, 1994); upper Maurice River Basin (Lacombe and Rosman, 1995); Mullica River Basin (Johnson and Watt, 1996); Salem River and Raccoon, Oldmans, Alloway, and Stow Creek Basins (Johnson and Charles, 1997); and Rancocas, Crosswicks, Assunpink, Blacks, and Crafts Creek Basins (Watt and others, 2003). The aggregated (observed) water-table contours are shown only in the areas where a surficial-aquifer study was completed before 2004. In Cape May County, New Jersey, as well as in Sussex County, Delaware, groundwater is confined.

The updip limit of the confined Cohansey aquifer in Cape May County is shown in figure 29. The simulated water-level contours in the confined Cohansey aquifer in New Jersey and the potentiometric surface of the 2013 New Jersey Coastal Plain water-level synoptic study (Cauller and Gordon, 2021) for this aquifer also are shown in figure 29. Simulated water levels in the confined part of the Cohansey aquifer (model layer 3, Cape May County) generally match observed water levels (table 11). The average water-level residual (simulated minus observed) for 17 observation points in Cape May County for the confined Cohansey aquifer is −0.2 ft and the root mean square error (RMSE) is 3.3 ft. The simulated water levels for the Upper Chesapeake aquifer in Delaware and the average water-level residual for four wells measured in Delaware in 2013 (table 11) also are shown in figure 29.

The hydrographs of simulated and observed water levels for well 090049 (fig. 30) are similar, with minima, maxima, and averages for the observed water levels of –21, –8, and –13 ft, respectively, for the observed water levels and –19, –7, and –11 ft, respectively, for the simulated water levels. Simulated water levels are, on average, 2 ft greater than observed. The hydrographs of simulated and observed water levels for well 090080 in south-central Cape May County (fig. 31) are similar, with minima, maxima, and averages for the observed water levels of –11, –1, and –6 ft, respectively, for the observed water levels and –18, –0.5, and –8 ft, respectively, for the simulated water levels.

Calibration of Water Levels in the Rio Grande Water-Bearing Zone

Simulated water-level contours and the potentiometric surface of the 2013 New Jersey Coastal Plain water-level synoptic study (Cauller and Gordon, 2021) in the Rio Grande water-bearing zone (model layer 5, parts of Ocean, Atlantic, and Cape May Counties) are shown in figure 32. The 2013 simulated water levels are an average 7 ft lower than the observed water levels at six wells in Cape May County and are 35 ft lower than the observed water levels at three wells in Ocean County. The simulated 2013 potentiometric contours indicate a cone of depression along the shoreline of Cape May and Atlantic Counties caused by drawdowns in the underlying Atlantic City 800-foot sand; a cone of depression also is observed in the simulated water-level contours. The differences between the simulated and observed water levels for 10 observation points in New Jersey are listed in table 12. The average water-level residual (simulated minus observed) for the 10 observation points (figs. 32 and 33) is −16 ft and the RMSE is 24 ft. The minimum; first, second, and third quartile; and maximum residuals are −45, −31, −13, 0, and 5 ft, respectively.

Calibration of Water Levels in the Atlantic City 800-Foot Sand

Simulated water-level contours in the Atlantic City 800-foot sand (model layer 7, parts of Cape May, Atlantic, and Ocean Counties) generally match the 2013 potentiometric surface in coastal Atlantic and Cape May Counties (fig. 34). The simulated water-level contours are similar in shape to the 2013 potentiometric surface (Cauller and Gordon, 2021), although the gradient of the simulated contours is less steep than that of the water-level contours drawn from measurements made in updip central Atlantic County. The simulated water levels are within 10 ft of observed levels in 13 of 34 wells (38 percent) in New Jersey and within 15 ft of observed levels in 24 of 34 wells (71 percent) (figs. 34 and 35 and table 13). The 2013 simulated water levels are an average 8 ft higher than the observed water levels at 13 wells in Cape May County and are an average 13 ft lower than the observed water levels at 14 wells in Atlantic County. The simulated water levels are an average 11 ft lower than observed levels at 7 wells in Ocean County. The average water-level residual (simulated minus observed) for 34 wells in Cape May, Atlantic, and Ocean Counties is −4 ft and the RMSE is 21 ft. The minimum; first, second, and third quartile; and maximum residuals are −57, −12, 3, 12, and 17 ft, respectively. The simulated water levels for the Calvert aquifer system in Delaware and the average water-level residual for one well measured in Delaware in 2013 (table 13) also are shown in figure 34.

The hydrographs of wells 010180 and 010578, in northeastern and southeastern Atlantic County, respectively (fig. 36), show the same trend and amplitude as the observed water levels, but the simulated water levels are, on average, 1 ft lower and 12 ft higher, respectively, than the observed. The hydrographs of simulated and observed water levels for well 010180 (fig. 36A) are similar, with minima, maxima, and average for the observed water levels of –57, –35, and –47 ft, respectively, for the observed water levels and –60, –36, and –48 ft, respectively, for the simulated water levels. The hydrographs of simulated and observed water levels for well 010158 (fig. 36B) are similar, with minima, maxima, and average for the observed water levels of –71, –49, and –60 ft, respectively, and –60, –38, and –48 ft, respectively, for the simulated water levels. The hydrograph for well 010578 (fig. 36B) shows simulated water levels that are closer to and more closely match the trend in the observed water levels than those simulated by Voronin (2004) for 1990–98.

Calibration of Water Levels in the Piney Point Aquifer

Simulated water-level contours in the Piney Point aquifer (model layer 9, parts of Monmouth, Ocean, Atlantic, Gloucester, Cumberland, and Salem Counties, New Jersey, and Kent County in Delaware) are generally similar in shape to the 2013 potentiometric surface (Cauller and Gordon, 2021) (fig. 37), including the contours defining cones of depression in Kent County in Delaware and Cumberland and Ocean Counties in New Jersey. The simulated water levels extend farther updip in Burlington and Ocean Counties than the 2013 potentiometric surface because the updip limit of the Piney Point aquifer in the new hydrogeologic framework used in this groundwater-flow model used new well information from that area which was not available at the time of the previous framework developed by Zapecza, (1989). However, the 2013 potentiometric surface used the updip limit of the Piney Point aquifer which was based on the earlier hydrogeologic framework presented in Zapecza, (1989) for that aquifer. Available Piney Point aquifer withdrawal data indicate anomalously low rates in 2004, 2008, and 2010, likely resulting from missing data for wells with reported withdrawals in other years. Therefore, the 2007 withdrawal rates were substituted for 2004, 2008, and 2010 withdrawals for wells with presumed missing data. The average water-level residuals (simulated minus observed) are shown in table 14 and figures 37 and 38 for 26 wells in New Jersey. The average water-level residual for 26 wells in Cumberland, Atlantic, Gloucester, Camden, Burlington, and Ocean Counties in New Jersey is 9 ft and the RMSE is 23 ft. The minimum; first, second, and third quartile; and maximum residuals are –22, –2, 4, 16, and 80 ft, respectively. The simulated water levels are within 10 ft of observed water levels for 14 of the 26 wells (54 percent) in New Jersey and within 20 ft of observed levels in 19 of 26 wells (73 percent) (figures 37 and 38 and table 14). The average water-level residuals also are shown in table 14 and figure 37 for four wells in Delaware.

Hydrographs for two wells screened in the Piney Point aquifer are shown in figure 39. Simulated water levels in well 050676, located in an area in southeastern Burlington County where withdrawals are minimal, have the same flat trend but a smaller amplitude than observed water levels and are, on average, 16 ft greater than observed water levels (fig. 39A). The water levels simulated for this well by Voronin (2004) also show a smaller amplitude than observed water levels and are 5 to 10 ft lower than the observed water levels during 1980–98. The simulated water levels in well 110096 in Cumberland County have the same decreasing trend and amplitude as the observed water levels (associated with increased nearby withdrawals), but the simulated water levels are on average 9 ft greater than observed water levels (fig. 39B). The water levels simulated for this well by Voronin (2004) show a pattern similar to that of the simulated water levels shown in figure 39B for 1980–98.

Calibration of Water Levels in the Vincentown Aquifer

Simulated water level contours in the Vincentown aquifer (model layer 11, present primarily in Burlington, Gloucester, Monmouth, Ocean, and Salem Counties) generally match the 2013 potentiometric surface (Cauller and Gordon, 2021) (fig. 40). The differences between the simulated and observed water levels are shown in table 15 and in figures 40 and 41 for wells measured in New Jersey in 2013. The average water-level residual (simulated minus observed) for nine wells in Burlington, Gloucester, Monmouth, Ocean, and Salem Counties is −2 ft and the RMSE is 11 ft. The minimum; first, second, and third quartile; and maximum residuals are −9, −6, −2, 8, and 29 ft, respectively. The simulated water levels are within 10 ft of observed water levels for seven of the nine wells (89 percent) in New Jersey. The simulated water levels for the Rancocas aquifer in Delaware and the average water-level residual for one well measured in Delaware in 2013 (table 15) also are shown in figure 40.

Hydrographs of simulated and observed water levels in wells 250636 and 290139 in Monmouth and Ocean Counties, respectively, show that simulated water levels are an average of 2 ft and about 3 ft lower than observed water levels, respectively, and have the same flat trend and similar amplitude (fig. 42). The hydrograph for well 290139 in Voronin (2004) is similar to the hydrograph shown in figure 42B for 1980–98, whereas the hydrograph for well 250636 (fig. 42A) in Voronin (2004) shows simulated water levels that are 20 to 25 ft lower than observed water levels during 1980–98.

Calibration of Water Levels in the Wenonah-Mount Laurel Aquifer

Simulated water-level contours in the Wenonah-Mount Laurel aquifer (model layer 13, Monmouth, Ocean, Burlington, Camden, Gloucester, and Salem Counties) generally match the 2013 potentiometric surface (Cauller and Gordon, 2021), but the simulated water levels are lower in updip areas of Monmouth County and central Gloucester County (fig. 43). The simulated water-level contours generally are similar in shape to the 2013 potentiometric surface over much of the extent of the aquifer and the simulated water levels are within 10 or 20 ft of the water levels that were measured at the cones of depression associated with pumping centers in coastal Ocean and Monmouth Counties and western Burlington and central Camden Counties. The differences between the simulated and observed water levels for wells measured in New Jersey in 2013 are shown in table 16 and in figures 43 and 44. The average water-level residual (simulated minus observed) for 48 wells in Salem, Gloucester, Camden, Burlington, Ocean, and Monmouth Counties is −8 ft and the RMSE is 23 ft. The minimum; first, second, and third quartiles; and maximum residuals are –57, –23, –4, 7, and 55 ft, respectively. The simulated water levels are within 10 ft of observed water levels for 19 of the 48 wells (40 percent) and within 20 ft of observed water levels for 33 of the 48 wells (69 percent) in New Jersey. The simulated water levels for the Mount Laurel aquifer in Delaware and the average water-level residual for two wells measured in Delaware in 2013 (table 16) also are shown in figure 43.