A SUTRA Model of Seawater Intrusion
in Western Long Island, New York

 P.E. Misut 1, W. Yulinsky 2, D. Cohen 3, D. St. Germain 3, C. I. Voss 1, J. Monti Jr. 1

1 U. S. Geological Survey, USA

2 N. Y. C. Department of Environmental Protection, USA

3 Malcolm Pirnie Inc., USA

 

INTRODUCTION:  New York City obtains about 99 percent of its freshwater from a surface-water system. The supply could run short, however, during a prolonged drought or as a result of a system malfunction. The aquifers of eastern New York City (Brooklyn and Queens, fig. 1) could provide a substantial supplement for use in times of shortage, but excessive pumping of these aquifers in the past has caused seawater intrusion, resulting in the shutdown of muncipal supply wells throughout Brooklyn and in several parts of Queens.

The coastal-plain aquifers are recharged by abundant precipitation and currently provide about 1 percent of the total supply for New York City. This water is provided mainly to local communities that lack connection to the city’s surface-water system. Rehabilitation of the existing infrastructure could provide a supplemental supply of about 2,000 kg/s, or 3 percent of the city’s total supply, if the surface-water supplies are diminished by drought or system failures.

In 2001, the U.S. Geological Survey, in cooperation with New York City Department of Environmental Protection (NYCDEP; http://www.nyc.gov/dep), began a 3-year study to develop a model to simulate a wide range of pumping scenarios and predict the location and degree of seawater intrusion that would occur if ground water were used to make up a deficit in the city’s supply. Hypothetical pumping scenarios that could increase the ground-water yield to more than 10,000 kg/s are being evaluated (Malcolm Pirnie, 1999). This abstract summarizes the coastal-plain aquifer system, describes the model and the initial conditions, and presents some results of interface-movement scenarios.

HYDROLOGY: The ground-water system on western Long Island (Brooklyn and Queens) consists of four unconsolidated aquifers and two confining units and is underlain by gently southward sloping impermeable bedrock (fig. 2). The bedrock crops out in the northwest, is about 300 m below sea level in the southeast along the south shore, and about 250 m below sea level in the southwest along the south shore. The sequence of unconsolidated deposits, from the bedrock upward, consists of the Lloyd aquifer, the Raritan clay, and the Magothy aquifer, all of Cretaceous age, and the Jameco aquifer, the Gardiners Clay, and the upper glacial aquifer (with outwash and moraine zones), all of Pleistocene age. These six units vary in thickness locally and pinch out in some areas (Smolensky and others, 1989). At several locations, the Gardiners Clay contains erosional holes that provide vertical hydraulic connections. The water table is mostly in the upper glacial aquifer.

The distribution of freshwater and the position of the freshwater/seawater interface do not conform to the Ghyben-Herzberg relation as a result of local anisotropy and varying unit thickness. The response of the seawater interface to the sea-level rise since the Pleistocene has been delayed by hundreds to thousands of years in the Lloyd aquifer beneath the Atlantic Ocean as a result of confinement by the overlying Raritan Clay (Meisler and others, 1984). Pumping during the mid-20th century caused increased salinity in water pumped from wells screened in the upper glacial, Jameco, and Magothy aquifers and necessitated the shutdown of all public-supply wells in Brooklyn in 1947 (Soren, 1976).

Inflows to and outflows from the aquifers are currently near steady state. About half of the annual precipitation (1.1 m/yr) recharged the flow system before urbanization, but the present extensive stormwater- and sanitary-sewer systems probably decrease recharge to one-third of the annual precipitation (Misut and Monti, 1999) About 440 kg/s of fresh ground water enters the aquifers system laterally from Nassau County to the east. Natural outflow occurs mostly near the shore, but partly as freshwater seepage to ponds and streams and as subsea discharge. The water supply that is piped from upstate surface-water reservoirs discharges to the sanitary-sewer system, bypassing the ground-water system.

MODELING. Numerical modeling of the Long Island ground-water-flow system began in the 1970's (Getzen, 1974, Getzen, 1977, Gupta and Pinder, 1978). Three-dimensional transient-state MODFLOW models were developed in the 1980's and 1990's (Reilly and Harbaugh, 1980; Buxton and Smolensky, 1999; Misut and Monti, 1999), and cross-sectional SHARP models were developed in the 1990's (Kontis, 1999). Development of a three-dimensional SUTRA model began in 2001. The MODFLOW, SHARP, and SUTRA codes with supporting documentation are available free of charge at http://water.usgs.gov/software. The New York City MODFLOW models used fixed no-horizontal-flow boundaries in combination with vertical leakage across confining units into constant-head boundaries to represent seawater interfaces. The SHARP models provided only a limited physical mechanism to represent movement of seawater interfaces in response to pumping.

The SUTRA model code (Voss, 1984) and graphical user interface (Voss and others, 1997) were chosen for this study for its ability to simulate variable fluid density driven by solute mass and/or temperature on three-dimensionally deformed finite-element meshes. A raster format was used to store hydrogeologic and hydrologic data at a resolution of 30 m to match available satellite imagery. Raster-formatted data were interpolated onto points corresponding to finite-element meshes; these point-wise data were then imported into a graphical user interface for SUTRA.

Boundary Conditions:  Specified flux was used to represent the ground-water recharge mechanism. Five distributions of recharge were prepared to correspond to rates under the following conditions:  predevelopment, two simplified historical stages of urbanization, current (2000) conditions, and 1962-66 drought conditions. The predevelopment period was assigned a uniform recharge distribution of 0.56 m/yr in inland areas. The simplified historical stages were represented by two zones that change historically in size and shape—an undeveloped zone, which was assigned a recharge rate of 0.56 m/yr, and an urbanizing zone, in which recharge was decreased to 0.25 m/yr to account for consumptive water use and conveyance of storm-water runoff to streams. The current-conditions distribution of recharge  was based on satellite imagery. The drought-simulation period entailed a scalar reduction in the current-conditions distribution.

Hydrostatic-pressure boundaries with saltwater chloride concentration of 19,000 mg/l were used in offshore zones. Pressures at sea-floor nodes were calculated from NOAA bathymetric data (http://www.noaa.gov). Hydrostatic-pressure boundaries with freshwater chloride concentration of 0 mg/l were applied at some streams and ponds and to the entire eastern side of the model to represent the continuation of aquifers beyond the model domain. A no-flow boundary was used to represent the model bottom (bedrock).

Model Discretization: A graphical user interface facilitated rediscretization through a one-step interpolation process that allows the model results to be used as initial conditions for a later run with a different mesh (fig. 3). Meshes were designed with varying grid density to represent:
(1) inflows across the eastern (Nassau County) border;
(2) movement of the south-shore seawater interfaces in response to post-glacial sea-level rise;
(3) local conditions at NYCDEP wellfields; and
(4) the two holes in the Gardiners Clay.
The dimensions of all of these meshes were designed to meet the requirements of being small enough to be held in RAM, yet large enough to minimize element sizes and provide the greatest accuracy. The meshes each consisted of about 70 layers of 2,000 elements, and each element represented an average area of about 2.5 km
2 . Ten layers were needed to represent convection flow patterns in each hydrogeologic unit (fig. 4). Time discretization of model runs also varied, depending on the purpose; the maximum time step in the simulation of sea-level rise did not exceed 10 yr, and simulations of manmade stresses (pumping) generally used 0.1-yr time steps.

Initial Conditions: The position of the seawater interface in the Lloyd aquifer and Raritan Clay still reflects the lowered sea level of the Wisconsin glaciation, about 10,000 years ago. Exploratory SUTRA simulations were conducted to investigate the post-glacial interface movement and to generate initial conditions for pumpage-evaluation runs that start at the present and allow the interfaces to move slowly landward. Runs that started with a model domain saturated completely with freshwater and applied time-invariant boundary conditions corresponding to current sea level did not approach inferred current interface positions within a 10,000-yr post-glacial period; they took much longer. Therefore, simulation of pre-glacial conditions may also be required. Furthermore, the inferred current interface configuration was never reached during these simulations. One explanation for this is the presence of holes in the Gardiners Clay, which are slightly below the present sea level and thus become exposed to seawater at some undetermined time leading up to the present as sea level rises. The arrival of the rising sea at a Gardiners Clay hole opens a new path of least resistance for density-driven downward flow of seawater. The rate of sea-level rise is uncertain, however, and this process cannot be represented by time-invariant boundary conditions in which sea level is fixed. Finally, an analysis was done to calculate rates of interface movement in a pre-development simulation, with current sea level, in which the interfaces were started at their inferred current position. Only the interfaces in the Lloyd and Raritan moved significantly and it can be assumed that the interfaces of the Magothy, Jameco, and upper glacial aquifers are at steady state at the current time. The urbanization period (1900-2000) was then simulated and provided results that were satisfactory as initial conditions for the pumpage-evaluation runs.

Model Calibration. Parameters were adjusted during a sequence of about 50 model calibration runs to attain the best match of simulated pressures and concentrations to field data.  Final values of the most sensitive parameters are as follows:

Porosity: 15 percent

Storativity: 10 -8 {kg/(m s2)}-1;

Permeability:
outwash horizontal: 1 x 10
-10 m/s,                         outwash vertical: 1 x 10 -11 m/s,
moraine horizontal: 2 x 10
-11 m/s                           moraine vertical: 2 x 10 -12 m/s
Gardiners Clay: 2 x 10
-15 m/s (isotropic),
Jameco horizontal: 2 x 10
-9 m/s,                            Jameco vertical: 2 x 10 -10 m/s,
Magothy horizontal: 6.5 x 10
-11 m/s                      Magothy vertical: 6.5 x 10 -13 m/s,
Raritan: 2 x 10
-15 m/s (isotropic),
Lloyd horizontal: 2.7 x 10
-11 m/s,                          Lloyd vertical: 2.7 x10 -12 m/s.

Except for the upper glacial aquifer, hydraulic properties of hydrogeologic units are uniform, and entail extrapolation to the offshore part of the model domain, where little data was available to calibrate against.

SIMULATION RESULTS. The drought-emergency reactivation scenario entailed rehabilitating wellfields that were taken out of service and replaced by surface-water supply. If a drought were to decrease the surface-water supply, these wells would be reactivated to produce a total of about 2,000 kg/s, about 80 percent greater than the amount currently pumped. Each of these wells is screened in either the Magothy, the Jameco, or the upper glacial aquifer. The location of the seawater interface in the Lloyd aquifer is not critical for this particular simulation, but other scenarios are being considered that involve new wells in the Lloyd aquifer. Areas with a potential for seawater intrusion were identified through an analysis of contours of change in solute-mass fraction from current conditions to the simulation result after 2 yr of pumping (fig 5 and fig 6). The pumping period was limited to 2 yr because longer periods resulted in landward water table gradients from the coast to the pumping wells.  A change equivalent to a 50 mg/l or greater increase in chloride concentration was identified as an indicator of significant intrusion potential; areas with this potential were in the nearshore areas closest to pumping wells and closest to the hole of the Gardiners Clay. Changes smaller than 50 mg/l may be attributable to numerical inaccuracies and/or naturally occurring changes resulting from post-Pleistocene sea-level rise and, therefore, warrant site-specific investigation. The shortest distance between an area with significant intrusion potential and a pumping well exceeds 3 km; therefore, the risk of well contamination within 2 yr of supplemental pumping is small.

 

Keywords: SUTRA, numerical modeling, finite element method

Corresponding author: Paul Misut, Hydrologist, USGS, 2045 Rt. 112, Coram, NY 11727. Email: pemisut@usgs.gov

References

______2002, Ground-Water Software, accessed February 3, 2003 at http://water.usgs.gov/software.

______2002, New York City Department of Environmental Protection, accessed February 3, 2003 at http://www.nyc.gov/dep.

Buxton, H.T., Smolensky, D.A. (1999) Simulation analysis of the development of the ground-water flow system of Long Island, New York: U.S. Geological Survey Water-Resources  Investigations Report 98-4069, 57 p.

Getzen, R.T. (1974) The Long Island ground-water reservoir –a case study in anisotropic flow: Urbana-Champaign, Illinois, University of Illinois, Ph.D. thesis, 130p.

Getzen, R.T. (1977) Analog-model analysis of regional three-dimensional flow in the ground-water reservoir of Long Island, New York: U.S. Geological Survey Professional Paper 982, 49 p.

Gupta, S.K., and Pinder, G.E. (1978) Three-dimensional finite-element model for multilayered ground-water reservoir of Long Island, New York: Princeton University, Department of Civil Engineering, Research Report 78-WR-14 (unpaginated).

Kontis, A.L. (1999) Simulation of Freshwater-Saltwater Interfaces in the Brooklyn-Queens Aquifer System, Long Island, New York: U.S. Geological Survey Water-Resources Investigations Report 98-4067, 26 p.

Malcolm Pirnie (1999) Feasibility Study for use of the Brooklyn Queens Aquifer as an Additional Potable Water Supply: Main Report (unpaginated).

Meisler, H., Leahy, P.P, and Knobel, L.L. (1984) Effect of Sea-Level Changes on Saltwater-Freshwater Relations in the Northern Atlantic Coastal Plain: U.S.Geological Survey Water-Supply Paper 2255, 28 p.

Misut, P.E., and Monti, J., Jr. (1999) Simulation of Ground-water Flow and Pumpage in Kings and Queens Counties, Long Island, New York: U.S. Geological Survey Water-Resources Investigations Report 98-4071, 50 p.

Reilly, T.E., and Harbaugh, A.W. (1980) A comparison of analog and digital modeling techniques for simulating three-dimensional ground-water flow on Long Island, New York: U.S. Geological Survey Water-Resources Investigations Report 80-14, 40 p.

Voss, C.I., Boldt, D.R., and Shapiro, A.M. (1997) A graphical-user interface for the U.S. Geological Survey's SUTRA code using Argus ONE (for simulation of variable-density saturated-unsaturated ground-water flow with solute or energy transport): U.S. Geological Survey Open-File Report 97-421, 106 p.

Voss, C.I., (1984) A finite-element simulation model for saturated-unsaturated, fluid-density-dependent ground-water flow with energy transport or chemically-reactive single-species solute transport: U.S. Geological Survey Water-Resources Investigations Report 84-4369, 409 p.