15th ASCE Engineering Mechanics Conference
June 2-5, 2002, Columbia University, New York, NY
A THREE-DIMENSIONAL FINITE-ELEMENT MODEL OF SEAWATER
INTRUSION IN WESTERN LONG ISLAND, NEW YORK
Paul E. Misut1, Jack Monti Jr.2,
and Clifford I. Voss3
A three-dimensional version of the U.S. Geological Survey’s SUTRA (Saturated-Unsaturated TRAnsport)
ground-water-flow model was tested for its ability to simulate seawater intrusion into the ground-water
system of western Long Island, N.Y. This model solves two equations–a fluid-mass balance for
unsaturated and saturated ground-water flow, and a solute mass balance. The solute concentration
corresponds to the mass fraction of total dissolved solids and accounts for density dynamics driven by
solute concentration. An extensive geographic information system of the region is the basis for the
aquifer-geometry and boundary-conditions representation. The direct-banded matrix solver that is
typically used with two-dimensional SUTRA models has large computer-memory requirements and was
replaced with an iterative solver. Model simulations indicate that the hydrologic system has not reached a
steady-state configuration in response to pumping, changes in recharge, and postglacial sea-level rise.
U.S. Geological Survey computer programs have been widely used to simulate problems that
entail seawater intrusion in ground-water systems. (Many of these codes with supporting
documentation are available free of charge at http://water.usgs.gov/software.) The SUTRA code
(Voss, 1984) can simulate variable fluid density, a critical feature in seawater intrusion problems.
The majority of SUTRA applications developed to date use a cross-sectional approach; however,
use of cross-sectional models to predict the possibility of seawater intrusion in response to
pumping on western Long Island is not advisable because the system has three-dimensional
flow patterns (Kontis, 1999). This paper describes a test of a new three-dimensional version of
the SUTRA code as applied to seawater intrusion in western Long Island.
U.S. Geological Survey, 2045 Rt. 112 , Coram, NY 11727. E-mail: email@example.com
U.S. Geological Survey, 2045 Rt. 112 , Coram, NY 11727. E-mail: firstname.lastname@example.org
U.S. Geological Survey, 12201 Sunrise Valley Dr. , Reston, VA 20192. E-mail: email@example.com
Development of the model required: (1) delineation of the extent, thickness, and hydraulic
characteristics of aquifers and confining units, and (2) definition of hydrologic-boundary
conditions, among which are recharge (from precipitation and lateral inflow of ground water),
discharge (to streams, the shore, subsea parts of the model, and wells), and no-flow boundaries
The sequence of aquifers and confining units in the modeled area consists of four aquifers
and two confining units. These have been described in detail by Misut and Monti (1999) and
many others; they are, from the bedrock (the bottom of the model) upward, the Lloyd aquifer, the
Raritan clay, the Magothy aquifer, the Jameco aquifer, the Gardiners Clay, and the upper glacial
aquifer. These six hydrogeologic units vary in thickness and, in some areas, pinch out. In
pinch-out areas, the model mesh was continued with a minimal thickness, and hydraulic
properties were extrapolated from the overlying layer. Representation of the spatial distribution
of thicknesses and hydraulic properties was facilitated through use of digital maps of
land-surface elevation, bathymetry, and the bottom elevations of geologic units. The
Environmental Systems Research Institute (ESRI, http://www.esri.com) GRID format was used to
store these data at a resolution of 30 meters to match the available LANDSAT satellite imagery.
(See http://edcwww.cr.usgs.gov/webglis/ for information on LANDSAT.) An interpolation of
these GRID data onto finite-element mesh nodes in ESRI shapefile format was then imported
into a graphical user interface to SUTRA (Voss, Boldt, and Shapiro 1998) as modified for
three-dimensional simulation by Winston and Voss (written communication, 2001). This
graphical user interface and the mesh-to-shapefile converter are plug-in extensions to the
ArgusONE software (http://www.argusint.com).
The top of the three-dimensional mesh and its boundary conditions are depicted in figures
1A and 1B, respectively. Topographic data were used to delineate the specified-flux boundary
where ground-water recharge is received from precipitation. Hydrostatic-pressure boundaries
with seawater concentration in offshore zones were calculated from bathymetric data. Streams
and ponds where discharge occurs are represented as specified-pressure boundaries.
FIG.1. Block diagrams of western Long Island, New York showing: A.
Land-surface area and offshore model area, B. Top of mesh and its boundary
The vertical sides of the mesh and the stratigraphic representation are depicted in figure 2.
The thickness of the model increases dramatically from north to south. The mesh contains 30
node layers and a total of about 54,000 elements. Five node layers were used to represent each of
the four aquifers, and three node layers were used to represent the two confining units to allow
curvature in simulated saltwater interfaces within hydrogeologic units. Three node layers were
used to represent the unsaturated zone. Outwash and moraine zones were used to represent
permeability differences in the upper glacial aquifer and the unsaturated zone (total of 9 node
layers); permeability is uniform throughout each of the other node layers and varies only
depending on the hydrogeologic unit that is represented. Horizontal placement of nodes was
facilitated by maps showing the locations of observation and pumping wells, surface-water
features, deeper geometric features such as an ancestral river valley, and refined discretization for
numerical purposes near the saltwater-interface locations. Hydrostatic-pressure boundaries
corresponding to the seawater chloride concentration of 19,000 mg/l were used in offshore parts
of the model domain. Hydrostatic-pressure boundaries with freshwater chloride concentration
of 0 mg/l were used at some streams and ponds and in parts of the eastern side of the model to
represent a connection to eastern Long Island.
FIG. 2. Interior and lateral boundary conditions and the 30 node layers
representing aquifers and confining units.
The bottom of the mesh represents the top of the crystalline bedrock and is depicted from above
in figure 3. It is considered a no-flow boundary because its permeability is negligible.
Bedrock-surface elevation ranges from about sea level in the northwestern part of the model to
600 meters below sea level in the southeast corner.
FIG. 3. Bottom of mesh (bedrock surface, a no-flow boundary).
Calibration of simulated pressures and chloride concentrations to their observed values
entailed adjusting specified pressures (along the eastern side of the model), porosity, storativity,
permeability, pumpage, and dispersivity. Final values of the constants were as follows: porosity:
15 percent; storativity: 10 -8 (kg/(m ⋅ s 2 )) −1 ; longitudinal dispersivity in the direction of
maximum permeability, 2,000 m; longitudinal dispersivity in the direction of minimum
permeability, 20 m; and transverse dispersivity, 5 m.
Statistics on residuals of head and chloride concentration in relation to 1999 steady-state
conditions are given by geologic unit in table 1. Negative values indicate that simulated values
were greater than observed values. Heads were measured in March 2000; simulated pressure and
density were converted to head at the elevation of the observation-well screen. Chloride
concentrations were measured throughout the 1990s.
Table 1. Steady-state calibration residuals (head and chloride concentration) for aquifers and
confining units in simulations of freshwater/saltwater interface on western Long Island, N.Y.
[RMS, root mean squared]
Mesh Head residual Chloride concentration residual
Unit layers (feet above sea level) (milligrams per liter)
Mean RMS Population Mean RMS Population
Moraine aquifer 4-9 -2.3 6.7 28 -1798 2857 5
Outwash aquifer 4-9 -2.5 8.3 28 2458 6112 11
Gardiners Clay 10 - 12 -0.6 5.3 7 543 3112 17
Jameco Gravel 13 - 17 -4.9 6.3 5 658 1652 8
Magothy aquifer 18 - 22 -3.2 5.8 11 189 2838 106
Raritan Clay 23 - 25 1.5 5.1 7 -1838 6185 23
Lloyd aquifer 26 - 30 -0.5 6.7 8 -7790 11,556 34
TOTAL 1 - 30 -2.0 6.9 97 -1247 5821 204
Negative chloride residuals in the Lloyd aquifer and Raritan Clay are due to a deviation of the
simulated steady-state configuration from conditions measured in the field. At present(2002),
seawater intrusion is occurring as a result of post-Pleistocene sea-level rise. The
freshwater/saltwater interface in deep parts of the northern Atlantic Coastal Plain has not yet
reached a steady-state position (Meisler and others, 1984); therefore, exploratory simulations
were conducted to investigate the movement of the interface over the last 10,000 years. No
definitive location of the interface could be obtained, however, because offshore data were
insufficient; therefore, the residuals shown in table 1 reflect simulated steady-state conditions,
not actual conditions. The simulated chloride concentrations provided by the best steady-state
calibration adjustments are too high in the Lloyd aquifer and Raritan clay because seawater
continues to intrude into these zones.
An important geometric feature of the model is three holes within the Gardiner’s Clay,
each several miles in diameter near the south shore of Long Island. The Gardiners Clay is the
uppermost confining unit (table 1), and these holes allow flow from the upper glacial aquifer into
the Jameco and Magothy aquifers below. The exploratory simulations that were conducted to
investigate the movement of the freshwater/saltwater interface over the last 10,000 years
indicated that seawater can move downward into parts of the freshwater flow system through
these holes, possibly hastening the regional movement of the interface from the south, where
conditions are confined. The sea above these holes is fairly shallow; therefore, this effect depends
in part on the rate of sea-level rise.
Two surfaces of equal simulated pressure are depicted in figure 4. The top surface, which
mimics the land and coastline, represents the water table, where pressure is equal to zero. Parts
of the model domain are above this surface–for example, in the northeastern corner. The flat,
triangular bottom surface is where pressure is equal to 3 x 10 6 kg/(m ⋅ s 2 ) . Simulated pressure at
the southeastern bottom corner of the model domain, where the bedrock surface is deepest, is
about 5 x 10 6 kg/(m ⋅ s 2 ) .
FIG. 4. Two surfaces of uniform simulated pressure, western Long Island, New
The surface of equal simulated chloride concentration (2,000 mg/l) is depicted in figure 5.
This seawater interface is near the coast and extends landward at depth in response to natural
gradients and pumping. Three breaks interrupt the continuity of the interface. The
northwestern break (B1 in fig. 5) is where the bedrock surface extends above sea level. The
northeastern break (B2 in fig. 5) is caused by the extension of the freshwater flow system beyond
the model domain into a north-shore peninsula (not modeled). The eastern break (B3 in fig. 5)
is caused by the extension of the freshwater flow system beyond the model domain into eastern
Long Island. Localized curvature of the interface surface is caused by simulated pumping and
by the model geometry. For example, a shelf is present beneath Jamaica Bay, where the Gardiners
Clay confines freshwater below saltwater. The interface becomes more vertical east and west of
the Jamaica shelf as a result of the three holes in the Gardiners Clay.
An analysis was done to calculate the amount of time necessary to reach steady-state
conditions from an initial state in which the interfaces are offshore, near their inferred present
position (Misut and Monti, 1999). Development of this initial state, in which the pressure field
is consistent with the concentration field and with boundary conditions, entailed simulation of the
interfaces of the Raritan Clay and Lloyd aquifer with a second preliminary, transient simulation.
This preliminary simulation had no pumping stresses and started the deep interfaces near the
model boundaries, then continued forward for 6,000 years, stopping at the inferred 1900’s
predevelopment location. The simulation from the 1900’s predevelopment location to a
predevelopment steady-state (with no pumping or changes in recharge) requires an additional
several thousand years of hypothetical time. This hypothetical time was sensitive to the amount
of cross-flow that was allowed through the holes in the Gardiners Clay.
FIG. 5. Surface of equal simulated chloride concentration (2000 mg/l), and location
of breaks (B1-B3) in surface continuity, western Long Island, New York.
The effect of pumping and changes in recharge on the movement of the interfaces at depth
was also analyzed. Incorporating these stresses accelerated the landward movement of the
interface and caused seawater intrusion beyond the predicted prevelopment steady-state
configuration in some areas. In other areas, however, simulation of pumping and recharge
starting with the interface at the 1900’s predevelopment location did not result in a close match
with measured increases in chloride concentration since 1920. The position of the interfaces at
depth is highly uncertain at present and may be farther inland and, thus, closer to steady state
than inferred previously by Misut and Monti (1999).
The solution of the western Long Island model requires large computational resources.
Although the new three-dimensional version of the SUTRA code used in this test leverages
modern PC computer technologies such as double-data-rate RAM, a gigahertz-speed CPU, and a
10,000-rpm hard drive, SUTRA’s traditional numerical methods were nevertheless found to run
slowly for a 30-layer mesh with about 54,000 finite elements–about 3 hours per time step. A
faster run of about 5 minutes per time step was attained by replacing the direct-banded matrix
solver with the Sparse Linear Algebra Package (SLAP) (Seager, 1988) of routines for
approximately solving large, sparse systems of linear equations. SLAP does not store the entire
matrix; preconditioned iterative methods store only the nonzero elements and their row and
column numbers. The extra effort at optimizing SLAP parameters and the errors introduced by an
approximation method were outweighed by the benefits of faster model-run times.
A three-dimensional version of the U.S. Geological Survey’s SUTRA ground-water-flow model
code was tested for its ability to simulate seawater intrusion into the ground-water system of
western Long Island, N.Y. The SUTRA code seems to accurately represent ground-water flow
and solute transport in the area, including the effect of pumping and recharge on seawater
intrusion. The amount of simulation time required for the interfaces to reach their present
position depends on factors that remain uncertain, namely–the rate of sea-level rise since the
Pleistocene glaciation and the amount of flow that moves through offshore holes in the Gardiners
Clay. Given reasonable assumptions about these and related factors, the model can reliably
predict the movement of the interfaces in response to anticipated future stresses.
The authors thank the New York City Department of Environmental Protection
(http://www.nyc.gov/html/dep/) and Argus Interware LLC for their cooperation in this study.
Any use of trade, product or firm names is for descriptive purposes only and does not imply
endorsement by the U.S. Government.
______2002, Ground-Water Software, accessed February 26, 2002 at
______2002, Environmental Systems Research Institute, accessed February 26, 2002 at
______2002, LANDSAT, accessed February 26, 2002 at http://edcwww.cr.usgs.gov/webglis/.
______2002, Argus Interware, accessed February 26, 2002 at http://www.argusint.com.
______2002, New York City Department of Environmental Protection, accessed February 26,
2002 at http://www.nyc.gov/html/dep/.
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.
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.
Seager, M. (1988) A SLAP for the Masses: Lawrence Livermore Nat. Laboratory Technical
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.