WATER RESOURCES RESEARCH, VOL. 23, NO.

10, PAGES 1851-1866, OCTOBER 1987

Variable Density Flow and Solute Transport Simulation

of Regional Aquifers Containing a Narrow
Freshwater-Saltwater Transition Zone

CLIFFORD I. VOSS

U.S. Geological Survey, Reston, Virginia

WILLIAM R. SOUZA

U.S. Geological Survey, Honolulu, Hawaii

Variable density flow and solute transport simulation of aquifer systemscontaining narrow transition
zones between freshwater and saltwater requires particular attention to certain aspectsof the numerical
method and its application to be successful.Typically, only casesinvolving wide transition zones have
been simulated with variable density transport models, possibly because of inaccuraciesin the modeling
approaches used. The major components of a successfulapproach are threefold. First, functionally
consistent approximation of terms involved in fluid velocity calculations is necessary.In the case of
Galerkin finite element methodology, a significantmodification is required to the standard approach in
order to achieve consistency.Second, the simulator must be verified in a particular series of tests. The
usual tests using Henry's problem for verification of density-dependent transport simulators are inad-
equate to check for consistencyof the velocity approximations and for the accuracy of simulating flow
driven by bouyancy forces.Third, adequately fine spatial discretizationis required when applying the
simulatorfor spatial stability of the numericaltransportsolutionand to allow accuraterepresentationof
narrow transition zones and the effectsof low transversedispersivity.The effectivenessof this approach is
demonstrated through simulation of the flow of fresh and saline groundwater in the layered basalt
aquiferof southernOahu, Hawaii. The transitionzone in this regionalflow systemis narrow exceptnear
the dischargearea where it is broadly dispersed.Simulation of this common situation with an inconsis-
tent approximation gives grosslyincorrect results,while simulation with a consistentmodel provides a
robust tool for analysisof systemhydrology.

INTRODUCTION freshwater-saltwater systems containing portions of narrow

The problem in analyzing a regional aquifer systemcontain- transition zones ar• skirted by basing the analysison a variety
ing a narrow transitionzone betweenfreshwaterand saltwater of sharp interface approximations. The various areal and
with a variable density fluid flow and dissolved solids trans- cross-sectionalapproachesare reviewed by Reilly and Good-
port model is twofold. man [1985]. While not minimizing the great value of such
1. The adequacy of the Fickian dispersiontheory is not approaches in representingfundamental aspects of aquifer
certain, and values of dispersionparameters are neither well- systemdynamics,they do not account for the often significant
known nor may they be directly measured at the regional effect of the mixing processon the structure and position of
scale. The physical basis of a particular concentration transi- the transition zone, nor do they allow a means of direct pre-
tion zone structure may be studied only after careful fitting of diction of the dissolved solids concentration of water arriving
simulated hydraulics to the field conditions.This requires that at supply wells. The impact of theselimitations is clear when
the transport simulationaccuratelyrepresentthe physicalgov- consideringthat water of potable quality contains less than
erning equationsbeforeevaluation of the adequacyof old and about one percentseawater.On the other hand, regional scale
new dispersion models and the determination of parameter analysisbased on variable density flow and transport simula-
values at the field scale can be undertaken. tion is capable of representingsystemdynamicsincluding the
2. Transport codes,however,have considerabledifficulties effectsof mixing, as founded on current theories of dispersion
simulating the movement of narrow concentrationfronts on a in aquifers, and is capable of making predictions of con-
centration.
regional scale.The difficultiesare (1) spatial instabilitiesin the
numerical concentration solution when tracking a narrow Often in regional aquifers, the transition zone is relatively
front, (2) inability to maintain the narrownessof a front due to narrow except in areas that undergo pumping, recharge, dis-
insufficient discretization, (3) significant numerical errors in charge,or tidal stresses. The ability to simultaneouslysimulate
calculating fluid velocity within narrow transition zones caus- the dynamics of both narrow and broad portions of the transi-
ing broadening of the zones, and (4) inaccuraciesin repre- tion zone is vital to hydrologic analysis and water supply
sentation of bulk fluid movement driven by density differences. prediction in suchcases.However, transport simulationshave
Classically,the problemsinherent in transport simulation of almost exclusively dealt with transition zones that are quite
broad on the regional scale(on the order of aquifer thickness
or greater).These include analysesreported by Lee and Cheng
This paper is not subjectto U.S. copyright. Publishedin 1987 by [1974], Segol and Pinder [1976], Volker and Rushton [1982],
the American Geophysical Union.
Frind [1982b], and Lebbe [1983]. Some authors [Volker,
Paper number 6W4664. 1980; Frind, 1982a] have shown results for narrower transi-
1851
 1852 Voss AND SOUZA:SIMULATIONOF NARROWTRANSITIONZONES

tion zones. Speculation may be warranted, however, on the in the aquifer,assumingthe contributionof solutedispersion
question of whether broad transition zones are often chosen as to the massaverageflux of fluid is negligible,may be givenby
candidates for transport simulation analysis becauseof the [Bear, 1979]
numerical problemsin representingnarrow zones,and further,
8p 8p 8C
whether in somecases,inaccuraciesin variable densitytrans-
port modelsgive broad transitionresultsthat are acceptedfor
pSøP
'• -•-E(•'•(•W
-•-V. (Epv)
= Qv (4)
lack of definitivefield data showingotherwise. into which Darcy's law (3) for mass average velocity may be
This paper describesa numerical modeling approach that
substituted.
The specificpressure storativitySopfor a rigid
dealswith the above-mentionedproblemsin representationof solidaquifermatrix(M/(L r2)) - 1 is givenby
freshwaterand saltwater systems.The approch guarantees
adequateregionalscaletransport simulationof narrow transi- Sop= + e0 ' (5)
tion zoneswith minimalnumericaldispersion.Throughappli- where• is porousmatrix compressibility (M/(LT2)) - •, and fi
cation of the modelingtechniqueto the southernOahu aqui-
is fluid compressibility
(M/(LT2)) - •
fer, Oahu, Hawaii, it is demonstratedthat transport simula-
A general boundary condition for the fluid mass balance
tion can provide an excellenttool for hydraulicinvestigation
that appliesat stationaryboundariesas well as approximately
even when narrow transition zones are involved.
at a water table (seeBear [1979, pp. 98-99] and Appendix B
for developmentof similar conditions)is
THEORY

The solute mass balance per unit aquifer volume at a point Syc•p
m • + epv. n + QIN* = 0 (6)
in an aquifer crosssection with variable density fluid is given
by [Bear, 1979] where
3C
Sy(x,z) watertablespecificyield(volumefluid released/
ep•- + epv
. VC- V. [eP(Oml
+ D)-VC]
= Or(C*
- C) (1) aquifer volume) for unit drop in hydraulic head, 1;
where Igl magnitudeof gravitationalacceleration,
L/T2;
Qis* fluid mass source due to flow acrossboundaries
C(x, z, t) solute concentration as a mass fraction, (massfluid rechargedper unit area of boundary/time)
(masssolute/massfluid) in units Ms/M, M/L2T.
where M s are units of solutemassand M are n unit outward normal to the boundary, 1.
units of fluid mass;
e(x, z) aquifer volumetric porosity, 1; This condition is accuratewhen at a water table, IVpl • plgl,
v(x, z, t) fluid velocity, in units L/T, where L are length that is, when pressuregradientsare nearly hydrostaticat the
units and T are time units; water table. This condition is usually satisfied becausecom-
Dm the moleculardiffusioncoefficientof solutein monly vertical pressuregradients are nearly hydrostatic and
pure fluid including aquifer material horizontal gradientsare much smallerthan vertical gradients.
tortuosityeffects,L2/T; The mechanicaldispersiontensor for an isotropic porous
I identity tensor, 1; medium in two spatial dimensionsis given by
D(x, z, t) mechanicaldispersiontensor,L2/T;
D= -DxxDxz]
Qv(x,z, t) fluidmasssource(massfluid/aquifer
volume/time)M/L 3 T;
_D:xD:•] (7)
C*(x, z, t) concentration of solute as a mass fraction in the where

sourcefluid, Ms/M;
p(x,z, t) fluiddensity,
M/Lf 3,where
Lf 3 isfluidvolume. (8)
Density is given as a linear function of concentration:

= +
8p
(c- Co) (2)
D•=(l-•2•)(drv•,
2+dtY•
2) (9)

1
where P0 is fluid density when C- Co; CO is a base solute
concentration; and c•p/SCis a constant coefficientof density
Oxz
=Ozx=(['•)(dl_.
- dr)(VxV:) (lO)

variability. dœ= %lvl (11)

Darcy's law gives the mass average fluid velocity at any
dT= •rlv[ (12)
point in a crosssectionas
whereIvlis the magnitudeof velocityand dL and dT are longi-
v= - ß(Vp - pg) (3)
tudinal and transversecoefficientsof mechanical dispersion
where L2/T,and0•L(X,
Z)is longitudinal
dispersivity,
L, and0•T(X,
Z)is
transversedispersivity,L.
k(x, z) permeabilitytensor,/?;
# fluid dynamic viscosity,M/LT; NUMERICAL METHOD WITH CONSISTENT VELOCITIES
g gravityvector,L/T2;
The first component of a successfulapproach to simulation
p(x, z, t) fluid pressure,
M/L T 2.
of variable density fluid movement is a consistent velocity
The massbalanceof fluid per unit aquifer volume at a point approximation. While the general conclusionsof the following
 MOSSAND SOUZA: SIMULATIONOF NARROW TRANSITIONZONES 1853

discussion on this topic are applicable to any numerical discontinuous from element to element [e.g., Segol et al.,
method, the ideas are developed on the basis of the finite 1975]. This inconsistencyalso would cause problems when it
element method. A weighted residual numerical method com- exists in simulators based on other methods and on other
bining Galerkin finite elements with integrated finite differ- hydraulic variables such as equivalent freshwater head, or a
encesdescribedin Appendix A is used to solve the two gov- quasi-streamfunction.
erning equations, (1) and (4). The Galerkin finite element The inconsistencyis most clearly demonstratedby consider-
method allows great geometric flexibility in mesh design and ation of hydrostaticconditionsfor a region of spacein which
gives robust direction- and anisotropy-independentrepre- pressure and concentration are approximated with a linear
sentation of fluid and solute fluxes.An integrated finite differ- variation from the top to bottom of the region. If the pressure
ence approximation for the spatial integration of all nonflux at thetopofa region withheigh• H iszero,thenthebottom
terms in the governing equations provides an economical pressure, under hydrostatic conditions, is given by PB=
alternative to the Galerkin method while giving accuracysuf- PAVGIgIH,where PAVGis the average density of fluid in the
ficient for any mildly nonlinearsimulationproblem. A back- region. If linear relation (2) holds between density and con-
ward finite differenceapproximation discretizesthe time de- centration, then when the concentration varies linearly from
rivatives.Upstream weightedmethodsfor the transport equa- top to bottom, the densitydoes as well. If the densitychanges
tion solution are not used, for reasonsdiscussedin Appendix linearly, the pressuremust change quadratically in order to
A.
maintain hydrostatic equilibrium. However, the approxi-
The variables,pressurep(x, z, t), and concentrationC(x, z, mation for pressureallows only a linear change across the
t), are discretizedas a piecewisecontinuousfunction in space region and a truly linear pressurechange exists only under
using bilinear basisfunctionson quadrilateral finite elements. conditions of constant fluid density. Combining a linear
Their value at any point (x, z, t) is assumedto be constantin changein density with a linear changein pressurein relation
the y direction, as the model is two dimensional.This is ex- (15) would resultin an upward velocitycalculatedat pointsin
pressedfor a meshwith NN nodesas the upper half of the region and a downward velocity in the
lower hall A zero velocity exists at the region centroid, how-
p(x,z, t) • • p,(t)cki(x,
z) (13) ever, as the densityand pressuregradient are consistentat this
i=l
point.The artificialvelocities calculated Withina regionof
NN
spacebecauseof inconsistentapproximationof pressuregradi-
C(x,z, t) • • C,(t)ck,(x,
z) (14) ent and density would dispersea sharp transition zone even
i=l

under hydrostatic conditions. Only employing PAVGat all

where p•(t) and C•(t) are the values of pressureand con- points in the region, in lieu of p(C), would yield a calculated
centration at node i, and tp•(x,z) is the basisfunction defined value of zero vertical velocity throughout.
at node i. The basis function also is defined to take on a
As a result of this inconsistency,artificial velocitieson the
constant value for all y at a point (x, z) and is given by the order of hundredsof metersper year may easilyarisein simu-
expressionsin Appendix B. The weighted residual numerical lation of common field situations.In the particular casewhere
method is describedin Appendix A. the transition zone is narrow, artificial velocities may have a
Velocities must be evaluated throughout the simulated significanteffect on the transport solution. This is especially
regionfor anynumerical method.For themethodin Appen- true for long-term or steady state simulationsin which very
dix A, velocities must be evaluated within each finite element small velocitiesstrongly impact the solution. Using the con-
in order to carry out the integrations that arise. In particular, vention of isoparametriccoordinates(•, r/) definedin Appen-
when using Gaussian integration, values of velocity are re- dix B, the magnitude of artificial velocity in a finite element
quired at each Gausspoint to completethe integrationof the may be calculated at a Gauss point (• = r/= -T-0.577,as de-
solute advection term. Velocity values are also required to fined in Appendix A). The value is v = 0.2885(k[g[Ap/e#) or
calculate the velocity-dependentdispersion tensor. Fluid ve- v -* 2830 (kAp/e)m/s,where Ap is the differencein densityfrom
locities basedon the piecewisecontinuouspressureof relation top to bottom of the element.For example,when permeability
(13) are continuouswithin each finite elementand discontinu- corresponds to that of siltysand,k = 10-• x m2 with porosity
ous at element boundaries.This is due to the discontinuity of e = 0.1, and fluid in an element changesfrom freshwaterat the
discretizedpressuregradientsat elementboundaries.Velocity top to seawaterat the bottom(Ap = 25 kg/m3),the artificial
vi at a point (x•, z) in an elementis based on Darcy's law velocity is about 220 m/yr. In simulating such a systemwith
(equation 3) and is usually given by [e.g., Pinder and Gray, an inconsistentmethod, a sharp transition zone initially con-
1977,pp. 171-173; Huyakorn and Pinder, 1983,pp. 197-205] tained in one element could artificially begin to broaden at a
4. rate of 440 m/yr, when no other natural velocitieswere super-

vi=--(e•).
(k•pkVtPk
+p,I•/IVE)
(15) imposedon the system.Moreover, the velocity,dependentdis-
persion coefficients,relations (7)-(12), take on large values
where the k subscript ranges over the four nodes in the ele- based on the artificial velocitiesand could dispersethe sharp
ment and E is the elevation of the point above datum. transition zone. This is clearly unacceptablefor steady state
A functional inconsistency in this approximation which simulation, and, depending on the geometry and scale of the
gives rise to artificial vertical velocities occurs in the case system, could cause serious errors in simulations lasting a
where both pressureand density are allowed to vary linearly month or more.
(or have the same order in space) in the vertical direction Two examplesof the effectsof inconsistentvelocityapproxi-
acrossan element. Thisinconsistency maybe the reasonfor mations on simulation results are offered in this paper. One
previous difficultieswith finite element simulators solving for example of a regional aquifer crosssectionin southern Oahu,
pressure and concentration which employ velocities that are Hawaii, is discussedlater. The other example, discussedhere,
1854 Voss AND SOUZA' SIMULATION OF NARROW TRANSITION ZONES

120
relates to local scale simulation of saltwater upconing beneath I I I I I I I i I

a pumping well. Reilly and Goodman [1987] compared sharp

interface and dispersed interface methods for analysis of the 100 ,,• -

upconing problem. Their "comparison one" relates sharp and

dispersedsteady state interface results for a problem of radi- • 80","',,
ally symmetric upconing in a homogeneous,anisotropic aqui-
fer cross section with steady pumpage from a partially pen- O
> 60- •', "",,
etrating well. The sharp interface solution used is that of Ben-
nett et al. [1968, caseC-5] from electrical analog analysis.The
solution with dispersion was obtained using the SUTRA
model [Voss, 1984]. Results from a version of SUTRA with
inconsistentvelocity approximations are compared here with
those of Reilly and Goodman[1987] using the consistentap-
proximations. 00 100 200 300 400 500 600 700 800 900 1000
HORIZONTAL DISTANCE, IN METERS
A 20-m long partially penetrating well centered 15 m below
the aquifer top (Figure 1) pumps out 56 kg/s of water. Re- Fig. 2. Comparison of saltwater upconing simulations at steady
state. Solid curve, consistent velocity approximation [from _Reillyand
charge is distributed along the top of the section(1.78 x
Goodman,1987]; dashed curve, inconsistentvelocity approximation.
10-5 kg s-• m -2) and also occursalong the outsideradial Concentrations arc given as percent saltwater.
boundary. Along this vertical boundary, a constant hydro-
static pressuredistribution is specifiedcorrespondingto fresh- An artificial velocity would be calculated in any numerical
water above salt water with a 16-m-thick transition zone cen-
method whenever the functional approximation in space of
tered 20 m above the aquifer bottom. Along the bottom, a pressure gradient and density are not consistent.A discrep-
constant pressurecorrespondingto that at the bottom of the ancy would arise both in the velocity term of the fluid mass
outer radial boundary is specified.No flow occursacrossthe balance (equation (4)) and in evaluation of fluid velocity for
axial boundary. Significantsystemparametersare as follows' transport; moreover, the discrepancywould generateartificial
kh= 2.56 x 10- • •m2 velocity-dependentdispersioncoefficients(11) and (12). If den-
sity is allowed to vary spatially in the vertical direction as an
kv= 1.00x 10-•m 2 Nth power polynomial, then pressuremust be allowed to vary
0.20 as an N + l th power polynomial in order that the pressure
gradient be consistentwith density and vary to the Nth power.
•o0 = 1000kg/m3 For finite elements and a linear density-concentration re-
lationship, this would be satisfiedby a choice of linear basis
= 1024.99kg/m3
•Osaltwater
function discretization for concentration and quadratic basis
% -- 1.0 m function discretization for pressure.Quadratic basisfunctions,
however, significantlyincreasecomputational expenseand an-
0.5 m
other approach may be preferable.
where kh and kv are horizontal and vertical permeabilities. A consistentmethod has been developed by Frind [1982a]
Vertical discretization in the finite element mesh (930 nodes, in a relatively economical finite element model that restricts
870 elements) is 4.0 m. Horizontal discretization increases the mesh to simple rectangular linear elementswith one set of
from a spacing of 0.4 m to a constant spacing of 50 m past a parallel sides aligned with both the gravity direction and the
radius of 100 m. The well radius is 0.35 m. principal permeability directions. Consistency under these
The steady state concentrationdistributionsfor the upcon- conditionsis obtained by choosingto assignonly one velocity
ing problem with consistentand inconsistentvelocity approxi- per element. This value is calculated at the element centroid
mations are shown in Figure 2. The inconsistentvelocity ap- based on the average density for the element. Further, the
proximation widens the transition zone somewhat and causes term in the mass balance derived from the density-gravity
a spatial instability in the numerical solution. Thus use of the term of Darcy's law is approximated using an average con-
inconsistent velocity approximation may lead to incorrect centration for each element which again givesan averageden-
analysisof sucha problem. sity for the entire element. While consistencyis indeed ob-
tainedin this way, the order of approximationin the horizon-
tal direction of density and velocity is also decreased.This
1000 M
may be a significant limitation in horizontally long elements
as described below.
A desirablemesh design for cross-sectionalsimulation of
SCREEN regional fresh-water-saltwater aquifer systems often requires
FRESHWATER
finer discretization in the vertical direction than in the hori-
120 M

INITIAL AND BOUNDARY zontal direction. In order to adequately discretize con-

TRANSITION ZONE

U16 M• centrations in horizontally trending transition zones between

freshwaterand saltwater in a regional system,more nodes are
SAL TWA TER
usually required in the vertical direction than in the horizon-
tal. This is true even when the aquifer is, for example, 10 km
Fig. 1. Geometry of saltwater upconing simulation (adapted from wide and only 100 m thick. Thus to maintain computational
Reilly and Goodman[1987])(not to scale). economy, elementswith high aspectratios (e.g., 1000 m wide
 VOSS AND SOUZA: SIMULATION OF NARROW TRANSITION ZONES 1855

and 1 m high) must often be employed in such regional analy- finite element meshesand on irregular meshescontaining non-
sis. Numerical tests of a centroid-only velocity consistent rectangular finite elements.This is a significantmodification
method, similar to that of Frind [-1982a], show that horizontal to standard Galerkin finite element methodology that results
element-to-elementoscillationsin velocity can occur in meshes in a consistentvelocity approximation.
with large element aspect ratios. The oscillations are
VERIFICATION OF VARIABLE DENSITY MODELS
characterized by a velocity which points upward in one
column of elements and downward in the next, giving rise to The secondcomponent of a successfulapproach to simula-
significant errors in concentration. The instability is exacer- tion of variable density fluid movement is verification that a
bated in regional simulations where the transition zone be- code accurately and stably representsthe physicsimplied by
tween fresh and saltwater is only few elements wide in the the set of governing equations. In order to simulate narrow
vertical direction. While this behavior may be managed or transition zones between freshwater and saltwater, verification
even eliminated by useof a fine horizontal discretization,some of a consistent velocity approximation must be obtained. In
of the advantage of the economical finite element approach particular, a simulator must accurately represent bulk flows
that uses one density value per element would be lost by driven by the buoyancy forces that arise in variable density
requiring the use of both fine vertical and horizontal dis- fluids. The discussion below considers the popular "Henry
cretization on regional scaleproblems. problem" [Henry, 1964] verification, presentstwo test prob-
The restriction requiring fine horizontal discretization may, lems to verify consistency,and presentsthe "Elder problem"
however, be circumvented in a consistent method. The re- [Elder, 1967] for verificationof purely bouyancy-drivenflows.
quirement of consistencygoverns the choice of approxi- Varible density transport models are typically verified by
mations only in the coordinate direction parallel to the direc- comparisonwith the Henry [1964] approximate analytic solu-
tion of gravity. In other directions the density-gravity term tion for steady state saltwater intrusion. However, no model
drops out of both the fluid massbalance (equation (4)) and the to date has successfullymatched the Henry solution, and the
velocity calculation based on (3). Thus within a region of test is not, in general, a sufficientverification for variable den-
spaceor finite element,a higher-order approximation for den- sity transport models.The Henry problem involvesfreshwater
sity may be employed horizontally rather than vertically. in a confined aquifer dischargingto a vertical open sea bound-
The implementation of vertical-only consistency is best ary over a diffuse wedge of saltwater that has intruded the
demonstrated for finite elements on a rectangular element aquifer. A number of numerical models based on significantly
with two sides aligned with the gravity direction. The ap- different methods give nearly identical results for the Henry
proach would be to first calculatea density which varies bilin- problem. These include a particle tracking model by Pinder
early within an element accordingto C(x, z, t) and (2) and then and Cooper [1970], finite element models by Segol et al.
to decreaseits variability only in the direction of gravity by [1975], Huyakorn and Taylor [1976-], Desai and Contractor
averaging the density values vertically in the element. The [1977], and Frind 1-1982a], a finite difference model by
averageddensity in an element varies linearly transverseto the IN TERA [1979], and the U.S. Geological Survey finite ele-
gravity direction and is constant along the gravity direction. ment model SUTRA [Voss, 1984]. None of these results
The vertical velocity based on this density would vary linearly match Henry's analytic solution. This may indicate some inac-
transverseto the direction of gravity and would be constant in curacy in Henry's results,possiblydue to missinghigher-order
the element along the direction of gravity. Horizontal veloci- terms which were originally dropped for the sake of reducing
ties would be unaffected by the density averaging. Such an computation time. While verification of a model by exact
approximation gives a higher-order functional representation comparisonwith Henry's resultsis thus not presentlypossible,
for the velocity field in a finite element mesh than does allow- some confidence in the accuracy of a particular model for
ing just one velocity (and one density)value per element. problems with highly dispersed transition zones may be
In the most general case of a bilinear quadrilateral finite gained if it matches the results of the above-listedsuite of
elementoriented arbitrarily to the gravity vector,the nature of models.
a consistent approximation is necessarilymore complex. One Results of the SUTRA model are compared with results of
complexity results from the variability of density in two di- other authors for the Henry problem. Figure 3 describesthe
mensions within each element. A second complexity results
because the gravity vector is not necessarilyparallel to any
sides of a given finite element. In the local coordinate system 1.0•VIETER
within an element (definedin Appendix B), the gravity vector
may change in both magnitude and direction from point to
point, although it is invariant in global coordinates.A general
finite element method may be developed by beginning with
the assumption that if the numerical pressuregradient and
density-gravity terms are consistently approximated in the
local coordinate systemin which basis functions are defined,
then the terms will remain consistent after transformation to
the global coordinate system.Thus consistencymust be en- //////// • ///////// ///////////////////
forced separately in local coordinates for each vector compo- o.o o.1 1.o 2.0 METERS

nent of the density gravity term. A general finite element

method describedby Voss [1984] is given in Appendix B and
is implemented in the SUTRA model [Voss, 1984]. The con- Fig. 3. Boundary conditionsand SUTRA finite elementmesh for
sistent method has been verified on both regular rectangular Henry [1964] solution.
 1856 Voss AND SOUZA' SIMULATION OF NARROW TRANSITION ZONES

physical systemand finite element mesh (231 nodes,200 ele-

ments) employed by the SUTRA simulation. Velocity- o • SUTRA (100.0min, /' /,,• ]
dependentdispersionwas neglectedby Henry [-1964]in favor • 0.8 D=6.6x10-Sm2/s)
'-- • I
of a constant dispersioncoefficientin order to make a semi- z• 0.6 Huyakorn andTay,or •/ • ///'
(1976) (Steady state,
._•
•

analytic solution tractable. Dispersion is modelled in the nu- <•

• 0.4
D=6.6
SUTRA
x 10-Sm2/s)
(1•.0min,
//
• • •
I
merical simulationsby usinga large constantvalue of molecu-
lar diffusivity and zero dispersivity. o 0.2

The comparisonbetweenvarious models and the Henry

0.0 I I I • I I I I I > x
solution[Henry, 1964] must be consideredin differentgroups 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
DISTANCE, IN METERS
as the variousauthorshavepublishedeithera differentrepre-
sentation of the solution or have used different parameter Fig. 4. Match of isochlorsalong bottom of aquifer for numerical
values. Following most authors, the parameterschosen to resultsof Huyakorn and Taylor [1976] and SUTRA.
match Henry's dimensionlessvaluesare
by 0.1 m for both the SUTRA finite element solution and
0.35 INTERA solution. INTERA simulation was done with cen-
tered in time and centeredin space finite differenceapproxi-
rkg(dissolved
solids)]
0.0357
L kg(seawater)
a mations. SUTRA and INTERA resultsare nearly identical but
do not compare well with Segol et al. and Henry. Moreover,
1024.99kg/m3 the SUTRA results do not approach the Henry solution even
for simulations lasting thousandsof minutes.
-- seawater density
Assigning
a valueof 6.6 x 10-6 m•/s to the moleculardiffu-
sion coefficient would give a total dispersion coe•cicnt of
8p
8C ['(kg
_700 kg(seawater)
(dissolved2m3i]
solids) 2.31 x 10-6 m•/s when the porosityis 0.35. This dispersion
value, rather than that which Henry used appears to have
P0 = 1000kg/m3 been employed by Pinder and Cooper [1970], Segol et al.
Q•v = 6.6 x 10- 2 kg/s [1975], Desai and Contractor [1977], and Frind [1982a]. The
discrepancyseemsto have begun in the work by Pinder and
k = 1.020408x 10-9 m2 (basedon K = 1.0 x 10-2 m/s) Cooper [1970], who used a form of the solute mass balance
101- 9.8 m/s2 equation divided by porosity, whereas in his semianalytical
solution, Henry E1964] did not divide by porosity. Published
0•L -'- 0•T -'- 0.0 results for these models and SUTRA using the molecular dif-
B= 1.0m fusionvalue6.6 x 10-6 m2/sare comparedto the Henrysolu-
tion in Figure 6. Simulatedresultsare nearly identical includ-
Dm-- 6.6 x 10-6 m2/s ing those of Desai and ContracIor El9??] who employed a
two cases
coarsermesh than the others. Becausea total dispersionvalue
Dm--18.8571 x 10-6 m2/s
other than Henry's were used, these solutions should not and
C•v = 0.0 do not comparewith Henry's;however,the solutionsusing
Henry's total dispersionof Figure 5 are an even poorer match
The total freshwaterrechargeis chosenas 6.6 x 10-5m/s with those of Henry.
per meter of cross-sectionthickness(perpendicular to plane of
Steady conditions in saltwater intrusion problems are ap-
section).Two different values of a total dispersioncoefficient
proached rather slowly; therefore Pinder and Cooper E1970]
appear to have been used by the various authors. For the
suggestedthat the solution would convergeto Henry's E1964]
chosentotal rechargerate and Henry's [1964] nondimensional
from opposite sideswhen different initial conditions are used
factor of total dispersioncoefficientdivided by rechargeequal
for the problem. Figure ? confirmsfor the caseof total disper-
to 0.1 the total dispersion
coefficient
mustbe 6.6 x 10-6 m2/s. sioncoe•cient equalto 2.31 x 10-6 m2/sthat the samesolu-
Because the total dispersion coefficient for the SUTRA model
(defined as the sum of molecular diffusion and mechanical
dispersionin equation (1)) in the caseof no velocity-dependent 1.0YMETER
dispersion,is given by a product of porosity and the molecular
diffusion coefficient, the diffusion coefficient should be set to
18.8571x 10-6 m2/swhen the porosityis 0.35. This givesa
simulated result for the concentration at the bottom of the
aquifer after 100 min as shown in Figure 4, where a compari-
son is made with steady state results from three different
models of Huyakorn and Taylor [1976]. Tests of longer simu-
lations with the SUTRA model show that concentrations do
not change appreciably after 100 min. Thus the curves are in
excellentagreement. x

0.0 1.0 2.0 METERS

The same parameters give a concentration distribution in
space as shown in Figure 5. Here, the Henry [1964] solution
Fig. 5. Match of isochlor contours for Henry [1964] analytical
for the 0.5 isochlor is compared with results of $egol et al. solution (for 0.50 isochlor) (long-dashedcurves); INTERA code solu-
[1975] and results from new simulations using SUTRA and tion (short-dashed curves); Segol et al. [1975] (dotted curve); and
the INTERA [1979] transport code. Mesh blocks were 0.1 m SUTRA solution for (0.25, 0.50, 0.75) isochlors{solid curves).
 VOSS AND SOUZA' SIMULATION OF NARROW TRANSITION ZONES 1857

y
1.0 METER layer of saltwater' and second, steady state simulation of the
same system with open vertical sides and uniform horizontal
flow. The first additional test simply checks for consistency
under hydrostatic conditions with a stable density configura-
tion and no flow allowed acrossany boundaries. The choice of
geometric scale and hydraulic parameters is not important
because the system goes to steady state. Longitudinal and
transverse dispersivitiesshould be set to values equivalent to
the length of the largest finite element, and diffusivity should
be set to zero. The correct solution is obtained only if the
pressure gradient and density-gravity terms are consistently
0.0 1.0 0.50 2.0dETERS
approximated. In the correct solution, the transition zone
should remain fixed in one row of elements. The second ad-
Fig. 6. Match of 0.50 isochlor contours for Henry [1964] prob-
lem with simulatedresultsfor Dm= 6.6 x 10-6 m2/sof Pinderand ditional test checks for consistencyin a system where flow is
Cooper [1970] (short-dashed curve); Segol et al. [1975] (dotted parallel to the transition zone and the mesh. The simulation
curve); Frind [1982] (long- and short-dashed curve); and Desai and should be carried out with a longitudinal dispersivity of half
Contractor [1977] (long-dashed curve). SUTRA results at isochlors
the horizontal mesh spacingand a zero transversedispersivity
(0.25,0.50,0.75)(solidcurve).Henry [1964] solutionfor Dm= 18.8571
x 10-6 m2/s(0.50 isochlor,dashed-dottedcurve). and diffusivity. The flow may be set up by creating a horizon-
tal pressuredifference between the hydrostatic pressure con-
tion is approached by the SUTRA model from both sides; ditions that are specifiedalong the vertical slides. The correct
however, this solution is again not Henry's. The initial con- steady state solution maintains the transition zone in the same
dition for the curve approaching from the left is an aquifer single row of elements. Inconsistent approximations result in
filled with saltwater; curvesfrom the right begin with an aqui- spreading of the sharp interface in both cases.
fer filled with freshwater. The numerical resultsof Elder [1967] for a problem of com-
Through the preponderanceof simulators that have closely plex natural convection are further suggestedas a basis for
matched results for the Pinder-Cooper version of the Henry verification of transport simulators.The primary value of this
problem [Henry, 1964] this has becomea standard test for all test is to verify the accuracy of a simulator in representing
variable density transport models. However, because of the bulk fluid flow driven purely by fluid density differences.Elder
unrealistically large amount of dispersion introduced in the dealt with thermally driven convectionin a cross-section.The
solution by the constant total dispersion coefficient,this test solute analog to his "short-heater problem" is a closed rec-
does not check whether a model is consistent or whether it tangular box, modeled in crosssection,with a sourceof solute
accurately representsdensity driven flows, nor does it check at the top implemented as a specifiedconcentration boundary
whether a model can representfield situations with relatively condition with a value of one (Figure 8) and concentration
narrow transition zones.In fact, Henry problem resultsare the specified with a value of zero along the entire base. A zero
same for both inconsistent and consistent approximations. value of pressure is specified at each upper corner, and the
Thus when the transition zone is broadly dispersedand con- pressureis initially hydrostatic.Solute enters the initially pure
centrationgradientsare low, the artificial velocitiescreatedby water by diffusion, increasesits density, and thereby begins a
the inconsistentapproach are insignificantcompared with the circulation process. The following data is used to match
field velocities,and the artificial components of the velocity- Elder's dimensionless results'
dependentdispersionare small relative to the total dispersion. •=0.1
The Henry problem is thus a verification of a variable-density
transport model only for simulation of highly dispersedtransi- k = 4.845 x 10 -13 m 2
tion zones.
#= 1.0 x 10-3 kgm -1 s-•
Two additional tests are suggestedto check for consistency:
first, steadystate simulation of a completelyclosedhorizontal Cinitia
I= 0
aquifer containing a horizontal layer of freshwater above a
p- 1000 kg/m3 + (200)C
y I1- 9.81 m/s2
1.o METER

Dm= 3.565x 10-6 m2/s

•ZL--- •ZT= O

NotethatPmax(C
= Cmax)
= p(C= 1) = 1200kg/m3.
Elder [1967• employed a standard finite difference repre-
sentation of the governingequationsfor vorticity, stream func-
tion, and thermal energy balance. The equations were solved
sequentiallywith a forward temporal difference.The matrix
equations were solved sequentially with a forward temporal
0.0 1.0
x

2.0 METERS
difference.The matrix equationswere resolvediteratively and
an overrelaxed alternating direction method [Elder, 1966• was
Fig. 7. Approach to solutionfrom initially fresh and initially salty used for the vorticity equation. Eldefts grid consistsof 80
conditions for Henry [1964] problem using SUTRA. nodeslaterally and 40 nodesvertically,and time discretization
1858 Voss AND SOUZA: SIMULATION OF NARROW TRANSITION ZONES

< 300
m

•
,,"P=O.O
IIIIIII IIII • • • ' • • • • • • • ''•? t t !! / t '' •1//I/I/I//11/•
P=

T
150 rn

IIIIIIIiIIIIIIIIIIII1/I•1 /IIIIIIIIIIII/I//I/I/
C=O.O

t 600
rn
Fig. 8. Boundary conditionsfor Elder [1967] problem.

was such that it took 20 steps to reach the first point of the complex density-driven flows and solute transport behav-
comparison (1 year elapsedtime) discussedbelow. ior. As with the Henry problem, comparison of results of a
The problem is simulated with the SUTRA model, using a given simulator with other numerical solutions lends confi-
mesh coarser than Elder's [19671, consistingof 44 finite ele- dence in the accuracy of the simulator.
ments horizontally and 25 elementsvertically (1170 nodes and
PRACTICAL DISCRETIZATION FOR NARROW TRANSITION ZONES
1100 elements). The simulation is carried out with time steps
of 1 month (2.6298x 10a s), and iterationsfor eachtime step The third component of a successfulapproach to obtaining
are considered complete when the maximum absolute change accurate and stable simulations of variable density transport
for all nodes from the previous iteration falls below 500 Pa for problems is the careful choice of spatial and temporal dis-
pressure and 0.01 for concentration. While the simulation cretization. Spatial discretization is sufficient if it (1) ad-
could theoretically be carried out on only the left or right half equately describesspatial variations in aquifer parametersand
of the system,modeling the entire box checksthe ability of the representsflow and transport processesand (2) provides nu-
simulator to produce symmetric vortices. merical stability. Relatively fine discretization must often be
Figure 9a compares nondimensional concentration results used in directions both along flow lines and acrossflow lines.
of simulation with SUTRA to nondimensional temperature A desirablemeshdesignoften employselementsthat are much
results of Elder [1967]. Figure 9b compares the flow field longer horizontally than vertically, as is discussedin a pre-
simulated by SUTRA with the steamline results of Elder vious section.Practical mesh designis a somewhat qualitative
[1967]. The concentration field is strongly coupled to the subject and appropriate choice of element size is discussedin
complex flow field which evolves through a seriesof transient the following.
anticlockwise and clockwise vortices. Although SUTRA uses Along flow lines a mesh Peclet number Pem,which varies
somewhat coarser spatial and temporal discretization, the re- from point to point, may be defined depending on local mesh
suits compare very well, spatially and through time, showing size,fluid velocity, and dispersioncoefficients.The mesh Peclet
that both numerical solutions give similar representation of number is given as
IvlasL IvlasL
m 1 year 10 years Pe• - - (16)
DL (Dm+ z,.Ivl)
0 , , I 1• / I .
5O
rn 1 year 10 years
•.4' - ' •...-.• •

x/.... '••lf'•x..
-_-
xi)
100
50 -
150 • i • i i i
2 years 15 years 100

150 , , ' ' '

50
15 years
0 2 years

100
t
150 • I I I I I
4 years 20 years •00 ' • •7
150' " '• '
4 ye•s 20 years
0

5O

15o' , , / '•- -',.--•o

•_ ,; / 100
100 j 300 400 500 100 200 300 400 500 m
15O
100 2• 300 400 500 100 200 300 400 500 rn
Fig. 9a. Comparisonof nondimensionalresultsfor Elder [1967]
problem, showing 20 and 60% of maximum concentrationsfor Fig. 9b. Comparison of flow fields for Elder [1967] problem,
SUTRA (solid curves) and of maximum temperature from Elder showing every seventhelement centroid velocity for SUTRA and
(dashed currves). steamlinesgiven by Elder.
 VOSS AND SOUZA: SIMULATION OF NARROW TRANSITION ZONES 1859

where AsL is the local distancebetween sidesof an element flow, the numericalsolutionis quite insensitiveto the value of
measuredalong the local flow direction,Iv[is the magnitudeof transverse dispersivity,aT. The true effectof a low valueof aT
local velocity; and DL is the total longitudinal dispersioncoef- is almost always lost in the artificial dispersioncausedby
ficient as defined in (16). This Peclet number is a measure of large transversemeshspacing.
the amount of local advective transport relative to the local Choiceof discretizationin time, followinga changein stress
amount of diffusiveand dispersivetransport.When, as in most on a system,shouldgiveearly time solutionsa smalltime step
fieldcases,moleculardiffusionDmis smallrelativeto mechani- and late time solutions a larger time step. Time steps for
cal dispersion,D L • %[v[,and the meshPecletnumberis density-dependenttransport simulation must be small when
beginninga simulationfrom poor initial conditionswith sharp
Pem• (17) concentration gradients.A small time step for the transport
\%/ solution may be definedas that which would allow a particle
When the mesh Peclet number is high, centeredfinite differ- of fluid to traverse on the order of a tenth of the distance
ence and Galerkin finite element approximations can give across any element. A large time step at early time would
values of concentration which oscillate in space. For the nu- allow fluid to traverseon the order of one elementper step.
merical method described here, oscillations do not occur when Restrictionson time step size for stability of transport solu-
the mesh Peclet number is approximatelyfour or less(Pe,,,< tions are most likely to be found in meshareasof converging
4). Thus in designing a mesh, discretization must be chosen so flow, suchas near wells.This is becausevelocitiesare high and
that along the direction of flow, 4aL > AsL. Should this re- mesh sizes likely small in those areas. At late time, when con-
quirement be satisfied everywhere in a mesh except in small centration gradients tend to be less sharp, much larger time
localized regions,the solution will likely oscillate only in these stepsmay be taken without sacrificingeither numerical accu-
regions and may not spoil the rest of the solution, resulting in racy or stability.
an acceptable solution to the problem. If oversized elements
are used where longitudinal concentration gradients are low, EXAMPLE: SOUTHERNOAHU, HAWAII
oscillations may not occur at all and the criterion may be Simulation of groundwater flow and movement of dissolved
safely violated by more than an order of magnitude. This solids in the southern Oahu aquifer near Honolulu, Hawaii, is
economy is often possible in cross-sectionalanalysis of seawa- a difficult problem, as the transition between freshwater and
ter intrusion in aquifers, as is shown in the Hawaii example saltwater is narrow inland, although it is broad near the coast.
that follows.
Becauseimportant hydraulic boundaries are relatively closeto
In the direction perpendicular to flow lines, only dispersive areas of interest near pumping centers,the entire aquifer must
transport takes place. Thus spatial discretization in the direc- be simulated, including both narrow and broad portions of
tion perpendicular to flow lines must be sufficientto represent the transition zone. The SUTRA model [l/oss, 1984] is ap-
a solution to the equation describing a diffusive process in plied to the southern Oahu system on a well-discretized but
that direction.The total transversedispersioncoefficientD T is practical mesh in order to demonstrate the efficacy of the
made up of a molecular diffusion component and a mechani- suggestedmodeling approach. Comprehensive description of
cal dispersion component, D T = (Dm+ aT[V[).In the steady hydrologic modeling of the southern Oahu aquifer is given by
transport casewhere D T is zero, with any finite value of total Souza and l/oss [1987]. In this section, a consistentvelocity
longitudinal dispersioncoefficientDL the analytical solution approximation is demonstrated to be a vital factor in suc-
gives constant concentration along each streamline. In partic- cessfulsimulation of the southern Oahu aquifer.
ular, for freshwater and saltwater under these conditions, the The southern Oahu basaltic aquifer is composed of thinly
fluids do not mix at all and are separated by a perfectly sharp bedded lava and rubble layers (thicknessesfrom 1 to 10 m) on
interface. In the simulation of these conditions with Pe,,,< 4 the southwestflank of the Koolau range. Regionally, the aqui-
and a flow field not parallel to the sides of a mesh composed fer is very thick (-• 1800 m) and quite permeable (K-• 500
of regular rectangular elements,the numerical solution gives a m/day). The aquifer fabric may be expected to have vertically
transition zone with a thickness of a few elements. Thus the anisotropic permeability between 10:1 and 10000:1 due to
sharp analytical solution for zero aT is not possibleto obtain the layered structure, although the regional anisotropy value
numerically as transversediscretizationis never fine enough to has not been establishedin previous studies.The basalt dips
represent the sharp interface solution. This error disappearsas down toward the southwestat an angle of a few degrees.Near
AsT/D T decreases,where AsT is the transversespacing.Frind the coast and below Pearl Harbor (the involuted inland water
[1982c] has suggestedalleviating this problem through mesh body of Figure 10), the basalt is semiconfinedby a wedge of
orientation that parallels streamlines. interlaced sediment and coral formations, referred to as the
In order to determine the amount of transverse dispersion "caprock" (Figure 11).
caused by discretization, a steady state transport simulation Primary recharge occurs over the Koolau mountains where
for any given problem may be run with aT = 0 to obtain the significant intrusions of vertical dike walls into the horizontal
narrowest transition zone possiblewith a given flow field and beds make horizontal flow difficult. The dike-intruded region
given mesh. A rule of thumb for accurate transverse dis- is considered a no-flow boundary for the sake of the aquifer
cretization in cases with velocity-dependent dispersion is to model that contributes a large specifiedrecharge to the basalt
choose the mesh spacing such that AsT • aT, although nu- aquifer due to spills over and around dike walls. A general
merical experimentsshow that AsT • 10 aT may be adequate description of hydrogeology is given in the work by l/isher
for many situations. This stringent criterion has been violated and Mink [-1964-[.
in most published transport simulations,calling into question Discharge in the undeveloped state of the aquifer occurs
the validity of simulation-fitted values of transverse dispersi- primarily at springs along the inland edge of the caprock.
vity. For a mesh that is coarse in the direction transverse to Some diffused discharge is expected as upwards leakage
1860 Voss AND SOUZA' SIMULATION OF NARROW TRANSITION ZONES

Approximate
orientalion
ofcros s-sectional

w•TEIa
Boo'/
HiGH-LEVEL....•....e_
•'
-1'

•[ sOUTHERN
O

ool c^PRoc,<

\
\

pectftc
0
Oee•q
10 KILOMETERS

Fig. 10. Southern Oahu aquifer map with boundaries and 1958 hydraulic heads adapted from Mink [1980].

through the caprock into Pearl Harbor. Total throughflow in and a portion of the caprock similar to that shown in Figure
the aquifer is on the order of 106 m3/day (255 Mgal/day) 11. The finite element mesh for the cross section (702 nodes,
[Mink, 1980]. 646 elements) has a portion with very fine vertical dis-
A cross-sectional model is highly appropriate for repre- cretization (15.24 m) in the region of the aquifer where the
sentation of systembehavior becauseof the regionwide nearly transition zone between freshwater and saltwater is expected
uniform flow field in which lateral convergenceis negligible. to reside (Figure 12). This allows sufficient discretization to
Moreover, pumping is spreadin a narrow band approximately represent the thin part of the transition zone. Below this
parallel to both contours of equal hydraulic head and to the region, concentrations, velocities, and pressureschange very
surface-caprockcontact near the edge of Pearl Harbor. Initial slowly in space and coarse discretization is sufficient. Sharp
hydraulic head in the system (predevelopment, 1890) was pressurechangesand converging high fluid velocitiesoccur at
10-13 m, leading to a 400- to 520-m-thick freshwater zone the springsnear the edge of the caprock; thus horizontal dis-
above the seawater interface estimated from the equilibrium cretization is finest in this region (304.8 m). Although regional
Ghyben-Herzberg principle. Present-dayheads are lower as a longitudinal dispersivity is expected to be less than 100 m,
result of significant pumping stressand range from 5 to 8 m. horizontal spacing may greatly exceedthe stability limit of 400
The present-day freshwater lens would be only 200 to 320 m m based on relation (17), because lines of constant con-
thick if it were in static equilibrium with the sea. In fact, the centration will, for the most part, parallel the flow directions.
interface is not sharp everywhere in the aquifer (Figure 11) Here, on the regional scaleunder steady conditions,longitudi-
and it is not clear whether static equilibrium exists under nal transport is less significant than transverse transport in
current pumping conditions. Thus a variable density model defining the concentration distribution in the transition zone.
capable of simulating transient conditionsis required to evalu- All boundaries of the mesh are closed to flow except where
ate the system. (1) freshwater recharge from the dike zone is specifiedat the
The simulated region includes the dike-free basalt aquifer back of the aquifer, (2) pressureis specifiedto be zero below

Sea Level

•Fi•EIS•W,•rE•?•['*i
?•TWEEN
DIKES
• FRESHWATER
Y CK-L-----=----•----=•
• ....

Fig. 11. Typical crosssectionof southernOahu aquifers(not to scale).

 VOSS AND SOUZA' SIMULATION OF NARROW TRANSITION ZONES 1861

WATERTABLE ' p= 0

-lOO
.....ii:::•::::
........................
:::;::;i.
LOW PERMEABILITY::::::?::::::::::::::::::::::

-200 HIGH PERMEABILITY

::::::::::::::::::::::::::::::::::::::::::::::
::::::::::::::::::::::::::::::::::::::::
Ki;11111•;11:;i;11:,11;i;i;111;11111:::
.....
•:::•:•:•:•:•:•:•:•:•:•:•:•:•:•:•:•:•:•:•:•:_::•:•:•:•:•:•:•:•.
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
[ ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
•
'" -300
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

z -400
:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
o_ -500

c3 -600 ?':::•';•'"•:r:'•
.............
•
- 1000

-1400

-1800

DISTANCE, IN KILOMETERS

Fig. 12. Mesh for cross-sectionalfinite element model of southern Oahu aquifer (not to scale).

Pearl Harbor, and (3) pressureis specifiedto be hydrostatic, where K n and Kv are horizontal and vertical hydraulic con-
based on a column of seawater where the mesh is arbitrarily ductivity,K• is the hydraulicconductivityof the seaboundary,
cut off on the seaward end. A higher impedance is specified and C is the dissolved solids concentration as a mass fraction.
along the seaward end below the caprock to represent the The simulated concentration distribution for predevelop-
effect of the aquifer seawardof the boundary. A parsimonious ment steady state conditions is shown in Figure 13 to scale
model identification uses only six parameters to characterize and in Figure 14 with vertical exaggeration. The transition
the entire system behavior [Souza and Voss, 1987]. These zone is very sharp inland, as expected,and gradually broadens
values, as well as known system constants, define the 1-m- near the caprock. The predevelopment velocity field (Figure
thick (in direction perpendicular to plane of section) repre- 15) illustratesthe saltwatercirculation,includingspring dis-
sentation of the system: charge near the caprock edge, and diffuse dischargethrough
the caprock.
K n = 457.2 m/day
These solutions are impossible to obtain with an inconsis-
Kn/K v = 200 tent velocity approximation (standard Galerkin method) as
demonstratedby comparisonin Figure 16 (for Kn/K• = 20).
K caprock= K n 10-'* Experiments show that neither low values of dispersivitynor
ks= Kn 10-2 high anisotropy (Kn/K• allow a relatively narrow transition
zone to be simulatedwith the inconsistentapproximation.The
Sy= e = 0.04 artificial velocity and resultant dispersion generated by the
• = 2.5 x 10 -9 Pa- x inconsistentapproximation give highly incorrect results,over-
whelming the physical dispersionin the system. However, ac-
Q inflow = 0.404 kg/s curate simulation based on a consistent velocity approxi-
P0= 1000kg/m3 mation provides a fundamental tool for hydrologic analysisof
the systemas describedin the work by Souzaand Voss[1987].
t3p/t3C
= 700 kg/m3
DISCUSSION
/• = 4.47 x 10-xø Pa-x
Analysis of a freshwater-saltwateraquifer system that con-
/• = 1.0 x 10-3(kg m- x s-x) tains a relatively narrow transition zone is a difficult problem
for variable density transport simulation. However, such
% = 80 m analysis is desirable as narrow transition zones are often en-
•r = 0.20 m countered in the field. Most published analyses either deal
directly with the case where data indicates a broad transition
C seawater -- 0.0357
zone, or give simulation resultsthat show a broad zone. Based
p seawater
- 1024.99kg/m3 on the typical difficulties expected in modelling variable den-

5 5O 95
o• 0 I I I• , • .............................
"•.'•2.1
::ii::::!•i::
i::•::i::i::i:.::i::
i::
iiii•l:.::ii
:::::::::::::::::::::::::::::::::::::
:::::::::::::::::::::::
i::::i
!iii
:•i
:•i
iiii:•iii!ii•iii:•iiii:•iii:•i
iiiiii:.•;•:•!ii!•iiiiiiiii::i!11
iliiiiiiiii
ilfiiiiiiil•iiiiiii
iiiiiiiiiiii
ii

a. ,,, 900

_z 1800
0
I I 9I I I 12I I I I I I 15 18
DISTANCE, IN KILOMETERS

Fig. 13. Simulated predevelopmentdistribution of total dissolvedsolids (as percent seawater)for southern Oahu aquifer
(to scale).
1862 Voss AND SOUZA' SIMULATION OF NARROW TRANSITION ZONES

1 O0 - 100

200 - 200

300- ..• 300

25
400 50 • 400
500 90 500

600
t
700
0
I
1
I
2
I
3
I
4
I
5
I
6
I
7
I
8
95
I
9
i
10
DISTANCE, IN KILOMETERS
i
11
i
12
i
13
i
14
i
15
i
16
i
17
t 600
18
700
0 1 2 3 4 5 6 7 8 9 10 11
DISTANCE, IN KILOMETERS
12 13 14 15 16 17 18

Fig. 14. Transition zone as simulated for predevelopment con- Fig. 16. Comparison of simulated predevelopment transition
ditions (concentration as percent seawater) for southern Oahu aquifer zone for Southern Oahu aquifer between model based on inconsistent
(not to scale). (dashed curves) velocity discretization and consistent (solid curves)
discretization (concentration as percent seawater)(not to scale).
sity transport, speculationmay be warranted with regard to
whether these analyses have accurately represented variable herein) gives a consistent velocity approximation for bilinear
density transport as described by the standard governing quadrilateral elements. The modified Galerkin method em-
equations. Errors in the simulation of the dynamics of a ploys pressure gradient and density gravity approximating
narrow transition zone stem from any inconsistencies,inaccu- functions that are consistentin the local coordinate system of
racies or instabilities inherent in a given model. The likely each element.
sources of simulation error are threefold. First, inconsistent 2. Three tests in addition to the Henry [1964] seawater
approximationsof terms involved in the fluid velocitycalcula- intrusion problem should be carried out to verify a simulator.
tion can lead to large artificial velocity and dispersioncompo- The Henry problem serves only as a basis for verifying a
nents. Second,inaccuracycan occur in representationof bulk variable density transport code in the case of highly dispersed
fluid flows driven by density differencesin the fluid. Tests transition zones. The test gives nearly identical results with
typically used to verify variable density flow and solute trans- inconsistent and consistent formulations; thus the inaccuracies
port models do not verify that a given model is free of the of inconsistent simulators can not be detected. Further, the
above sources of error. Third, discretization is typically too test only weakly verifies accuracy in simulating bulk flow
coarsefor the desiredlevel of transversedispersion. driven by density differences.Two of the additional tests are
The following modeling approach alleviates these difficul- simple checks for consistencyof the velocity approximation.
ties. These simulate a systemwith a narrow transition that should
1. A consistentvelocity approximation must be employed. remain narrow for all long-term and steady state simulations.
The standard Galerkin finite element method gives an incon- The third additional test is the solute analog of the Elder
sistent velocity approximation that can generateartificial ve- [1967] natural convection problem. This verifies the accuracy
locities and dispersionwhich overwhelm physicaldispersionin of simulation of flow driven purely by density differencesin
a simulation. Finite difference methods may or may not be the fluid.
consistent,dependingon the particular interpolations used.A 3. Spatial and temporal discretization must be sufficientto
modified general Galerkin finite element method (presented guarantee accuracy and stability. In particular, vertical dis-
cretization should be on the order of the transverse dispersi-
I__. I
vity value when flow is predominantly horizontal. This re-
lOO quires unrealistically fine discretization in some caseswhich
__.._.- _.-.-..•..--..-'//I ....•....... must be compromised. Thus transport simulation studies
•i::i::iiii•iiii::iiiiiiiiiiiiiiii?•iii?:iiiiiiiiiiiiiii:•l
200 _. _•..•..• •-///•- I :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
should always begin with a steady state simulation for the case
::::::::::::::::::::::::::::::::::::::::::::::::::::::

of zero transverse dispersion. This checks for the narrowest

300
• • ......... ::::::::::::::::::::::::::::::::::::::::::::::::
transition zone possiblewith a given meshand flow field.
400 - The fine gridding required to obtain stable and accurate
solutions when transition zones are narrow leads to expensive
500 - computer runs on field scale problems. While this is unfortu-
nate, the expense of transport simulation is a natural and
600 -
possibly unavoidable feature of the problem. Even numerical
700 i i i solutions by particle tracking or a method of characteristics,
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
which do not have stringent stability restrictions for problems
DISTANCE, IN KILOMETERS
with high Peclet numbers,require fine discretization or a large
Fig. 15. Velocity field as simulatedfor predevelopmentconditions number of tracking particles for accurate approximation of
in southern Oahu aquifer. (Anglesare scaledto show true direction the dispersion term when the dispersivity is small but nonzero.
on exaggeratedmap.) Largestvector= 7.3 x 10-s m/s at direction
7.5ø above horizontal, occurs just below inland edge of caprock.
Thus computational expensemay be unavoidable becausea
Vector lengths are proportional to log of magnitude. Vectors less large amount of information must be borne by any numerical
than 0.25% of maximum velocityare not plotted (not to scale). method to simultaneously define both the concentration dis-
 VOSSANDSOUZA:SIMULATION
OFNARROW
TRANSITION
ZONES 1863

tribution within and movement of narrow finite thickness in whichthe definitionfor fluid masssourceQi(x,z, t), M/T:
transition zones.

In order to guaranteeuseful results,a variable density flow

and solute transport simulator must be based on a consistent i= Qv
=• (A4a)
velocity approximation, must be verified with the before-
has beenemployedand Qm,(x,z, t), M/T, is the fluid mass
mentioned tests, and must be used with proper spatial and sourceat node i along the boundary F of the model region
ß

temporal discretization.An illustrative application of such a

simulator is the analysisof the regional aquifer of southern
Oahu, Hawaii, that contains both narrow and broad sections i:lkAi/
QiN,
• • (QIN,• (A4b)
of a freshwater-saltwater transition zone. Successful evaluation
where A• is the area of model boundary containedin cell i.
of the effectsof newly proposeddispersionmodels in such a
The general boundary condition employed in (A3) has the
system would require analysisbased on this modeling ap- form
proach. Simulation based on this approach makes possible
state-of-the-art analysis of a range of difficult problems at
various scales.
- fr(Spv'n)
•idF=
frQiN*•idF+
fr(•3•)•idF
(A5a)
where n is the outward unit normal at the boundary and the
APPENDIX A: WEIGHTED RESIDUAL METHOD
last term has been evaluated as constant for each cell:
Applying a weightedresidual method to the fluid massbal-
ance, (4) resultsin NN relations: Vi

wherehi(x,z) is effectivecell height,bi = •/A i. In a horizontal

rectangularelementabuttingthe water table, bi is half of the
element height.
+ (V. dr= dr i=i,NN(A1) Discretization of the term involving permeability in (A3)
using (13) gives
where V is volume. An integrated finite difference technique
for the time derivative and source term may be implemented
by assumingthat the complete term involving the time deriva-
tive and the source has a constant value throughout a cell
surrounding the node. The cell extends to the boundary de-
fined by the connectedbisectorsof opposite sides of all ele- = ßv+,ßv+,av
ments contiguous to a node (Figure 17). The volume of the cell
at node i, V•, is givenby

V,=f•4,dV Discretization of the (pg) term is discussedin the consistent

velocity sectionof the text and in Appendix B.
By clearing the integrals of constants,applying Green's the- Discretizing the time derivatives using a backwards finite
orem to the term containing velocity, employing Darcy's law difference approximation
(equation (3)) and the general boundary condition (6), the
weighted residual relation becomes
dpi pin+1 __pin
(A7)
dt Atn+ 1
dCi Cin__Cin- 1
(AS)
dt Atn
+ .(Vp-pg) .V•b, dV=Q,+QiN, (A3) where

Atn=tn- tn-1 (A9)

P7= p,(t.) (AlO)
/z C7 = Ci(tn) (A11)
and employing relation (A9) results in a form of (A3) which
may be solvedfor pressure
at all nodes,pi"+x. All coefficients
except (3C/&) are evaluated by linear projection to the half
time level(n + «) whenonly a singleiterationis allowedper
time step and at the (n + 1) time level for time steps with
multiple iterations. The (3C/3t) term is usually small for sea-
water intrusion and thus is lagged one time step.
The integrations are carried out element by element after
x transforming the quadrilateral element to a local coordinate
Fig. 17. Cells, elements,and nodes for a two-dimensional finite ele- system in which each element is a square with sides of length
ment mesh composedof quadrilateral elements. two. The numerical Gauss integration uses two Gauss point
1864 Voss AND SOUZA' SIMULATION OF NARROW TRANSITION ZONES

locations in each local coordinate direction (Figure 18). At The dispersiveflux acrossboundaries has been assumedto be
each of the four Gauss points in an element,density is evalu- much less than the advective flux in (A13). Integrations are
ated as a function of the bilinearly interpolated concentration carried out as described above, using densitiesbased on con-
C(tn+1), at that point accordingto relation (2). At the nodes, centration,C(tn+1) in relation (2).
density is evaluated based on the nodal concentration value The methods of upstream weighting the advective solute
Ci(tn+1) from relation (2). In the term involving gravity, how- flux term or of employing asymmetric rather than Galerkin
ever, a specialevaluation of the density valuesat Gauss points weighting functions [e.g., Huyakorn and Pinder, 1983] are not
must be made in order to be consistentwith the local pressure employed. These methods have been applied to transport
gradient, as is describedin the main text and in Appendix B. modelsin order to damp spatial oscillationsof concentration
Verification of the water table form of the cross-sectional in mesheswhich are too coarse. They increase the effective
fluid mass balance was made by comparison of numerical local value of longitudinal dispersivityby a factor proportion-
solutions with the general solution of Neurnan [1974]. Al- al to local mesh spacing in the direction of flow. Increasing
though fluid storageis accountedfor at a moving water table dispersivity through these complicated methods is compu-
by this formulation, the position of a water table is held fixed tationally expensive.Further, having dispersivityas a function
in spacethroughout a simulation. This is a practical approach of the numerical discretization size rather than of a physical
as long as true movement of the water table is only a small parameter is confusing when applying a model. A direct ad-
fraction of the total height of the aquifer unit, as is discussed justment of the longitudinal dispersivity has the equivalent
by Neurnan 1-1974]in relation to analytic solution. damping effect and is preferable; however, the correct mod-
The weighted residual method applied to the solute mass eling approach for both finite element and finite difference
balance (equation (1)) yields NN relations models is to design the mesh fine enough to give nonoscillat-
ing solutions at the specifiedphysical value of longitudinal

;v(ePr3-•--)qbi
rSV
+j•v
(epv
'VC)qbi
dV dispersivity.

APPENDIX B: CONSISTENT VELOCITY APPROXIMATIONS

- ;v{V'[eP(DmI
+D)'
VC]}qbi
dV=fv[Qt,(C*
- C)]qbi
dV The requirement of consistencyis applied in each local di-
rection of isopatametric coordinates, • and r/, (Figure 18) as
i= I, NN (12) follows.First, the • componentof the density-gravityterm in
Assumingcellwisediscretizationfor the time derivative and local coordinatesmust have the same spacedependencyas the
source terms, clearing the integrals of constants,employing pressure derivative term in local coordinates c•p/c• in that
Green's theorem on the dispersion term and employing a direction. Second,the r/component of the density-gravityterm
backwardsfinite differenceapproximation for the time deriva- must vary in spacefunctionally as does c¾/c•rl.The basis func-
tions and derivatives are as follows:
tive and discretization(14) for C, resultsin
fll(•, r/)= E_H_ (Bla)
At. + I ._1
-- Cin)'+
NN
•1 n+lfV
j=
Cj [(8p¾)
ßV•)j]•i dg
f12(•, •/) = E +H_ (B lb)

f13(•, r/) = F.+H+ (Blc)

+
j=•1C.•n+
1 [eP(Dml
q-O)-Vlp./]
ßVlPi
dV fl½(•, •/) = E_ H + (Bld)
= Qi(Ci* -- Ci)n+' i = 1, NN (A13) where

whichmay be solvedfor concentration

at all nodes,Cin+l E_(•) = «(1 -- •) (B2a)
E+(•) = «(1 + •) (B2b)

NODE 4 NODE 3
H_ (•) = «(1 - •) (B2c)
(-1 (1,1)
H +0/) = «(1 + •/3 (B2d)
The two-dimensionalbilinear basisfunctionsq•i(x,z) in (13)
GAUSS POINT 3 GAUSS POINT 4 and (14) when defined in the local element coordinate system
x x
(•, r/) are denoted as fli(•, r/), i- 1, 2, 3, 4. There is one basis
(-0.577, 0.577) (0.577, 0.577)
function defined for each node in a quadrilateral element.The
derivatives of the bilinear basis functions depend on only one
spacecoordinate'

(-0.577, -0.577)
x
(0.577, -0.577)
x
.... «H_ •f•l
•r/
= _ñ=
2--
(B3a)
GAUSS POINT 1 GAUSS POINT 2
c•fl2
c•fl2
_ _x= (B3b)
(-1,-1) (1,-1)
NODE 1 NODE 2
-- +«H+ -- +•+1= (B3c)

Fig. 18. Finite element in local coordinate system with Gauss - +«F._ (B3d)
points.
 VOSS AND SOUZA: SIMULATION OF NARROW TRANSITION ZONES 1865

The basisfunctionfor node i is f•i(•, r/). The pressuregradi- where(/5•),,and(/5•)z aretheconsistently

discretized
density-
ent within an element in local coordinates is defined in terms gravity term componentsin global coordinates,and [J-•] is
of the derivativeswith respectto the local coordinates: the inverse of the Jacobian matrix.

c•p 4 c•f•
i (B4a)Acknowledgments.
The authorsgratefullyacknowledge
T. Reilly,
i--1
U.S. Geological
Survey,for helpingto bringthe consistency
problem
to light;F. Lewis,formerlyof U.S.Geological Survey,for analyzing
the Elder [1961] problem;and the research groupheadedby E.
&/(•,q)=• p,•
i= (B4b) Buetow,Technical University
lem as a test case.
of Berlin,for suggesting
theElderprob-
The essence of this consistentmethodis in findinga local
discretization of pg, with a spatialfunctionalitythat is consis- REFERENCES
tent with the local pressurederivatives. This functionmay be Bear, J., Hydraulics of Groundwater, 567 pp., Mc-Graw Hill, New
found by inspectionto be York, 1979.
Bennett, G. D., M. J. Mundorff, and S. A. Hussain, Electric-analog
studies of brine coning beneath fresh-water wells in the Punjab
(pg)½*(•,
if)=• p,g½,
•i=1 (B5a) Region, West Pakistan, U.S. Geol. Surv. Water Supply Pap., 1608-J,
pp. J1-J31, 1968.
Desai, C. S., and D. N. Contractor, Finite element analysis of flow,
(pg),*(•,
rl)=Y',p,g,,
•
i=l (B5b) diffusion and salt water intrusion in porous media, in Formulation
and ComputationalAlgorithmsin Finite Element Analysis,edited by
where the vertical bars indicate absolutevalue; pi is the value K. J. Bathe et al., pp. 958-983, MIT, Cambridge, Mass., 1977.
of p at node i in the element based on the value of C at the Elder, J. W., Numerical experimentswith free convection in a vertical
slot, J. Fluid Mech., 24(4), 823-843, 1966.
nodethroughrelation(2); g½,is the • component
of g at node Elder, J. W., Transient convection in a porous medium, J. Fluid
i; gn,is ther/component
of g at nodei; and(pg)½*
(•, if) is the Mech., 27(3), 609-623, 1967.
consistentrepresentation of the term for the • direction and Frind, E. O., Simulation of long-term transient density-dependent
(pg),* (•, •) for the ff direction.Note that ff dependence
is transport in groundwater, Adv. Water Resour.,5, 73-78, 1982a.
removed for the ff term at a given •, and • dependence is Frind, E. O., Seawater intrusion in continuous coastal aquifer-
aquitard systems,Adv. Water Resour.,5, 89-97, 1982b.
removed from the • term for a given if. Frind, E. O., The principal direction technique: A new approach to
This functionality is equivalent to that of the pressuregradi- groundwater contaminant transport modeling, in Finite Elements in
ent term in each local direction (equations (B4a) and (B4b)). Water Resources,edited by Holz et al., pp. 13/25-13/42, Springer-
No particular significanceshould be attached to the absolute Verlag, New York, 1982c.
values of basis Mnction derivatives, except that these satisfy Henry, H. R., Effects of dispersionon salt encroachmentin coastal
aquifers, U.S. Geol. Surv. Water Supply Pap., 1613-C, pp. C71-C84,
the requirementsof consistency.This discretization of the (pg) 1964.
term (equations (B5a) and (B5b)) is robust in that it allows Huyakorn, P.S., and G. F. Pinder, ComputationalMethods in Subsur-
both the density and local gravity vector components to vary face Flow, 473 pp., Academic, Orlando, Fla., 1983.
Huyakorn, P.S., and C. Taylor, Finite element models for coupled
overan element.The localgravityvector(g½,gn)requiredby
groundwater flow and convective dispersion, Finite Elements in
relations (B5a) and (B5b) may be obtained for each element by Water Resources,edited by Gray et al., pp. 1.131-1.151, Pentech,
a coordinate transformation: London, 1976.
INTERA, Revision of the documentation for a model for calculating
effects of liquid waste disposal in deep saline aquifers, U.S. Geol.
g• Surv. Water Resour. Invest. Rep., 79-96, 73 pp., 1979.
Lebbe, L. C., Mathematical model of the evolution of the fresh water
where [J] is the Jacobian matrix defined by lens under the dune and beach with semi-diurnal tides, Geol. Appl.
Idrogeol., •8(2), 211-266, 1983.
Lee, C-H., and R. T. Cheng, On seawater encroachment in coastal
aquifers, Water Resour. Res., 10(5), 1039-1043, 1974.
Mink, J. F., State of the groundwater resourcesof Southern Oahu,
•x •z • •: • • x• z•
technical report, 83 pp., Board of Water Supply, City and County
of Honolulu, Hawaii, 1980.
Neuman, S. P., Effect of partial penetration on flow in unconfined
•B7• aquifers considering delayed gravity response, Water Resour. Res.,
•0(2), 303-311, 1974.
where(x•,z•) are the globalcoordinatesof the nodesin the Pinder, G. F., and H. H. Cooper, A numerical technique for calculat-
element. ing the transient position of the saltwater front, Water Resour.Res.,
The eightgravityvectorcomponents at the nodesin each 6(3), 875-882, 1970.
elementneed be calculatedonly oncefor a given meshand Pinder, G. F., and W. G. Gray, Finite Element Simulation in Surface
and SubsurfaceHydrology, 295 pp., Academic, Orlando, Fla., 1977.
maybe savedat thestartof a simulation. The valuesof the Reilly, T. E., and A. S. Goodman, Quantitative analysis of saltwater-
consistent density-gravity
termsand pressure derivativesin freshwater relationships in groundwater systems--A historical per-
globalcoordinates areobtained throughanothertransforma- spective,J. Hydrol., 80, 125-160, 1985.
tion: Reilly, T. E., and A. S. Goodman, Analysis of saltwater upconing
beneath a pumping well, J. Hydrol., 89, 169-204, 1987.
Segol, G., and G. F. Pinder, Transient simulation of saltwater intru-
sion in southeastern Florida, Water Resour. Res., 12(1), 65-70,
1976.
(B8)
Segol, G., G. F. Pinder, and W. G. Gray, A Galerkin finite element
technique for calculating the transient position of the saltwater
front, Water Resour. Res., ••(2), 343-347, 1975.
Souza, W. R., and C. I. Voss, Analysis of an anisotropic coastal
(B9) aquifer system using variable density flow and solute transport
simulation, J. Hydrol., 92, 17-41, 1987.
1866 Voss AND SOUZA: SIMULATION OF NARROW TRANSITION ZONES

Visher, F. N., and J. F. Mink, Ground-water resourcesin southern unsaturatedfluid-density-dependentground-waterflow with energy
Oahu, Hawaii, U.S. Geol.Surv.Water SupplyPap., 1778, 133 pp., transport or chemically-reactive single-speciessolute transport,
1964.
U.S. Geol.Surv.Water Resour.Invest.Rep.,84-4369,409 pp., 1984.
Volker, R. E., Comparison of immiscible and miscible fluid models for
sea water intrusion in aquifers,paper presentedat Proceedings,7th W. R. Souza, U.S. Geological Survey, Box 50166, Honolulu, HI
Australian Hydraulics and Fluid Mechanics Conference, Inst. of 96850.
Eng., Brisbane,Qld., 1980. C. I. Voss, U.S. Geological Survey, 431 National Center, Reston,
Volker, R. E., and K. R. Rushton, An assessment of the importance of VA 22092.
some parametersfor seawaterin trusion in aquifersand a compari-
son of dispersive and sharp-interface modelling approaches, J. (ReceivedSeptember 10, 1986;
Hydrol., 56, 239-250, 1982. revised June 10, 1987;
Voss, C. I., SUTRA: A finite-element simulation model for saturated- acceptedJune 19, 1987.)