United States
Environmental Protection
Agency
Office of Research and
Development
\Afeshington DC 20460
EPA/600/R-99/110
January 2000
Analytic Element Modeling of
Coastal Aquifers
-2ft
/ i >,
/ i \
-80
TimeRSOOdays
2ft
40
80
ft
-2ft
-80
TimeF=4000days
40
|q|=0.002k
80
-------
EPA/600/R-99/110
January 2000
Analytic Element Modeling
of Coastal Aquifers
by
Mark Bakker
University of Minnesota, Minneapolis
Stephen R. Kraemer
U.S. Environmental Protection Agency
Ada, Oklahoma
Willem J. de Lange
Institute for Inland Water Management and Waste Water Treatment (RIZA)
Lelystad, The Netherlands
Otto D.L. Strack
University of Minnesota, Minneapolis
CR-823718
Project Officer
Stephen R. Kraemer
Subsurface Protection and Remediation Division
National Risk Management Research Laboratory
Ada, OK 74820
National Risk Management Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
Cincinnati, OH 45268
-------
Notice
The U.S. Environmental Protection Agency through its Office of Research and Development partially
funded and collaborated in the research described here under Cooperative Agreement No. CR-823718 to the
University of Minnesota, Dr. Otto Strack, Principal Investigator. It has been subjected to the Agency's peer
and administrative review and has been approved for publication as an EPA document. Mention of trade
names or commercial products does not constitute endorsement or recommendation for use.
All research projects making conclusions or recommendations based on environmentally related mea-
surements and funded by the Environmental Protection Agency are required to participate in the Agency
Quality Assurance Plan. This project did not involve physical measurements and, as such, did not require
a QA plan.
-------
Abstract
Four topics were studied concerning the modeling of ground-water flow in coastal aquifers with analytic
elements: (1) practical experience was obtained by constructing a ground-water model of the shallow aquifers
below the Delmarva Peninsula USA using the commercial program MVAEM; (2) a significant increase in
performance was obtained by implementing the theory for variable density flow in a computer program that
ran on a supercomputer using vectorization; (3) a new representation for the density variation was developed
that can simulate the change from brackish to fresh water more accurately: and (4) it was shown that for a
specific example of a bell-shaped transition zone a Dupuit model gives accurate results unless the bell-shape
is too narrow compared to the thickness of the aquifer.
in
-------
Foreword
The U.S. Environmental Protection Agency is charged by Congress with protecting the Nation's land,
air and water resources. Under a mandate of national environmental laws, the Agency strives to formulate
and implement actions leading to a compatible balance between human activities and the ability of natural
systems to support and nurture life. To meet these mandates, EPA's research program is providing data, and
technical support for solving environmental problems today and building a science knowledge base necessary
to manage our ecological resources wisely, understand how pollutants affect our health, and prevent or reduce
environmental risks in the future.
The National Risk Management Research Laboratory is the Agency's center for investigation of techno-
logical and management approaches for reducing risks from threats to human health and the environment.
The focus of the Laboratory's research program is on methods for the prevention and control of pollution to
air, land, water, and subsurface resources; protection of water quality in public water systems; remediation
of contaminated sites and ground water; and prevention and control of indoor air pollution. The goal of
this research effort is to catalyze development and implementation of innovative, cost-effective environmen-
tal technologies; develop scientific and engineering information needed by EPA to support regulatory and
policy decisions, and provide technical support and information transfer to ensure effective implementation
of environmental regulations and strategies.
Salt water intrusion is a potential threat to drinking water supplies in the coastal areas of the USA
due to over-pumping. In addition to the pumping, groundwater flow in coastal aquifers is affected by
the difference in density between fresh and salt water. Computer models provide a tool for predicting
the movement of salt water under past, present, and future pumping conditions. However, the simulation of
three-dimensional variable density flow is computationally expensive. This project investigates innovations in
algorithm formulation and computing architecture for problem solving involving variable density groundwater
flow, with the aquifers beneath the Delrnarva Peninsula, USA, providing context.
Clinton W. Hall, Director
Subsurface Protection and Remediation Division
National Risk Management Research Laboratory
IV
-------
Contents
Introduction 1
Chapter 1. An Analytic Element Model of the Upper Chesapeake Aquifer, Delmarva Peninsula (USA) .. 5
Chapter 2. A Dupuit Formulation for Variable Density Flow 17
Chapter 3. Implementation on the Supercomputer 26
Chapter 4. The Dupuit Approximation for Variable Density Flow in Coastal Aquifers 34
Chapter 5. The Specific Discharge Vector 47
Chapter 6. Results and Conclusions 54
References 57
Appendix Al 61
Appendix A2 68
Appendix B 76
-------
Acknowledgements
Funding for this project came through the EPA High Performance Computing and Communication pro-
gram, which is managed by Ms. Joan Novak, EPA Atmospheric Modeling Division, Research Triangle Park.
North Carolina. Computer time on the CrayC916 was obtained through a grant from the Minnesota Super-
computer Institute to the University of Minnesota. Dr. Erik Anderson and Dr. David Steward contributed to
the project while at the Department of Civil Engineering of the University of Minnesota, Minneapolis. The
technical content of the report benefited from reviews by Dr. Henk Haitjema. Indiana University, Bloom-
ington, Dr. Willem Zaadnoordijk, I WACO and Delft University of Technology, The Netherlands, and Dr.
Gualbert Oude Essink, Utrecht University, The Netherlands. Ms. Joan Elliott provided the editorial review.
Contact information (as of January 2000):
Dr. Mark Bakker, University of Nebraska, Dept. of Civil Engineering, 6001 Dodge Street, Omaha, NE 68182,
nibakkerCciunomaha.edu
Dr. Stephen Kraemer, USEPA. 960 College Station Road, Athens, GA 30605, kraemer.stephen@epa.gov
Dr. Wirn de Lange, RIZA, P.O. Box 17,8200 AA Lelystad, The Netherlands, W.dLange@riza.rws.rninvenw.nl
Dr. Otto Strack, University of Minnesota, 122 Civil Engineering, 500 Pillsbury Drive S.E., Minneapolis, MN
55455, strac001@tc.urnn.edu
VI
-------
Introduction
Motivation
Salt water intrusion is a potential threat to drinking water supplies in the coastal areas of the USA. Cities
along the coast are increasing their pumping of groundwater to support a rising population. The increased
pumping may result in an increase of chlorides in the well water due to the upconing of brackish groundwater.
Even a small concentration of chlorides will give the water a salty taste; water tastes salty to most people
if the concentration of chlorides is 0.25 g/1 or greater. Sea water, on the other hand, contains about 18
grams of chlorides per liter. The maximum guideline concentration set by the World Health Organization
is 0.25 grams of chlorides per liter. The potential upconing of brackish groundwater ma}' be studied by the
simulation of groundwater flow with a numerical model.
Groundwater flow in coastal aquifers is affected by the difference in density between frcs.li and salt water.
The fresh water is separated from the salt water through a brackish transition zone, in which the salinity
(and thus the density) of the water varies from that of salt water to that of fresh water. If the transition zone
is relatively thin, the transition from fresh to salt water may be modeled as an abrupt one: the fresh water is
separated from the salt water by an interface. If this is not the case, the effect of the variation in density on
the flow in the transition zone must be taken into account. The modeling of the flow generated by variations
in density, the variable density flow7, is the subject of this report. Specifically, it is investigated whether
groundwater flow in coastal aquifers can be modeled under the Dupuit approximation in combination with
analytic elements.
Background
Numerous numerical models are available to simulate variable density flow for both two-dimensional flow
in the vertical plane and three dimensional flow. The flow field may be modeled with the finite element
method or the finite-difference method while the solute transport equation ma}' also be solved by the random
walk method or the method of characteristics. Two-dimensional models include SUTRA (Voss. 1.984), and
MOCDENSE (Sanford and Konikow, 1985), the variable density version of MOC (Konikow and Bredehoeft,
1978). Three-dimensional models include H.ST3I.) (Kipp, 1986) and SWJ.CHA (Lester, 1991). Maas and
Emke (1988) developed a procedure to simulate variable density flow with numerical models for single
I
-------
density flow. Olsthoorn (1996) presented a method to use MODFLOW (McDonald and Harbaugh, 1984) for
variable density flow. The three-dimensional version of MOC, called MOC3D (1996), has been modified to
include density driven flow by Oude Essink (1998).
The construction of numerical models of three-dimensional variable density flow is limited by the avail-
ability of density data and by the speed of digital computers (Oude Essink and Boekelrnan, 1996). The
Dupuit approximation may be adopted to reduce the computation time and for the usual reason of simplic-
ity (Strack, 1995). The resistance to flow in the vertical direction is neglected in Dupuit models, and vertical
flow is governed by continuity of flow only; this leads to a hydrostatic pressure distribution in the vertical
direction. The Dupuit approximation is reasonable for flow in aquifers of great horizontal extend, also re-
ferred to as aquifers with shallow flow or regional flow. It is possible to construct regional groundwater flow
models of coastal aquifers by adopting the Dupuit approximation. The Dupuit theory for variable density
flow may be combined with any method to model Dupuit flow in a single density model. In this report the
flow field is modeled with analytic elements (Strack, 1989; Ilaitjema, 1995). Analytic element models have
been constructed successfully to simulate the fresh water head in the coastal aquifers of The Netherlands
(e.g., Minnerna and van der Meij, 1997). A model of a coastal aquifer beneath the Delmarva Peninsula is
described in this report.
The Dupuit theory for variable density flow
The Dupuit theory for variable density flow, as formulated by Strack (1995). assumes that the density distri-
bution in the aquifer is known at some time. In practice, the density is known only at a number of isolated
points. Strack (1995) proposes to represent the density distribution with a three-dimensional interpolator
function that interpolates between the points of known density; he suggests to use the multiquadric radial
basis interpolator (Hardy, 1971) for this purpose. It is noted that the multiquadric interpolator is a form of
Kriging; the interpolator is identical to Kriging with a linear variogram if the shape factor in the interpolator
is set equal to zero.
The density distribution (and thus the flow field) will change over time as the salt moves with the
groundwater flow. The flow is incompressible in Strack's model and the flow field thus represents the flow
at a given time. The evolution of the density distribution in the aquifer may be simulated by numerical
integration through time. During each time step, the velocity field is fixed and the salt is moved with the
groundwater flow. At the end of a time step the velocity field corresponding to the new density distribution
is computed and the process is repeated. This procedure is also known as successive steady-state solutions;
transient solutions are obtained with successive steady-state solutions throughout this report. Second order
processes that affect the salinity distribution, such as diffusion, are neglected. It may be expected that such
a procedure gives reasonable results for relatively short times (on the order of 25 years, as is of interest for
-------
most engineering problems).
Straok's theory has been implemented, prior to this study, in the proprietary computer program Multi-
layer Variable density Analytic Element Model (MVAEM). The computation time involved in the construc-
tion of large regional models of the fresh water head with MVAEM is significant. Multi-layer models that
include thousands of points where the salinity is specified take on the order of hours to compute a solution
on a high end PC (in 1996). A large portion of the computation time is used to compute the effect of
the variation in density on the flow. The computation times of the transient simulations involve a repeated
computation of the solution at different times plus a large number of evaluations of the velocity in the aquifer
and is of the order of days or more, as compared to the hours it takes to obtain one steady state solution.
These computation times diminish the practicality of the modeling approach.
The use of a Dupuit model to simulate variable density flow raises a number of additional issues. Transient
simulations require an accurate representation of the velocity field (not just of heads) and the numerical
integration requires an analysis of stability arid convergence. Furthermore, it must be investigated how
accurate a Dupuit model is to simulate the change of a salt distribution over time, especially in the case
of upconing of salt, or brackish water. Dupuit models are general!}' used for regional (shallow) flow and
may become inaccurate when the shallow flow assumption is not appropriate, for example near a partially
penetrating well where the flow field changes rapidly over a distance of several times the aquifer thickness.
And finally, the assessment of the upconing of salt water below a pumping well needs analysis of the physical
stability of the transition zone itself.
Objective
The objective of this report is to investigate the performance of the Dupuit theory for variable density flow
combined with analytic elements to model groundwater flow in coastal aquifers. It is impossible to answer
all the questions raised in the foregoing in one project. This project concerns four areas of study:
1. The analytic element modeling of groundwater flow in the first confined aquifer beneath the Delmarva
Peninsula.
2. The reduction of computation time by the use of a supercomputer.
3. The accurate representation of the density distribution.
4. The implications of adopting the Dupuit approximation for variable density flow.
-------
Report structure
The report is structured as follows. In Chapter 1, a model is presented of the fresh water head in the first
confined aquifer on the Del mar va Peninsula using the program MVAEM. The salt distribution (arid thus the
density distribution) was modeled separately prior to the groundwater flow simulations; use was made of a
three dimensional visualization package.
The Dupuit theory for variable density flow is derived following the paper by Strack (1995) in Chapter 2.
The density distribution is represented by a three dimensional multiquadric interpolator, as is done in the
program MVAEM, which was used for the modeling study. The second half of this chapter (starting with
the section "Head and potential") is highly mathematical and intended for the reader who is interested in
implementing the theory in a computer program. Reading this section is not required for understanding the
subsequent chapters.
It is investigated in Chapter 3 whether the use of a supercomputer will reduce the unpractically long
computation times to such a level that interactive modeling becomes possible. The Dupuit theory for variable
density flow was implemented in the analytic element code Variable Density Single Layer Wells Line-sinks
(VDSLWL), which is based on the public domain program SLWL (Strack. 1989). VDSLWL was written to
run on a vector machine, the CrayC916 at the Minnesota Supercomputer Institute. A brief description of
the density module is provided and some implementation issues are addressed.
The performance of the multiquadric interpolator to represent the density distribution is investigated
in Chapter 4. The multiquadric interpolator includes a shape factor that controls the smoothness of the
interpolator. In practice, this shape factor is often set close to zero to obtain a reasonable representation
of the density distribution. This leads to an acceptable approximation of the fresh water head, but the
resulting velocity field appears to be physically unrealistic. A new representation for the density distribution
is proposed to overcome this problem. This representation is better controlled but also less flexible. A
new exact solution for variable density flow in a vertical cross section is presented and compared to the
Dupuit solution. The derivation of the expressions for the specific discharge vector corresponding to the
new representation of the density distribution is length}' and is presented separately in Chapter 5. Finally,
results are summarized and conclusions drawn in Chapter 6.
-------
CHAPTER 1
An Analytic Element Model of the Upper Chesapeake Aquifer,
Delmarva Peninsula (USA)
Introduction
Salt water intrusion is a potential threat to drinking water supplies relying on groundwater in coastal aquifers
of the USA. The World Health Organization has set a guideline concentration of chlorides at 250 rrig/1; water
tastes salty to most people when the concentration of chlorides is greater than this value. The city of Ocean
City. Maryland, on the mid Atlantic coastal plain, has experienced a rise in the chlorides in one of its wells
from 70 rng/1 in 1975 to 215 mg/1 in 1988. The well field is probably experiencing upconing of brackish
water from the underlying aquifer due to increased pumping (Achmad and Wilson, 1993). The increased
pumping is needed to supply water to the growing population in the region.
The analytic element method for modeling groundwater flow has been extended to represent the influence
of a variation in density of the water (due to a variation in salinity) on the groundwater flow (Strack, 1995).
In this chapter we explore the application of analytic element modeling within a coastal aquifer where fresh
and sea water meet. The shallow fresh water aquifer systems of the Delmarva (Delaware-Maryland-Virginia)
Peninsula are stratified and multi-layered with alternating sand and clay layers, and are wedge-shaped.
thickening to the east and subcropping or pinching-out in the west (Vroblesky and Fleck, 1991). The aquifers
are bounded below by bedrock; fresh water meets sea water in the lower aquifers, directly underneath the
Chesapeake Bay, and offshore beneath the Atlantic Ocean.
The upper Chesapeake aquifer is considered a single geohydrologic unit at the regional scale, and is
bounded below by the St. Mary's confining unit, and above by the upper Chesapeake confining unit. The
upper Chesapeake aquifer subcrops into the surflcial Columbia aquifer in the western part of Delmarva, and
then gently dips to the southeast at a slope of about 0.01 (15 rti per 1600 in) (Vroblesky and Fleck, 1991)
(See Figures 1.1 and 1.2). The upper Chesapeake aquifer contains three major sand bodies of Miocene and
Pliocene age, which are, from lowermost to uppermost, the Manokin, Ocean City, and Pocomoke aquifers,
respectively.
-------
Rich
studyreg
Figure 1.1: The subcrop of the upper Chesapeake aquifer (shaded area) beneath the Delmarva Peninsula
USA
Delmarva
/ upper Chesapeake aquifer
vertical exaggeration 100X
Figure 1.2: Perspective view of the upper Chesapeake aquifer using GMS software (WES, 1997) (vertical
exaggeration lOOx)
-------
Approach
Groundwater flow in the upper Chesapeake aquifer was simulated using the Multi-layer Variable-density
Analytic Element Model (MVAEM Version 1.1 ©1995 Strack Consulting, North Oaks, MN). A description
of the mathematical basis of the point, line, and area elements used in MVAEM can be found in Strack
(1989), while the extension to include variable density flow is discussed in Chapter 2 of this report, and
in Strack (1995), Strack and Bakker (1995). MVAEM solves for the steady-state flow field, and uses the
Dupuit approximation; that is, resistance to vertical flow is assumed negligible. The range of applicability
of the Dupuit approximation for variable density flow is explored in Chapter 4.
MVAEM computes the influence of variable density flow, or density driven flow, using an estimate of the
continuous three-dimensional distribution of density in the aquifer system. The density must be specified at
a number of points in the aquifer (referred to as density points): MVAEM interpolates between these points
to obtain a continuous density distribution throughout the aquifer. These density points are inferred from
measurements of chloride concentration using an empirical relationship. For the upper Chesapeake aquifer
model, we used an empirical relationship between chloride concentrations and density (Van Dam, 1973).
Assuming the water temperature is 15 degrees Celsius, the density p may be written as
p = 1000 + 1.455[C7]/1000 - 0.0065(11 + 0.4[C7]/1000)2 (1.1)
where chloride concentration [Cl] is in mg/1 and the resulting density in kg/m3.
If is essential to evaluate flic interpolated density distribution and make sure it is a reasonable represen-
tation of what is believed to be the density distribution in the aquifer. If the interpolation of MVAEM does
not seem reasonable, additional points must be added where the density is specified to better control the
interpolator, for example below the ocean bottom.
The distribution of data points used in the model is shown in Figure 1.3. Many of the chloride data
points came from the QWDATA database of the USGS Water Resources Division in Baltimore, MD. Other
sources include Meisler (1989), Phelan (1987), Richardson (1992), and Woodruff (1969). Only fairly recent
observations (taken after 1940) were used assuming that the chloride (density) distribution did not change
significantly during this time period.
The continuous three-dimensional distribution of density was created using a rnultiquadric interpolator
in a procedure described in Chapter 2. It is noted that this interpolation technique is identical to Kriging
with a linear variogram. A series of chloride concentration values were added in order to obtain a more
realistic density distribution. Specifically, points were added in the upper Chesapeake aquifer beneath the
Chesapeake Bay and were given values similar to the lower surface waters, noting that resistance layers in
the Bay are partly absent. Data points below the upper Chesapeake aquifer were not used in the calculation.
because if is unlikely that the density distributions are related across the separating St. Mary's confining
-------
O estimated
O estimated
+ Meisler, 1989
O USGS QWGW database
A Phelan, 1987
0 Richardson, 1992
* Woodruff, 1969
densplot
Figure 1.3: Density points used by the interpolator
unit. Another series of points was added in the upper Chesapeake aquifer bounding the coastline of the
Atlantic Ocean based on the section model of Achmad and Wilson (1993). These points are shown as
"estimated" in Figure 1.3. A total number of 668 data points were used, 148 estimated, see Appendix Al.
A mesh of constant strength area elements was used to simulate the leakage between the surficial un-
confined aquifer and the upper Chesapeake aquifer. The leakages of the resistance elements are equal to
the difference between the given head above (in the surficial aquifer) and the head in the upper Chesapeake
aquifer, divided by the resistance of the upper Chesapeake confining layer. This condition is enforced at the
center of each element. The resistance values are based on estimates of the vertical hydraulic conductivity.
The resistance is equal to the thickness of the resistance layer divided by its vertical hydraulic conductivity,
and has the units of time. (The resistance is the inverse of the conductance, a parameter used in many other
models.) The layout of the elements and table of resistances and heads assigned can be found in Appendix
A2.
Inhomogeneity polygons were used to represent variable aquifer thickness, variable hydraulic conductivity,
and sloping base, as shown in Figure 1.4. Aquifer properties are constant within each polygon. These
polygons are composed of line doublet elements which create the appropriate jump in the discharge potential
and maintain continuity of flow. The base within each polygon is horizontal, and steps in the base were
limited to half the aquifer thickness. The spatial pattern of doublet elements is based on the transmissivity
distribution shown by Leahy and Martin (1993), and the estimations of the aquifer base on well logs reported
in Vroblesky and Fleck (1991). The two small doublet polygons represent the local increase in thickness and
transmissivity reported near Ocean City, Maryland.
-------
Figure 1.4: MVAEM doublet polygons representing heterogeneous upper Chesapeake aquifer (B, H, k are
aquifer base elevation above mean sea level (m), thickness (m) and hydraulic conductivity(m/day)
-------
4400000.00
4350000.00
4150000.00
Piezo
'"OOOOtJsttoO.OO 400000.00 450000.00 500000.00 550000.00
600000.00
Figure 1.5: Locations of observations wells lor fresh water heads
MVAEM was run lor two scenarios: (1) variable density - fresh and salt water; and (2) single density -
all fresh water. The predictions ol fresh water heads were compared to monitoring wells located in Sussex
County, Delaware, and Wicomico and Worchester Counties in Maryland (see Figure 1.5). The fresh water
head is defined as the elevation to which water rises in a standpipe il the standpipe is filled with fresh water
only. The selected observation wells are screened in the upper Chesapeake aquifer, and do not report an
influence ol nearby pumping wells.
Results
The upper Chesapeake model representing variable density groundwater flow beneath the Delmarva Penin-
sula has not been extensively calibrated, and results are considered preliminary. Figure 1.6 shows a fence
diagram ol the density distribution generated with the multiquadric interpolator; notice that the density
distribution is continuous in three-dimensions. The transition from fresh to sea water can be seen clearly
along the Delmarva shoreline and coastline in the contour map ol Figure 1.7. A contour map ol the MVAEM
fresh water heads (at elevation 50 m below sea level) is shown in Figure 1.8.
The model predicted fresh water heads were compared to observed water levels in nine wells screened in
the upper Chesapeake aquifer (See Figure 1.5). The observed water levels are based on a 5 year average over
the period 1987-1992 (James et al., 1992) and are based on the density ol water at the well screen. It is noted
10
-------
1010 1030
1000 1020
Density (kg/m^)
denxyz
Figure 1.6: Fence diagram of three-dimensional density distribution (kg/m3) inferred from discrete measure-
ments of chloride concentrations and the multiquadric interpolator
elevation z = -50 m
Figure 1.7: Contours of water density at elevation z=-50 m
11
-------
(5.2e5.4.3e6)
elevation z = -50 m
headsjw
(3.7c5,4.lrf)
Figure 1.8: Contours of MVAEM fresh water heads at elevation z=-50 m
12
-------
Table 1.1: Comparison of observed heads with simulated heads. Unit of heads is meters.
Well ID
Nf44-01
M52-11
Oh54-01
Oi 24-06
Pf24-03
Qh54-04
Rj 22-05
WOAe23
WODe-36
UTM-x
468900
487572
484053
490558
468686
484210
495517
474254
474211
UTM-y
4292706
4290620
4280763
4285068
4275014
4262639
4257631
4254373
4233291
z (m msl)
-30
-40.8
-81
-70.0
-37.5
-90.8
-120
-71.6
-89.9
obs. head
8.9
1.8
2.7
1.8
13.4
4.5
1.0
8.5
5.2
Variable density
fw head
9.17
2.71
3.01
0.032
13.8
6.38
0.61
12.7
3.51
difference
0.27
0.91
0.31
-1.8
0.42
1.9
-0.39
4.3
-1.7
Single density
fw head
9.17
2.69
3.03
0.0083
13.8
6.39
0.56
12.7
3.53
difference
0.27
0.90
0.34
-1.8
0.43
1.9
-0.43
4.3
-1.7
that, due to the small number of observations, no analysis has been performed to determine whether these
values represent regional flow conditions. The difference between the MVAEM model predicted fresh water
head and the observed head, for both the variable density and single density (all fresh water) simulations.
are shown in Table 1.1.
MVAEM predicts the water exchange between the unconfined surficial aquifer and the upper Chesapeake
aquifer. The leakage (m/d) for the area elements is shown in Figure 1.9. The mean value of the leakage is
8.79E-5 m/d (1.3 in/yr), while the maxirnurn leakage entering the aquifer is 2.87E-3 m/d (41.2 in/yr) and
the maximum leakage leaving the aquifer is 3.89E-4 m/d (5.6 in/yr). Figure 1.9 shows most of the leakage
to the upper Chesapeake aquifer occurs in the uplands of Delrnarva. while most of the exchange with the
surficial aquifer occurs along the shore. The USGS reports an average flux of 122 ft3/s through the Upper
Chesapeake aquifer layer (4690 mi2) of their model giving a leakage of 2.457E-5 m/d (0.35 in/yr) (Fleck and
Vroblesky, 1996).
Discussion
The definition of the continuous three-dimensional chloride concentration (density) distribution in the up-
per Chesapeake aquifer is problematic given the paucity of available data. Most of the observation wells
reporting a chloride concentration were of the fresh water beneath the Delrnarva Peninsula. Only a handful
of observation wells were available with a "salty" chloride concentration, and only a couple of these wells
reported a chloride concentration that varied with depth. The resulting density contours are essentially
13
-------
Logleak
leakage - in (m/d)
I I icr3to icr6
I I < icr6
leakage - out (m/d)
I I icr3to icr6
< 10
-6
Figure 1.9: Distribution of leakage into and out of the upper Chesapeake aquifer
14
-------
constant with depth, as shown in Figure 1.6. It is unlikely that this distribution is realistic, especially in
the few areas where the density decreases with depth. Also, given the sensitivity of velocity calculations
to the density variation, and given the rather questionable predicted density distribution based on so few
data points, no calculations of the velocity distribution in the aquifer are presented here. Therefore, no
predictions of the movement of the salt water transition zone are offered.
The MVAEM model did a reasonable job in predicting the fresh water heads in the aquifer, based on the
comparison to the observed heads in monitoring wells, and the contours of predicted heads. It should be
noted that MVAEM is a steady-state model, and the comparison was made to average heads. The heads are
known to vary up to 1 m over the seasons. Also, the modeled heads were constrained by the head specified
area elements, which were based on published water table contour maps. The solutions for fresh water heads
are also relatively insensitive to the density distribution, at least at the observation points, as evidenced in
Table 1.1. Further investigation of the reasonableness of the MVAEM fluxes is warranted.
Conclusions
The analytic clement method was used to build a groundwater flow model of the upper Chesapeake aquifer
of the Delmarva Peninsula, USA. The application of the MVAEM code demonstrated the potential of the
analytic element method for representing variable density flow and to increase the understanding of the
transition zone between fresh and sea water. No definitive site specific conclusions arc offered given the
uncertainties involved in this particular model application.
The analytic element method has the potential to simulate a reasonable representation of the fresh water
head distribution, given adequate investment in model calibration. However, analytic element models, as
with other numerical models, are only as good as the input data and adequacy of the conceptual model. In
addition, the appropriateness of using the analytic element method to represent the aquifer base elevation
in a piece-wise manner in order to approximate a smoothly sloping base elevation was not examined. More
investment is needed to better define the variation of density in space. Better definition of aquifer geometry
and properties is needed. More monitoring wells (as opposed to water supply wells) are needed to better
define the fresh water head distribution. Measurements of fluxes to tidal rivers and shorelines would be
useful. These observations would facilitate a complete water balance analysis.
Additional complexities of the aquifer system may be represented in the future with new developments of
the analytic element method. It is possible to build a multi-layer model and better represent the influences
of the surficial aquifer. For example, the major streams could be represented as curvilinear line elements,
while the minor streams may be lumped into an effective feeding resistance based upon their drainage density
(De Lange, 1996). Alternatively, leakage to the upper Chesapeake aquifer from the surficial aquifer could
be represented using advanced variable strength area elements (Strack arid Jankovic, 2000). Also, future
15
-------
developments in the analytic element method include functions to represent a continuously sloping aquifer
base and a transient aquifer response.
Until these challenges are met, it is too early to assess fully the practical advantages, and disadvantages,
of the analytic element method in comparison to numerical techniques such as finite differences and finite
elements.
16
-------
CHAPTER 2
A Dupuit Formulation for Variable Density Flow
Introduction
The Dupuit theory for variable density flow was presented by Strack (1995). The theory is reproduced here,
in a slightly different form, for the case of confined flow in a piecewise horizontal aquifer. A Cartesian x\,
2>2, z coordinate system is adopted with the z axis pointing vertically upward.
The specific discharge vector field for variable density flow is rotational, even if the aquifer is homogeneous
and isotropic. Strack (1995) showed, however, that the discharge vector field (the specific discharge integrated
over the saturated thickness of the aquifer) is irrotational. Thus, a potential 4> may be defined such that
the discharge vector is minus the gradient of this potential. As such, the three-dimensional rotational flow
problem is reduced to a two dimensional potential flow problem. Furthermore. Strack (1995) approximated
the pressure distribution in the vertical direction as hydrostatic (the Dupuit approximation) to obtain a
relation between the potential and the fresh water head. The fresh water head is defined as the elevation
to which water rises in a standpipe if the standpipe is filled with fresh water only (e.g., Lusczynski, 1961).
Basic Equations
Darcy's law in terms of pressure is
K dp
fj, dxi
K dp
fj, dz
where qi, q-2, qz are the components of the specific discharge vector in the x\, x-2, z directions, respectively,
i is the intrinsic permeability, /z is the dynamic viscosity, p is the pressure, p is the density of water and g
is the acceleration due to gravity. The fresh water head,
-------
where pj is the density of fresh water and Z is the elevation above an arbitrary datum. Combination of
(2.1) and (2.2) gives
dd>
'dxi (2_3)
where k is the hydraulic conductivity of fresh water
and // is the dimensionless density
(2.4)
Pf
The discharge vector is the total flow integrated over the saturated thickness of the aquifer and has
components
where Zf, and zt are the bottom and top of the aquifer, respectively. Integration and differentiation may be
reversed if zt, and zt are not a function of x\ and XT (as for a confined aquifer with horizontal base and top)
k / dz
where the potential $ is defined as
Continuity of flow gives
where Nt is the water leaving the aquifer at the top and Nt, is the water entering the aquifer at the bottom;
both Nt and N;, may be functions of x\ and x^. The partial derivative in the x,: direction is written as di
and the Einstein summation convention is adopted for repeated indices; only the index i is used to indicate
the components of a vector and summation is implied for the horizontal directions only. Substitution of (2.7)
for Qi in (2.9) gives
V2* = Nt - Nb
where V" is the Laplacian in the horizontal directions.
18
(2.10)
-------
The Dupuit approximation
The Dupuit approximation is adopted, which means that the resistance to flow in the vertical direction is
neglected and the pressure distribution is hydrostatic, so that dp/dz —pg and thus
Integration gi ves
where Ffai.x?) is an, as of yet. unknown function of x\ and XT. Substitution of (2.12) for 6 in (2.8) and
division by k gives
tidz
(2.13)
Performing the latter integration leads to an expression for F
ff
(2.14)
Substitution of (2.14) for F into (2.12) gives an expression for the head as a function of the potential
(2.15)
or vice versa
k.H(j>+ kl-1 / vdz ~~k vdzdz
Expressions (2.15) and (2.16) are identical to expressions (40) and (39), respectively, in Strack (1995).
The specific discharge vector
The horizontal components of the specific discharge vector are obtained from differentiation of (2.15)
= 1,2
(2.17)
where integration and differentiation are interchanged for the first integral in the last term. The vertical
component of flow may be obtained from continuity
/-)„
(2.18)
19
-------
which gives
The divergence of q,. may be obtained from (2.17) which gives for (2.19)
The first and last terms on the right-hand side of (2.20) may be combined, using (2.9)
II
For the case that Art = Nf, = 0 integration of (2.20) gives
This equation may be simplified by interchanging differentiation and integration and by rearranging terms
H
H
H J j ' H .,
Combination of (2.20) through (2.23) gives the general expression for qz
It may be verified from equation (2.24) that qz(z = zt) = Nt and qz(z = z;,) = Ni,, as asserted. The
expressions for the specific discharge vector. (2.17) and (2.24), are the same as equations (42) and (50) in
Strack (1995).
The three-dimensional multiquadric interpolator
The dimensionless density distribution ;/, see (2.5), is represented by a multiquadric interpolator (Hardy.
1971). which is written as follows
= Es~
0
- ct
-------
where
and S3 is the horizontal scale factor. The M + 1 constants a (m = ()...., M) are determined from M + 1
conditions. M conditions are obtained by requiring that v equals a specified value v at M collocation points
n m. m,
Z_v
m=l
'in
The constants A (m = 1,..., M) may be chosen arbitrarily, but do afreet the shape of the interpolator
function. The smaller the value of A, the sharper the change of the interpolator function at a collocation
point. The multiquadric interpolator is a convenient interpolator for the density distribution, especially when
m
the constants A are chosen small (preferably zero) relative to the size of the model domain (VanCervcri and
Maas, 1994).
Head and Potential
Expressions for the head in terms of the potential and vice versa (equations (2.15) and (2.16)) may be
m
obtained by integration of the dimensionless density distribution (2.25). The variable 6 is introduced as
1
where
and
These integrals may be checked by differentiation. Alter application of the limits of the double integral, q>
becomes
1 I
—
-------
where
/ -.or- w' - f, / ''ff!- \ m / W>>\ .", T'o
rmt = V P i(xi - a;i) + (x'2 - x^i] + (zt - z)1 + A2
The expression for the head in terms of the potential becomes
(2.33)
(2.34)
and the expression for the potential in terms of the head is
M
= kH4, - kH a + kHa z - ^
2H
The specific discharge vector
The specific discharge vector depends on the first and second derivatives of the density distribution (see
equations (2.17) and (2.24)) which may be written as
*~—^
m 0
m dr
OJ-T— i = 1.2
Differentiation of r (2.26) is straightforward and gives:
Expressions for the components of the specific discharge vector may be obtained from equations (2.17)
through (2.24) by working out the integrals. The following two integrals are used (these can again be verified
by differentiationl
dr
Til Til
-------
The contribution of multiquadric point m, with unit strength, to qi will be called qi and is obtained by
combination of (2.1), (2.23), and (2.24)
H
Using that
and rearranging terms gives
H
H
H
' rmt
rmt
(2.44)
(2.45)
An expression for the vertical component of the specific discharge vector may be obtained if the double
integral of the Laplacian of the density distribution is carried out. The following integrals are used
.9 . , 77),,
so that
r r
i i T~?2rn 7 7 o O'2 If '""\ 1 / t!!f : t!!f\ t!!f
I I V rdzdz = 2p \(z - z) m(z — z | r) - r
which may be simplified, after some algebra, to
(2.49)
The vertical component of the specific discharge vector may be obtained by substitution of (2.49) for the
double integrals in (2.24) and gives, after a considerable amount of algebra,
rn rn
H
H
-------
It is noted that the specific form of (2.50) is chosen so that the logarithms are the same in the expressions
for q\ and q z\ this will facilitate computations.
The three components of the discharge vector may now be obtained as follows
M
H
H
Rotation
As stated before, the specific discharge field is rotational. It will be shown that, although continuity of flow
is met exactly when making the Dupuit approximation (see (2.18)), the curl of the specific discharge vector
is represented approximately. Darcy's law may be rewritten as
% -diX i 1,2
_.
dz
where
The curl R of the specific discharge vector may be written as
R = (d'2qz — dzq-2- dzqi — d-\qz, diq% — fhqi)
where dz stands for partial differentiation in the ^-direction. Differentiation of (2.52) and substitution of
the result in (2.54) gives
R = (-kd-1v,kd]_v,Q) (2.55)
where it is used that % is single valued (d^d^x = d'l^iX]-
The curl of the specific discharge vector obtained with the Dupuit approximation may be computed by
differentiation of (2.17) and (2.24). Differentiation gives
dzqi = kdiv dzq-2 = kd^v diq\ = diq-2 (2.56)
so that the curl of the specific discharge obtained with the Dupuit approximation becomes
R = (-kd-2v + (hqz,kdvv + dvqz, 0) (2.57)
Equation (2.57) is only equal to (2.55) for the case that d^qz = d\qz = 0; this is the case if V2// = 0. For
almost all other cases, the curl is represented approximately.
-------
Implementation
The Dupuit formulation for variable density flow may be implemented in any groundwater code for Dupuit
flow of a single density fluid. Prior to this study, this theory has been implemented, in a slightly different
manner, in the commercial program Multi Layer Analytic Element Model (MLAEM; Strack, 1992) and
is called MVAEM (where the V stands for Variable density). The flow field in MLAEM and MVAEM is
modeled with analytic elements (Strack, 1989; Ilaitjema 1995). The implementation of the theory in the
analytic element code Single Layer Wells Line sinks (SLWL; Strack, 1989) is discussed in the next chapter.
-------
CHAPTER 3
Implementation on the Supercomputer
Introduction
The analytic element code Single Layer Wells Line-sinks, SLWL (Strack, 1989), is modified to run on a
CrayC916. The Dupuit theory for variable density flow is implemented and the resulting program is called
Variable Density SLWL (VDSLWL); VDSLWL is an experimental code.
Implementation
The flow field consists of a potential flow part plus a part due to the variation in density. This may be seen,
for example, from the equations for the horizontal components of the specific discharge vector (2.51):
M
^i . \~^ mm -10 fo 1\
qi=—+^iaqi « = 1, 2 (3.1)
rr i ..... 1
where a is the strength of multiquadric point m and qi (m = 1. .... M) depend on the density distribution
and aquifer parameters only (sec 2.45). It is noted here that this docs not mean that the flow field written in
the form (3.1) is the sum of single density flow plus variable density flow. The discharge vector Qi depends
on the boundary conditions, which in turn depend on the density distribution.
The flow field is modeled with analytic elements. The analytic elements used here are wells, constant
strength line ..... sinks, and constant strength, circular ponds. Each element may be written as the product of
a strength parameter (for example, the discharge of a well) and an influence function that depends on the
geometry only. The influence functions for the analytic elements used in this study may be found in Strack
n n
(1989). The strength of element n is called s and the influence function A(xi,X2) so that the potential due
to N elements is
The derivative of A in the Xi direction is written as A.; so that the discharge vector may be written as
N
-------
and the specific discharge vector becomes
1 N r, M
I r-—^ 11 " ^-—^ mrn
-------
COMPLEX FUNCTION COMLS(CZ,CZS,CZE)
IMPLICIT COMPLEX (C), LOGICAL (L)
DATA RPI /3.1415926/
CBZ=(2.0*CZ-(CZS+CZE))/(CZE-CZS)
COM1=(.0,.0)
COM2=(.0,.0)
IF(CABS(CBZ+1.).GT..0001) COM1=(CBZ+1.)*CLOG(CBZ+1.)
IF(CABS(CBZ-1.).GT..0001) COM2=(CBZ-1.)*CLOG(CBZ-1.)
COMLS=COMl-COM2+2.*CLOG(.5*(CZE-CZS))-2.
COMLS=.25/RPI*CABS(CZE-CZS)*COMLS
RETURN
END
The vectorization on the C916 is automatic, provided that the code is in a vectorizable form. Only inner
do loops can be vectorized and the do loop cannot contain calls to other subroutines or external functions.
In addition, a recurrence relation, like the line
CFLSU=CFLSU+RLSDIS(IAD)*COMLS(CZ,CLSZS(IAD),CLSZE(IAD))
in CFLSU, is not guaranteed to be vectorized correctly.
The do loop in CFLSU is not vectorizable because it has a call to an external function COMLS. Furthermore,
problems are anticipated because of the presence of a recurrence relation in the do-loop. The do-loop will
be vectorizable if the function COMLS is brought inline with the function CFLSU.
The program fpp (FORTRAN pro processor) on the C916 may be used to bring COMLS inline with CFLSU.
This results in a do-loop that is indeed vectorizable, but does not modify the recurrence relation at the end
of the do-loop. As such, the code does not give correct results.
The code was rewritten to bring COMLS inline with CFLSU, instead of using Jpp. The recurrence relation
at the end of the do loop was moved to a separate do loop. Compiler directives were added to obstruct cf'77
from vectorizing the latter do loop. The modified code is
COMPLEX FUNCTION CFLSU(CZ)
IMPLICIT COMPLEX (C), LOGICAL (L)
INCLUDE 'SLLS.CMN'
DIMENSION CSCR(IOO)
DATA RPI /3.1415926/
28
-------
Table 3.1: Bench .mark results
Machine
486PC.50MHz
C916 w/o vectorization
C916 w/ vectorization
Wallclock (seconds)
130
30.25
2.65
CPU time (seconds)
-
22.22
2.49
CFLSU=(.0,.0)
DO 100 ILSF=1,NLSF
IAD=ILSPTF(ILSF)
CZS=CLSZS(IAD)
CZE=CLSZE(IAD)
CBZ=(2.0*CZ-(CZS+CZE))/(CZE-CZS)
COM1=(.0,.0)
COM2=(.0,.0)
IF(CABS(CBZ+1.).GT..0001) COM1=(CBZ+1.)*CLOG(CBZ+1.)
IF(CABS(CBZ-1.).GT..0001) COM2=(CBZ-1.)*CLOG(CBZ-1.)
COM=COMl-COM2+2.*CLOG(.5*(CZE-CZS))-2.
COM=.25/RPI*CABS(CZE-CZS)*COM
CSCR(ILSF)=RLSDIS(IAD)*COM
100 CONTINUE
CDIR$ NOVECTOR
DO I=1,NLSF
CFLSU=CFLSU+CSCR(I)
ENDDO
CDIR$ VECTOR
RETURN
END
The vectorized code ran at 293.95 Mflops. A performance of over 300 Mflops may be obtained by replacing
the second do loop in CFLSU by a library call that sums up the CSCR array. Results of a benchmark for the
outlined problem are shown in Table 3.1.
It may be concluded from Table 3.1 that vectorization gives a significant increase in performance. To
enable vectorization. the code has to be rewritten to create do loops that do not call any external functions
-------
or subroutines and do not include recurrence relations. This influences the modular structure of the program.
The use of fpp to inline code has to be used with caution, because fpp does not consider the presence of
recurrence relations.
Performance
The change of the salinity distribution over time may be approximated using consecutive steady—slate ap-
proximations. Given the initial density distribution, points where the density is specified (referred to as
density points here) are moved with the flow over a certain time interval. The velocity field is fixed over the
time interval. As a first order approximation, the density at a point that moves with the flow is assumed
not to change; processes such as diffusion and dispersion are not taken into account. At the end of the time
interval, the new locations of the density points may be used to compute a new density distribution and thus
a new velocity field. This process is repeated until the desired time is reached.
The computational effort of transient simulations may be separated into two parts: (1) the computation
of the density distribution and (2) the advection of points during a time interval with a suitable particle
tracking technique. The computation of the density distribution consists of solving flic system of linear
equations presented by equations (2.27) and (2.28). Such a system may be solved using a standard routine
and is not used as a bench mark.
The advection of the density points requires the repeated evaluation of the three components of the
specific discharge vector given by (2.51). The parts of equations (2.51) that represent the influence of the
variation in density consist of simple sums, which are easily vectorizable. Computations have been performed
on a PC with an Intel Pentium Pro 200 MHz processor and a Cray C916 with 9 processors and a clockspeed
of 238 MHz. The computer program is written in FORTRAN??, using the Lahey F77L-EM/32 compiler on
the PC and the CrayCFT?? compiler on the C916.
A benchmark is performed for a hypothetical problem where the flow is caused by density differences
only (hence Qi = Q? = Nf, = Nt = 0). The components of the specific discharge vector were computed
ten times at each density point (representing a procedure that needs 10 evaluations of the velocity over
a time interval). The results are presented in Table 3.2. Column 1 is the number of points where the
density is specified, column 2 the computation time (in seconds) on the PC, Column 3 the computation
time on the CrayC916 without vectorization, and column 4 the computation time on the CrayC916 with
vectorization and autotasking, using I processor. The computational speed on the PC and the CrayC916
without vectorization are similar while the computation time on the CrayC916 with vectorization is an order
of magnitude smaller. The CrayC916 with vectorization runs at a speed of 363 Mflops for the case of 1600
density points.
It is of interest to note that the computation time to solve the system of equations is comparable to
-------
Table 3.2: Comparison of performance with time in seconds
Number of
points
200
400
600
800
1000
1200
1400
1600
3200
PC
2.31
9.01
20.38
35.87
55.97
80.63
109.74
143.52
CrayC916
w/o vectorization
2.39
9.22
20.48
35.95
CrayC916
w/ vectorization
0.39
1.20
2.37
3.95
5.88
8.32
11.08
14.43
54.92
the computation time of the evaluation of the discharge vector as discussed above. The system was solved
using LDLi-decomposition. This procedure is known to be inefficient and hard to vectorize. It may have an
advantage, however, if the procedure to simulate the change of the density is modified as follows. Instead of
following a density point with the flow, it is determined what density will arrive at the density point. Thus
the location of the density point is fixed, but the density will change. If the location remains the same,
so will the system of equations (2.27) and (2.28) that has to be solved; this system depends only on the
location of the points. Hence, the matrix equation has to be inverted only once and the LDU decomposition
stored. Back substitution using the LDU-decomposition consists of the two matrix multiplications and the
computation time involved is insignificant.
The evaluation of the specific discharge vector at the M density points is an M2 process (M evaluations at
M locations). Hence, a graph of the number of points versus the computation time should be a straight line
on log log paper; this graph is presented in Figure 3.1 for the PC and the CrayC916 with vectorization. The
slopes of the lines represent the real powers of flic process. The computations on the PC arc approximately
an M1-99 process and on the CrayC916 an M1-7"1 process.
The benchmark shows that vectorization results in an order of magnitude increase in speed. Use of a
vector machine makes it possible to solve problems with thousands of density points, as needed in regional
modeling, in a timely manner. Much larger systems of equations can be handled on the Cray than on the
PC.
-------
1000
100
10
1
0.1
X
Jl
100 1000
Number of density points
10000
Figure 3.1: Comparison of computation times; PC (triangles), CrayC916 (squares).
Instructions for use of the density module in VDSLWL
VDSLWL has a command line interface. Commands can be typed in or read from an input file. If VDSLWL
is run on the supercomputer, instructions are read from the file IN.SL, output is written to the file OUT.SL
(unless specified differently) and messages are written to the file MES.SL. The input of aquifer parameters
and analytic elements is described in Strack (1989). The only change that has been made is that when a
head is specified, either for a head specified well, a head specified line-sink or the reference point, both the
head and an elevation have to be specified. A horizontal grid of the density distribution may be obtained
with the command (N)(Z) where Z is the elevation in the aquifer where the grid is computed.
The check module includes three new commands: , (x,y,z), (x,y,z). These
commands will return input density information, specific discharge at a point and dimensionless density at a
point. The same convention is adopted for the density module. The command accesses the density
module. The following commands are available in the density module
(x.y.z.nu,delta)..(b)....
(tol)..
The first command specifies the data at a point. Specify the dimensionless density (v) not the density
(p). A separate A may be specified for each point. The factor (3 is the horizontal scale factor. All density
data should be entered consecutively. After all density data is entered solve the system of equations by
typing . The command provides data to check whether the solution is correct. The
32
-------
command finds all the points within a specified tolerance from each other. And finally,
returns command to the main menu. An example data file with the following three data points
x y z v A
1000 1000 -20 0.02 1
1000 2000 -30 0.01 1
2000 2000 -10 0.005 1
would look as follows
DENS
BETA .01
1000 1000 -20 0.02 1
1000 2000 -30 0.01 1
2000 2000 -10 0.005 1
SOLVE
CONTROL
QUIT
The CONTROL statement returns the following information
I,ALPHA,NUCOMPUTED,NUGIVEN 1 -6.52104E-04 2.00000E-02 2.00000E-02
I.ALPHA.NUCOMPUTED.NUGIVEN 2 2.57503E-04 l.OOOOOE-02 l.OOOOOE-02
I,ALPHA,NUCOMPUTED,NUGIVEN 3 3.94602E-04 5.00000E-03 5.00000E-03
NUO.SUM OF ALPHA-S 1.01553E-02 O.OOOOOE+00
Some notes on modification of the code
Details of the analytic clement part of the code are given in Strack (1989). The listing of the density module
of VDSLWL is presented in Appendix B. Two issues are of interest for VDSLWL.
1) In the top part of the file SLMN.FOR. the input and output are directed to either a file (for use on the
supercomputer) or the console (for use on a PC).
2) The maximum number of density points is specified in the common block vardens.cmn and in the routine
ludcmp (in the file ludcmp.for).
-------
CHAPTER 4
The Dupuit Approximation for Variable Density Flow
In Coastal Aquifers
Introduction
The implications of the Dupuit approximation for variable density flow in coastal aquifers are investigated.
The density distribution in the aquifer is represented by a number of surfaces of constant density. The ele-
vations of the surfaces are approximated with multiquadric interpolators; the density varies linearly between
them. A new exact solution is derived for two-dimensional flow in the vertical plane, and is compared to the
Dupuit solution. The problem used for comparison consists of a bell-shaped transition zone between fresh
and salt water (as may be expected from upconing under a pumping well).
The density distribution changes with time because the salt moves with the groundwater. The change
of the density distribution is simulated by computing the change of the surfaces of constant density through
time. A transient simulation is presented for a hypothetical problem.
Representation of the density distribution
In practice, the density distribution as a function of x\, x%, and z is unknown. Rather, the density is known
at a number of isolated points in the aquifer. The integrals in the expressions for the head, potential, and
specific discharge vector, derived in Chapter 2. may be carried out when a functional form is chosen to
represent the density distribution. In light of the expressions for the specific discharge vector, (2.17) and
(2.24), the function that represents the density distribution needs to have continuous first derivatives and
Laplacian in the two horizontal directions to ensure a continuous flow field.
Strack (1995) proposed to represent the density distribution using a three dimensional radial basis inter-
polator function; such functions have an infinite number of continuous derivatives (except for some special
cases). This approach has been implemented in the commercial software package MVAEM, and for this
project in the program VDSLWL, and has been used successfully to simulate the distribution of fresh water
heads in parts of The Netherlands (e.g., Van Gerven and Maas. 1994; Minnema and van der Meij, 1997).
Radial basis functions, such as the multiquadric interpolator (Hardy, 1971), are nicely behaved functions
34
-------
Function to be interpolated
Multi-quadrics interpolation
Figure 4.1: Overshoot and fluctuation of the interpolator function
that are suitable for the representation of continuously varying functions. However, radial basis functions,
as well as most other interpolation functions, are not suitable for the representation of a function that has
discontinuities in its derivatives; for example, a function that varies in one part of a domain and is constant
in another part. Consider, for example, a one-dimensional function that varies linearly for x < 0 and is
constant for x > 0 (the solid line in figure 4.1). When this function is approximated by a multiquadric
interpolator a problem arises near x = 0, because the interpolator cannot make a sharp bend (the derivative
of the interpolator is continuous). As a result, the interpolator will fluctuate around the function that it
represents over a considerable distance from the sharp bend (the dashed line in Figure 4.1). Such a behavior
is undesirable if the interpolator is to be used for the representation of the density distribution. Another
problem with the three-dimensional interpolator function is the shape factor A. In practice, this shape
factor is often set close to zero to obtain a reasonable representation of the density distribution (Van Gerven
and Maas, 1994). This results in a reasonable variation of the fresh water head, but the resulting velocity
field appears to be physically unrealistic. (It turns out that for A approaching zero, the Laplacian of the
interpolator function (V2^) tends to infinity, as will be explained later in this chapter.)
As an alternative functional representation, it is proposed to divide the aquifer up in a number, say N+l,
of regions (see Figure 4.2). The nih region is bounded below by surface n of constant density v = vn and
on top by surface n + 1 with density v = vn+i- The density in region n varies linearly from vn to vn+i in
the vertical direction. If the elevation of surface n is represented by the function Cn(zi, xi)> then the density
may be written as
+
z ~ >™ / \
-(Vn+l -Vn)
Cn+l — Cn
<
(4.1)
35
-------
v-u
Figure 4.2: Surfaces of constant density
The elevation of surface n, C,n(xi->x2)-> maY be represented by a two-dimensional interpolation function.
In practice, salinities are measured in nested observation wells, such that the density is known at a number
of elevations at one location. These measurements may be used to estimate the elevations of surfaces of
constant density; the interpolation function is fitted through these elevations.
The advantage of representation (4.1) is that the change from constant density (in the salt and fresh
water zones) to varying density (in the transition zone) is represented accurately. Also, it is convenient to
have explicit expressions for the elevations of surfaces of constant density while doing transient simulations,
as will become apparent in the final part of this chapter. A third advantage is that the integrations in the
expressions of the head, potential, and specific discharge vector may be carried out independent of the choice
of the functions £„, since they do not depend on z.
The new representation is, however, more restrictive, since the density is approximated as piecewise linear
in the vertical direction. This restriction may be overcome somewhat by approximating the vertical variation
by a second or third order polynomial, but that has not been explored. In addition, representation (4.1) will
have to be modified if surface n intersects the base (or top) of the aquifer. Such a modification is possible,
but falls outside the scope of this chapter.
The integrations and differentiations for the specific discharge vector are carried out in chapter 5. It is
noted that as a result of the choice of this particular representation, the vertical variation of qx is quadratic
in the transition zone and is constant in the fresh and salt water zones; the vertical variation of qz is linear in
the constant density areas and is a third order polynomial in the transition zone (where the vertical variation
of v is linear).
36
-------
Comparison with exact solution
Strack arid Bakker (1995) showed that solutions obtained with the Dupuit approximation for variable density
flow compare well with exact solutions for problems where the density varies in the x\ direction only. In this
section, a comparison is made to problems where the density varies in xi and z directions. Flow is considered
in the vertical x,z plane (the index 1 is dropped from xi for notational convenience in this section). The
problem used for comparison is a transition zone that has a bell-shape, representing, for example, the
upconing under a well. It will be shown that the Dupuit approximation overestimates the specific discharge
vector and becomes inaccurate when the horizontal size of the bell shape becomes small.
Consider the following hypothetical situation of an aquifer in which the density varies linearly from salt
(t/i = /./s) at z = d to fresh (1/2 = 0) at z = (2- The thickness of the transition zone is constant and equal
to h so that the density in the transition zone may be written as
h " >J- —bz
The elevation of the top and bottom of the transition zone are
Ci = Ae (x/<7>2 -\h C2 = Ci I h (4.3)
where A and a are chosen as: a = 2h, A. = 0.2h and v8 = 0.025.
An exact solution is derived for the case that the aquifer is infinitely thick. It will be shown that the flow
caused by the density variation only has a limited extent in the vertical direction if the resistance to flow in
the vertical direction is not neglected. The continuity equation may be written, with the aid of Darcy's law,
as
dqx dqz d^x &X '^v
or
/,._ (A tn
"• .-) 11-0;
dz
where x = k(/> (equation (2.53)). The term —kdi//dz, and thus the Laplacian of x, equals zero in the fresh
and salt water zones and kvs/h in the transition zone.
The function x is modeled with analytic elements (Strack, 1989). The transition zone is represented by
an area sink of constant extraction rate kvs/h\ the boundary of the transition zone is approximated with
a polygon. The transition zone is infinitely long. Far away, where the transition zone becomes horizontal,
the area sink is approximated by long line sinks of strength ki>s. When z approaches positive infinity, |^
approaches ^ki/s (half the "water" extracted by the area sink comes from the top, the other half from the
bottom of the aquifer); for z approaching negative infinity, -^ approaches —^ki//,. The desired behavior at
-------
. J i
\ I
t
Hq|=0.001k
Figure 4.3: Exact solution for variable density flow in an infinite aquifer with the boundary of the transition
zone consisting of straight segments
infinity is that qz equals zero. In terms of derivatives of x this becomes, with equation (2.52) for qz, and
using that v = 0 in the fresh water zone and v = vs in the salt water zone
lim ^ = 0 lim ^ = -kvs (4.6)
z^+oo dz z^-oo dz
A term x = — \kvsz (of which the Laplacian is zero) is added to the solution to obtain the correct behavior
at infinity. The flow field for the exact solution may now be obtained with Darcy's law (2.52). Note that
the specific discharge vector is not just the gradient of x, but that an extra term — kv must be added to qz.
The flow field and transition zone are shown in Figure 4.3.
An exact solution for an aquifer of finite thickness H may be obtained as follows. (The solution is exact
in the sense that the vertical resistance to flow is not neglected; the boundary condition along the base of
the aquifer will be met approximately.) The area-sink and the two line-sinks that represent the area-sink
far away are imaged through the line z = H/2; this will create an impermeable upper boundary. The
impermeable base is approximated by a long line-sink with a polynomial strength; this line-sink is also
imaged through z = H/2. The coefficients of the polynomial are computed such that dx/dz = —kvs along
the base, using the procedure proposed by Jankovic and Barnes (1997). For this solution the extra term
X = —\kvsz is not needed.
The flow field obtained with the Dupuit approximation for an aquifer of finite thickness H may be
obtained using the expressions presented in the next chapter with (4.2) for the density distribution. The
38
-------
• • ' / / ^ ^^L^MLJ^J*^^ ^ ^ v ' ' '
' f /
-...ft
' \ \
. , v \
t <
"\ ^"^
-It
^ ^ |
1
- > t
—*. — *• ^*
— k- k- — »—
\ ^ ^
ft.--
/ /
s* ' • •
~|q=0.002k
Figure 4.4: Flow field for Dupuit solution (a = H/2, h =
functions ^i and ^2 are given by (4.3) and the derivatives and Laplacians are
f?Cl_£C2_ 2A^-(x/^
^ ~~ ^ ~~ —//i^^e
ax ax a1
ax2 ax2
Solutions are obtained for two cases. Case 1 is characterized by the following values: a = H , h = H/2; case
2 is characterized by: a = H/2, h = H/4:. The flow field for case 2 is shown in Figure 4.4.
The value of qx at z = —h/2 is plotted versus x for the exact (solid line) and the Dupuit solution (dashed
line) for case 2 (see Figure 4.5). The vertical component of flow is plotted versus z at x = 0 (Figure 4.6). The
solid lines represent the exact solution and the dashed lines the Dupuit solution. The thin lines correspond
to case 1, and the thick lines to case 2. The dotted line is the exact solution for an infinite aquifer. It is
concluded from Figures 4.5 and 4.6 that a solution obtained with the Dupuit approximation overestimates
the specific discharge vector. Such a behavior has also been observed for Dupuit interface flow (see, e.g.,
Bear, 1972).
The value of a is a measure of the horizontal size of the upconing; the horizontal size of the upconing
is approximately 4 2a (see equation (4.3)). The
Dupuit solution becomes inaccurate when the horizontal size of the bell shape is smaller than twice the
thickness of the aquifer, as may be seen from Figure 4.6. Further research is needed to draw general
conclusions about the range of applicability of the Dupuit approximation.
39
-------
0.0012
0.0008 -
0.0004-
0.0000 -
-0.0004 -
-0.0008 -
-0.0012
-1
-0.5
0 0.5
x/H
Figure 4.5: Comparison of qx versus x for exact (solid) and Dupuit (dashed) solutions
0.0000
-0.0025
-2
-1
0
z/h
Figure 4.6: Comparison of qz versus z lor exact and Uupuit solutions; (thick lines), H = 4h (thin lines);
exact (solid lines) approximate (dashed lines) exact solution for semi-infinite aquifer (dotted line)
40
-------
Choice of the shape factor A
A procedure is outlined to simulate the change of the salinity distribution over time due to advection of
the salt with the groundwater. If one is interested in relatively small times (on the order of 25 years, as
is the case for most engineering purposes) it might be reasonable to neglect other processes that affect the
salinity distribution, such as diffusion and microscopic dispersion. The flow in the aquifer is approximated
as incompressible.
The elevation of surface n is represented by a multiquadric interpolator function (Hardy, 1971), which
may be written as
M ,
Cn(x1:x2) = ]T an\J(xi ~~ ?i)2 + (z2 - z2)2 + A2 + an (4.8)
m=l
The constant A controls the smoothness of the interpolator. The M I 1 constants an (m = 0,..., M) are
rn
determined from M + 1 conditions. M conditions are that („ equals a specified value Qn at M collocation
points (a: ], #2), and one condition is that the sum of the an (rn = f,..., M) equals zero.
To obtain accurate expressions for the specific discharge vector, the multiquadric interpolator must not
only represent the elevations of the surfaces of constant density accurately, but also the first and second
derivatives of the elevations of the surfaces. This puts a constraint on the choice of the smoothness parameter
A. Consider, for example, flow in the vertical plane where the multiquadric function is a function of one
variable (x). If A is chosen equal to zero, the multiquadric function reduces to a piecewise linear interpolator.
The derivative of a piecewise linear interpolator is discontinuous and the second derivative becomes infinite
at the collocation points and is zero between them. If A is reduced, the derivatives approach the behavior
of the derivatives of a piecewise linear interpolator.
It was found that an accurate representation may be obtained if A is chosen equal to the distance
between the (regularly spaced) collocation points. (In practice, A may be chosen equal to the average
distance between collocation points, for example.) Figure 4.7 shows the function Q\ (solid line) of equation
(4.3) and two multiquadric representations. The collocation points are spaced a distance a/2 apart and A
is chosen equal to a/2 (dashed line) and a/20 (dotted line), respectively. The first and second derivatives
of the function arid the two multiquadric representations are also shown in Figure 4.7. The case for which
A is equal to the distance between the collocation points is almost indistinguishable from the function it
represents; the interpolator with smaller A gives poor results for the derivatives. A similar behavior is
observed if the multiquadric interpolator is a function of two or three variables. It is rioted that a collocation
point is located at the maximum of the function; in practice the location of the maximum is not known and
the representation of the density distribution is less accurate.
41
-------
Transient simulations
Expressions for the movement of the surfaces of constant density through time are obtained by applying
continuity of flow in every region separately. The salt water region is bounded below by the aquifer base and
on top by surface 1. The salt water zone is called region 0 and the discharge vector in the salt water zone is
o
represented by Qi. so that continuity of flow in the salt region gives
(4.9)
dt
where 9 is the effective porosity. A superscript is added to the specific discharge vector; a superscript n
indicates evaluation at z fn. The horizontal components of the specific discharge vector do not vary with
z in the salt water region, so that
4 = 9i(Ci--6) (4.10)
n
Substitution of (4.10) for Qi in (4.9) and differentiation gives
(4.11)
The sum of the latter two terms equals qz(z = Ci) = qz, as may be seen from (2.19), and (4.11) may be
written as
A similar equation may be derived for the movement of surface 2. Continuity of flow in region 1 states
and the divergence of the discharge vector may be written as
(4.14)
i
where use is made of Leibniz's rule. Substitution of (4.14) for diQi in (4.13) and using that
(4.15)
Ci
gives, after rearrangement of terms
/~)A-> '>
-------
Function
x-axis
x-axis
x-axis
Exact
Delta=3igma/20
— — — — Delta^Sigma/2
Figure 4.7: Behavior of the multiquadric interpolator
43
-------
In general, the differential equation for the movement of surface n may be written as
(4.17;
where ~~ is the material time derivative for the two horizontal directions. Equations (4.16) arid (4.17) are
equivalent to the equations for a moving interface as obtained by. e.g., Bear (1972) and De Josselin de Jong
(1981).
The movement of a surface of constant density through time may be approximated by numerical integra-
tion of either (4.17) or (4.18), where the specific discharge vector is taken constant during a time step. At
the end of the time step, a new specific discharge field is computed, based on the new density distribution.
The accuracy may be improved by using a predictor corrector procedure and/or smaller time steps.
For general cases, (4.18) is probably preferred, especially in the presence of leakage or in the neighborhood
of inhomogeneities in the aquifer properties. It ma}' also be beneficial to choose the collocation points for
the multiquadric interpolator on a regular grid. Equation (4.18) should then be integrated by determining
what elevation will arrive at a grid point during a time step, instead of where the elevation of a grid point
moves to.
As an example, the change through time of the density distribution of a hypothetical problem is computed.
Consider an aquifer of hydraulic conductivity k = Irti/d and thickness H = 40m. The aquifer is divided into
3 regions. The density in the salt and fresh regions are vi = 0.025 and 1/9 = 0, respectively. The density in
the transition zone at t = 0 is
—^-i'i Ci
-------
20
Time=0days
20 •
o
-20
-80
Time=500 days
TimeMOOOdays
40
40
|q=0.002k
80
tfj
a
on
.
Figure 4.8: Results of transient simulation
-------
Conclusions
The implications of the Dupuit approximation for variable density flow in coastal aquifers were investigated.
The density distribution was represented by a number of surfaces of constant density; the elevations of the
surfaces were approximated with rnultiquadric interpolators and the density varies linearly between them.
It was shown that the smoothness parameter A in the multiquadric interpolator must be of the order of
the average distance between control points to obtain reasonable results for the gradient arid Laplacian
of the density distribution; this is required to obtain accurate specific discharges (and thus velocities). A
new exact solution was derived for two dimensional, variable density flow in the vertical plane. The exact
solution was compared to the Dupuit solution. The problem chosen for comparison consisted of a bell shaped
transition zone. The comparison showed that the Dupuit approximation overestimates the specific discharge
vector and that the flow field becomes inaccurate when horizontal size of the upconing is smaller than two
times the aquifer thickness for the specific case investigated. Additional research is needed to draw general
conclusions on the range of application of the Dupuit approximation for variable density flow. The new
density distribution was high!}' useful for the comparison of Dupuit solutions to exact solutions. Equations
were derived for the movement of surfaces of constant density through time. The equations were used for a
simple transient simulation.
Acknowledgment
The computations for the exact solutions in Figures 4.3, 4.5, and 4.6 were obtained with the program Split
written by Igor Jankovic.
46
-------
CHAPTER 5
The Specific Discharge Vector
for a Vertically PlecewIse—LInear Density Distribution
Introduction
Expressions are derived for the specific discharge vector resulting from the new density distribution intro-
duced in the previous chapter. First, expressions are derived for the simple case of a transition zone consisting
of two constant density planes (A' = 2). General expressions for arbitrary A' are derived in the second part.
The expressions for the specific discharge vector are reproduced here from Chapter 2 for completeness
kd.j.0
/""")
— 1 «9,:
ti
7±M £ £i
H * ff '
ti
ff
V2
A simple transition zone
The integrations in the expressions for the specific discharge vector may be carried out when a functional
form is chosen for the dimensionless density distribution. As an example, the density is taken to vary linearly
from t/i to t/2 between (,i(x,y) and {^(avjy). For z < Ci the density is equal to v\ arid for z > £2 the density
is i/o. In functional form this becomes
47
-------
The integrations in (5.1) may be carried out without knowing the functions Ci(x;?/) and C,'2(x^y}^ because
they are independent of z.
!/•,
l z 2
It may be verified from (5.4) that the function | vdz is continuous. Differentiation of (5.4) gives
C2
z < Ci
where it is assumed that 25 is constant. Expression (5.5) may he integrated from zi, to zt
(v\ — 4/2)
+ \(v\ - ^2)(d;Ci + diC,i}(zt - Qi]
= ±(z/i -t/2)[i(C2-Ci)(2diCi I d.^) I (-t-C2)(<3iCi I <%)]
The horizontal components of the discharge vector become, with (5.1), (5.5). and (5.6)
']
z>C,'i
Cl < z < C'2
^2 - mr
0 z < Ci
The vertical component of flow is obtained from (5.2) which may be written as
~H
zb
The derivation of V2 J vdz, which requires differentiation of (5.5), is messy for the region Ci <
vector fi(x,y,z) is introduced for convenience as
(5.4)
(5.5)
(5.6)
(5.8)
£2- The
r
9-i / vdz \(v\- i"2)fi Ci < z < C2
(5.10)
48
-------
where g.i is the first fraction in (5.10) and hi is the second fraction. Differentiation of g-i gives
d-igi M — ^-77^—7-1-7;———^' —
And differentiation of hi
This gives for V2 f i/d/:
Integration of (5.14) gives for d
l,S2 - SI
so that
(5.13)
(5.14)
(5.15)
(5.17)
49
-------
Combining the previous three equations gives
V2 A (^ r \ i \7'2 ^ / ~ A \ , 1 f r f \f\~Jil^ \7'2 r \
Ci(a - u; IV (,2\z - (,2) I 3(1,2 - CIAV 1,2 - v (.1)
3(C2-Ci)3
0
< Ci
(5.18)
Differentiation of (5.6) gives
V2 f vdzdz= l(i/| -1/2) [i(^C2-t)lCi)(2«9lCi I ^C2) I 3(C2-Ci)(2V2Ci I V2C2)-
(5.19)
= ~-(/./i - 1/2) [-|(9jC2^iC2 + d-itidiCi + <9iCidjCi)
|(C2 - Ci)(2V2Ci + V2C2) + (zt - C
and consecutive integration
The vertical component of the discharge vector may now be computed with (5.8), (5.18) (with (5.15). (5.16).
and (5.17)) and (5.20) which gives
S2
-------
A general transition zone
Consider a transition zone that varies linearly from v = v\ at z = Ci to v = 1/2 at z = (2, then varies linearly
from v 1/2 at z £2 to v vy, at z £3 and so on until the fresh water at z (j,v- Hence the density
varies linearly in the vertical direction over N — 1 sections, or written as an equation:
C»+ 1 — Cn
The integrations in (5.1) are carried out in what follows. Since the expressions become lengthy, they are not
combined into one expression.
(5.24)
C,v)
The expressions for qz include the integral V2 I vdz. The derivation of V2 J i/dz, which requires differen-
n
tiation of (5.24), is rncssy for the region (*„, < z < Cn+i- The vector fi(x,y,z) is introduced for convenience
-------
as
where /, is
Cn+ i ~ Ci
where gi is the first fraction in (5.27) and hi is the second fraction. Differentiation of g.i gives
™ (,r.
And differentiation of hi
Combination of (5.28) and (5.4) gives
This gives for V2 / vdz
//i¥2Ci - /AvV2(
V2 / vdz = < /yi '
J
0
Integration of (5.30) gives
"^ -i* — ' 1
Cn < z < Cn + 1
2 < U
-------
so that
I if,* , — f V9VV i VV . ,}
i 3vVm—I Sm/V~v ^"i ' v >?ri —I/
Combining equations (5.26) through (5.33) gives
?/-,\7V, f? — iM — VA/ 1 v \?'2r (f ,,—f}
"IV m^-4 my Z^?n=l m ST))ASm + l Sm^
C, (,',v
(5.34)
All integrals are now computed and the specific discharge vector may be computed. The expressions have
been implemented in a FORTRAN program, which has been used in Chapter 4 to assess the implications of
the Dupuit approximation.
-------
CHAPTER 6
Results and Conclusions
The objective of this report is to investigate the performance of the Dupuit theory for variable density
flow combined with analytic elements to model groundwater flow in coastal aquifers. Four areas of study
were identified in the introduction. The conclusions pertaining to these four areas of study are summarized
in this chapter.
Study area 1.
The analytic element modeling of groundwater flow in the first confined aquifer beneath the
Delmarva Peninsula.
An analytic element model was presented for groundwater flow in the shallow aquifers beneath the
Delmarva peninsula, in Chapter 1. The modeling of the Delmarva peninsula is greatly hampered by the
unavailability of measurements of the salinity of the groundwater, especially near the shore. Only six
measurements were found of brackish or salt water. All other measurements indicated fresh water. The six
measurements were insufficient to construct a three-dimensional picture of the salinity distribution below
the Delmarva peninsula. To overcome this problem, publications of salinity in the Chesapeake Bay and the
Atlantic Ocean were used to construct a conceptual model of the salinity distribution below the peninsula.
This distribution results in a horizontal variation from fresh water on the peninsula to salt water in the
bay and ocean. The variation in the vertical direction is assumed to be negligible. This is known to be
unrealistic, but there are not enough data to propose otherwise.
The resulting model of the fresh water head was compared to nine measurements of the head on the eastern
part of the peninsula (the counties of Sussex. DE, Wicomico, MD, and Worchester, MD). The comparison
showed that the simulation of fresh water heads in that area is reasonable, but general conclusions for the
entire model cannot be drawn. Additional data (head, discharge and chloride measurements) are necessary
to calibrate the model and to improve the conceptual model of the density distribution.
The modeling study demonstrated the many challenges in building groundwater models of flow systems
in coastal aquifers. The biggest constraints in the Chesapeake Bay area are the availability of density data
54
-------
and water table measurements. To assess the full capabilities of the approach, it should be applied to an
area where more data are available, perhaps an area outside the United States.
Study area 2
The reduction of computation time by the use of a supercomputer.
The second area of study was addressed in Chapter 3. The Dupuit approximation for variable density flow
was implemented in an analytic element model. The resulting program was written to run on the CrayC916, a
vector machine. The Dupuit formulation for variable density flow, as well as the analytic clement formulation
for groundwater flow, are suited ideally for implementation on supercomputers. Both formulations result in
large sums of complicated functions that are easily vectorizable. The performance of the vectorized code
was an order of magnitude better than the urivectorized code. The vccorizcd code performed at a speed of
over 300 billion floating point operations per second (300 Mflops), which is a significant improvement over
the non-vectorized code. It is noted that high performance computing seems to move away from vector
machines to massively parallel machines. This will have no adverse effect on the implementation of analytic
element codes on supercomputers. The large sums of complicated functions can just as easy be modified to
run on massively parallel machines as on vector machines. An analytic element code written specifically for
use of a massively parallel machine is presented by Haitjema et al. (1997).
Study area 3
The accurate representation of the density distribution.
In Chapters 1. 2 and 3, the density distribution was represented by the three-dimensional multiquadric
interpolator function. It was noted that the shape factor A of the interpolator was chosen close to zero.
The main reason for this choice was to improve control of the interpolator, especially at the transition from
brackish water of variable density to fresh or salt water of constant density. It was shown in Chapter 4 that
this will result in reasonable simulations of the density (and thus the fresh water head, as was the objective
of the modeling study in Chapter 1), but that the resulting velocity vector was unrealistic near the points
where the density is specified. This makes it difficult to simulate the change of the salinity distribution
through time; during such a simulation the salt moves with the groundwater and the velocities have to be
accurate.
A new function to represent the density was introduced to solve this problem. The representation
consists of a number of surfaces of constant density; the elevations of the surfaces are approximated with the
multiquadric interpolator and the density varies linearly between them. It was shown that the smoothness
-------
parameter A in the multiquadric interpolator must be of the order of the average distance between control
points to obtain reasonable results for the velocity field. It is noted that this was not practical for the original
three dimensional interpolator.
Study area 4
The implications of adopting the Dupuit approximation for variable density flow.
The new density distribution was used to compare a Dupuit solution with an exact solution, a solution in
which the vertical resistance to flow is not neglected. The problem chosen for comparison consisted of a bell
shaped transition zone. The comparison showed that the Dupuit approximation overestimates the specific
discharge vector and that the flow field becomes inaccurate when the width of the bell-shaped transition
zone is less than two times the thickness of the aquifer.
The new density distribution represents the transition from variable density zones (brackish water) to
constant density zones (fresh or salt water) accurately. The new distribution does introduce complications
that remain to be solved, however, such as the intersection of the transition zone with the base or top of the
aquifer. The new density distribution is highly useful for the comparison between Dupuit solutions and exact
solutions, because it is relatively easy to obtain exact solutions for flow in the vertical plane corresponding to
the new density distribution. A comparison between exact and Dupuit solutions may be used to determine
the full range of applicability of the Dupuit approximation for variable density flow.
-------
References
Achrtiad, G., and J .M. Wilson, 1993. Hydrogcologic framework and the distribution and movement of brack-
ish water in the Ocean City-Manokin aquifer system at Ocean City, Maryland. Maryland Geological
Survey Report of Investigations No. 57.
Bear, J. 1972. Dynamics of fluids in porous media. Dover Publications, Inc., New York. NY.
De Josselin de Jong. G. 1981. The simultaneous flow of fresh and salt water in aquifers of large horizontal
extension determined by shear flow and vortex theory. In: Flow and transport in porous media. A.
Verruijt and F.B.J. Barends, Editors. A.A. Balkema, Rotterdam, The Netherlands, pp. 75 82.
De Lange, W.J., 1996. Groundwater modeling of large domains with analytic elements, Ph.D. dissertation.
Delft University of Technology, published as RIZA note 96.028, Lelystad.
De Lange, W.J. 1997. On the effects of three-dimensional density variation in groundwater flow on the
head and flow distribution in a multi—aquifer system. In: Conference companion Part 1. International
conference Analytic based modeling of groundwater flow. April 1997. Nunspeet, The Netherlands.
Fleck, W.B., and D.A. Vroblesky, 1996. Simulation of ground-water flow of the Coastal Plain aquifers in
parts of Maryland, Delaware, and the District of Columbia, USGS Professional Paper 1404-J.
Haitjema, H., V. Kelson, and K. Luther. 1997. Modeling three dimensional groundwater flow on a regional
scale using the analytic element method. Report to the US EPA pursuant to project number CR
823637-01-1.
Haitjema, H. 1995. Analytic. Element Modeling of Ground/water Flow. Academic Press, San Diego, CA.
Hardy, R.L. 1971. Multiquadric equations of topography and other irregular surfaces. J. Geophys. Res.
76, 1905 1915.
James R.W., R.H. Simmons and B.F. Strain. 1992. Water Resources Data Maryland and Delaware,
Volume 2. U.S. Geological Survey Water-Data Report MD-DE-92-2, Towson. MD.
-------
Jankovic, I., and R.J. Barnes. 1997. High order line elements for two dimensional groundwater flow. In:
Proceedings of the Analytic-based 'modeling of groundwater flow conference. Nunspeet. The Nether-
lands. April 1997.
Kipp, K.L. Jr. 1986. HST3D. A computer code for simulation of heal and solute transport in three-
dimensional groundwater flow systems. IGWMC, USGS Water Resources Investigations Report 86-
4095.
Konikow, L.F., and J.D. Bredehoeft. 1978. Computer model of two dimensional solute transport and
dispersion in ground water. Techniques and Water Resources Investigations. Book 7.
Konikow. L.F., D.J. Goode, G.Z. Hornberger. 1996. A three dimensional method-of-characteristics solute
transport model (MOC3D), USGS Water Resources Investigations Report 96-67406.
Leahy, P. P., and M. Martin, 1993. Geohydrology and simulation of ground-water flow in the Northern
Atlantic Coastal Plain aquifer system, USGS Professional Paper 1404-K.
Lester, B. 1991. SW1CHA. A three-dimensional finite element code for analyzing seawater intrusion in
coastal aquifers. Version 5.05. Geo'Trans, Inc., Sterling VA.
Lusczynski, N.S. 1961. Head arid flow of ground water of variable density. J. Geophys. Res. 12(66),
4247 4256.
Maas, C and M.J. Rmke. 1988. Solving varying density groundwater problems with a single density
computer program. Natuurwetenschappelijk Tijdschrift 1989, Vol 70, p!43 154.
McDonald, M.G. and A.W. Harbaugh. 1984. A modular three-dimensional finite—difference groundwater
flow model. USGS Open File Report 83-875.
Meisler, !L, 1981. Preliminary delineation of salty ground water in the northern Atlantic Coastal Plain,
USGS Open-File Report 81-71.
Meisler, II., 1989. The occurrence and geochemistry of salty ground water in the Northern Atlantic Coastal
Plain, USGS Professional Paper 1404-D.
Minnema, B. and J.L. Van Der Meij. 1997. Three-dimensional density-dependent groundwater flow. In:
Conference companion Part 1. International conference Analytic based modeling of groundwater flow.
April 1997. Nunspeet, The Netherlands.
Olsthoorn, T.N. 1996. Variable density groundwater flow modelling with MODFLOW. In: 14th Salt Water
Intrusion Meeting, 16-21 June, 1996. Malrno, Sweden. Rapporter och meddelanden nr. 87. Geological
Survey of Sweden, Uppsala, Sweden.
58
-------
Oude Essink, G.H.P., and R.H. Boekelman. 1996. Problems of large scale modelling of salt water intrusion
in 3D. Proc. 14th Salt Water Intrusion Meeting, Mai mo, Sweden, June 1996. pp. 16-31.
Oude Essink, G.H.P. 1998. MOC3D adapted to simulate 3D density-dependent gronndwater flow. Paper
presented at the MODFLOW98 Conference, October 1998, Golden, CO.
Plielan, D.J., 1987. Water levels, chloride concentrations, and pumpage in the coastal aquifers of Delaware
and Maryland, USGS Water-Resources Investigations Report 87-4229.
Richardson, D.L., 1992. Hydrogeology and analysis of the ground-water flow system of the Eastern Shore,
Virginia, USGS Open File Report 91-490.
Strack, O.D.L. 1989. Groundwater Mechanics. Prentice Hall, Rnglewood Cliff's. N.I.
Strack, O.D.L. 1995. A Dupuit-Forchheimer model for three-dimensional flow with variable density. Water
Resources Research, 31 (12), pp. 3007 3017.
Strack, O.D.L., and M. Bakker. 1995. A validation of a Dupuit-Forchheimer formulation for flow with
variable density. Water Resources Research, 31(12), 3019 3024.
Strack, O.D.L., 1997. Principles of the analytic element method. In: Conference companion Part 1.
International conference Analytic-based modeling of groundwater flow. April 1997. Nunspeet, The
Netherlands.
Strack, O.D.L., 1997. In: Conference companion Part 2. International conference Analytic based modeling
of groundwater flow. April 1997. Nunspeet, The Netherlands.
Strack, O.D.L., and I. Jankovic. 2000. A multi-quadric area-sink for analytic element modeling of ground-
water flow. Journal of Hydrology, in press.
Strack, O.E., 1992. MLAEM User's Manual. Strack Engineering, North Oaks. MN.
Sanford. W.E. and L.F. Konikow. 1985. A two-constiuent solute-transport model for ground water having
variable density. USGS Water Resources Investigations Report 85-4279.
Van Darn, J.C., 1973. Lecture notes to Geohydrology (F15b), Delft University of Technology, Department
of Civil Engineering, Delft, The Netherlands.
Van Gerven, M.W.., and C. Maas. 1994. Testing the new variable density module of MLAEM in the dune
water catchment of the Amsterdam Municipal Waterworks (in dutch). KIWA-report SWO 94.212,
KIWA, Nieuwegein, The Netherlands.
-------
Voss, C.I. 1984. SUTRA A finite element simulation for saturated imsaturated , fluid density-dependent
ground-water flow with energy transport or chemically reactive single-species solute transport. USGS
Water Resources Investigations Report 84-4369.
Vroblesky, D.A., and W.B. Fleck, 1991. Hydrogeologic framework of the coastal plain of Maryland.
Delaware, and the District of Columbia, IJSGS Professional Paper 1404-E.
WES, 1997. The Groundwater Modeling System version 2.0, US Army Corps of Engineers, Waterways
Experiment Station, Vicksburg, MS.
Woodruff, K.D., 1969. The occurrence of saline ground water in Delaware aquifers. Delaware Geological
Survey Report of Investigations No. 13.
-------
Appendix Al
x (rn)
424108
425613
425613
425613
425613
410711
425613
431596
428688
428688
430174
381120
442072
433191
439144
439158
437674
439158
442139
439172
436205
436205
378318
440656
353163
474791
474791
482209
470356
479254
344373
335481
488145
488145
480740
480740
489629
482222
371109
445196
467414
480745
482226
482226
491113
492594
492594
480749
480749
449651
y (m)
4204046
4205912
4205912
4205912
4205912
4206061
4205912
4209651
4217197
4217197
4217184
4217752
4218971
4222799
4224633
4226513
4226524
4226513
422837 1
4228393
4228416
4228416
4229077
4228382
4233256
4231964
4231964
4233824
4235739
4237592
4239063
4239239
4237574
4237574
4239468
4239468
4239452
4239465
4242354
4241511
4241391
4241348
4241345
4241345
4241330
4241329
4241329
4243228
4243228
4243362
z (rn)
-347.5
-299.9
-295.4
-151.5
-341.7
-1.1
-106.7
-13.4
-173.7
-342.9
-8.7
-2.1
-9.1
-15.2
-62.5
-16.0
-12.2
-32.3
-7.2
-9.8
-18.3
-19.2
-4.7
-14.9
-19.8
-4.6
-13.7
-7.0
-11.7
-17.5
-1.5
-2.3
-14.3
-27.1
-4.6
-12.3
-14.6
-12.6
-2.9
-8.4
-4.9
-16.5
-16.3
-26.4
-17.5
-22.9
-12 2
-32.6
-32.8
-22.6
P (kg/rnj)
999.3
999.4
999.4
999.4
999.4
999.3
1001.2
999.3
999.2
999.2
999.3
999.3
999.3
999.2
999.3
999.2
999.2
999.3
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.3
999.2
999.2
999.2
999.2
999.2
999.2
999.3
999.2
999.2
999.3
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
x (in)
435778
437272
449394
455317
446405
458315
459804
411298
412806
411303
412806
411462
415888
412864
415964
417576
419020
420555
419166
429597
420484
428115
426700
428115
422266
432745
440205
434285
435764
435808
440286
440191
446245
452316
446282
449394
452393
459804
411298
420484
409815
409815
411318
409815
412747
412806
415888
417390
417390
417390
y (m)
4173864
4173853
4201970
4196295
4200109
4200040
4200032
4121447
4119552
4119567
4119552
4134606
4127040
4125191
4134560
4145824
4140170
4143915
4155209
4149474
4136396
4151366
4160779
4151366
4166460
4168248
4166311
4173876
4171984
4177624
41.77590
4164431
4175669
4 1 92552
4181309
4201970
4205712
4200032
4121447
4136396
4123343
4123343
4123327
4123343
4113912
4119552
4127040
4127025
4127025
4127025
z (m)
-6.7
-4.6
-7.9
-10.4
-7.6
-4.6
-11.3
-36.0
-30.2
-29.6
-31.4
-40.2
-34.1
-41.2
-38.7
-35.1
-34.8
-36.6
-29.9
-45.1
-47.6
-32.0
-41.5
-37.5
-28.0
-31.4
-44.8
-43.3
-41.8
-36.0
-35.1
-36.6
-43.3
-41.2
-36.6
-36.3
-33.8
-37.8
-60.4
-60.4
-57.3
-48.2
-54.3
-63.7
-58.8
-50.3
-64.6
-54.9
-55.8
-55.2
p (kg/m3)
999.2
999.2
999.2
999.3
999.2
999.2
999.2
999.2
999.2
999.2
999.3
999.2
999.2
999.2
999.2
999.3
999.2
999.2
999.2
999.3
999.2
999.2
999.2
999.3
999.3
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.5
999.6
999.5
999.5
999.2
999.2
999.2
999.2
999.2
999.2
61
-------
X
491115
491115
344482
491115
448170
449651
332633
486676
449663
449663
449663
449663
449663
448182
449663
448182
448182
448182
445221
451155
494079
445234
443767
349030
458565
473363
463005
451166
486682
458575
449698
443780
448219
448219
448219
448219
448219
448219
489644
457096
489646
488167
439356
451189
451189
480772
319490
485209
341741
495564
y
4243210
4243210
4244704
4243210
4243372
4243362
4244941
4245097
4245243
4245243
4245243
4245243
4245243
4245252
4245243
4245252
4245252
4245252
4245271
4247114
4246968
4247152
4249042
4250260
4248952
4248890
4248931
4248994
4248857
4250832
4250883
4250922
4250892
4250892
4250892
4250892
4250892
4250892
4250733
4250840
4252613
4252615
4252835
4252754
4252754
4252629
4254631
4252620
4256043
4254487
z
-22.9
-25.9
-4.0
-36.7
-23.2
-14.0
-6.1
-28.0
-7.6
-8.5
-25.6
-13.1
-10.7
-10.7
-9.3
-9.9
-14.0
-8.7
-11.7
-22.1
-82.6
-19.2
-12.2
-1.2
-1.8
-1.8
-28.3
-25.9
-24.4
-1.8
-20.1
-18.3
-45.7
-57.9
-38.7
-36.6
-24.4
-2.0
-12.3
-2.4
-10.7
-27.7
-10.4
-6.1
-21.0
-17.5
-5.8
-21.3
-4.6
-1 2
P
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.3
999.2
999.2
999.2
999.2
999.3
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.3
999.2
999.2
999.2
999.2
999.2
999.2
999.3
X
417390
412864
419166
429597
423588
420484
428115
423709
426634
425137
426734
426784
432837
440286
440218
434210
441845
443324
435969
435969
446245
452316
447810
449394
450837
455338
459804
459804
459804
459804
458315
411298
412806
415888
412864
412864
417576
422054
419166
341188
428115
431204
432715
432730
422266
432745
440205
440205
440205
435808
y
4127025
4125191
4155209
4149474
4147647
4136396
4151366
4160806
4153259
4153273
4164539
4170179
4179528
4177590
4168191
4164476
4186979
4185089
4198304
4198304
4175669
4192552
4186939
4201970
4194441
4200055
4200032
4200032
4200032
4200032
4200040
4121447
4119552
4127040
4125191
4125191
4145824
4143901
4155209
4150692
4151366
4162621
4164488
4166368
4166460
4168248
4166311
4166311
4166311
4177624
z
-54.6
-61.0
-61.9
-65.2
-63.4
-59.1
-54.9
-45.7
-59. 1
-49.4
-56.4
-54.6
-39.6
-56.4
-64.0
-56.1
-46.3
-56. 1
-33.5
-33.5
-64.6
-61.3
-58 8
-50.9
-60.4
-47.2
-72.2
-63.7
-71.0
-77.4
-67.4
-75.6
-77.4
-86.0
-83.8
-96.0
-66.5
-79.9
-81.7
-91.4
-86.3
-68.6
-67.7
-74.7
-51.2
-72.2
-78.9
-74.1
-76.8
-84.4
P
999.2
999.2
1001.2
999.7
999.2
999.2
999.3
999.2
999.3
999.2
999.2
999.2
999.3
999.2
999.2
999.2
999.2
999.2
999.2
999.3
999.2
999.2
999 2
999.3
999.2
999.2
999.5
999.3
999.4
999.4
999.3
999.7
999.3
999.4
999.2
1000.1
999.6
999.2
1007.8
1002.1
1001.9
999.2
999.2
999.2
999.6
999.2
999.2
999.3
999.3
1000.3
-------
X
486691
448243
495564
483737
319575
325514
324052
344772
309269
334463
427591
427591
426130
430595
325718
361163
361163
361163
361163
587112
587112
587112
587112
587112
587112
621912
621912
621912
621912
624253
624253
624253
624253
624253
624253
624253
624253
624253
624253
624253
624253
624253
624253
624253
624253
624253
624253
624253
624253
624253
y
4254497
4254653
4254487
4256383
4258392
4258260
4260174
4259746
426051 1
4261834
4260453
4260453
4262347
4266068
4267663
4266970
4266970
4266970
4266970
4394561
4394561
4394561
4394561
4394561
4394561
4364874
4364874
4364874
4364874
4302710
4302710
4302710
4302710
4302710
4302710
4302710
4302710
4302710
4302710
4302710
4302710
4302710
4302710
4302710
4302710
4302710
4302710
4302710
4302710
4302710
z
-33.2
-6.9
-11.0
-16.8
-5 f^
-2.9
-4.0
-2.0
-5.3
-6.1
-15.7
-4.4
-4.1
-19.8
-3.7
-7.9
-7.0
-9.4
-8.2
-31.4
-73.5
-93.3
-93.0
-102.7
-140.2
-46.9
-64.3
-73.2
-82.9
-64.9
-74.1
-83.2
-92.7
-102.1
-102.1
-130.8
-140.2
-140.2
-149.7
-159.1
-168.9
-187.8
-197.2
-216.1
-235.0
-253.6
-282.2
-310.3
-329.5
-338.9
P
999.2
999.2
999.3
999.2
999.2
999.2
999.3
999.3
999.2
999.2
999.2
1002.1
999.2
999.2
999.2
999.2
999.2
999.6
999.5
1021.6
1021.6
1010.3
1000.4
1001.7
1001 .4
1019.9
1016.9
1014.3
1004.8
1022.7
1024.3
1014.8
1008.0
1006,4
1004.4
1004.5
1004.1
1003.6
1003.9
1004.1
1004.5
1005.9
1006.1
1007.1
1009.4
1007.8
1011.8
1009.2
1014.4
1018.9
X
435808
440286
435969
443259
452316
447846
449394
446405
446405
461901
461901
463727
463727
461901
448528
472855
486014
487792
487792
487979
487792
487792
486665
486665
486665
486665
487297
487320
487320
487320
487320
447699
466391
488562
488562
489844
487987
492643
492100
492643
492100
492396
492990
493300
492863
492863
492863
493794
493794
493794
y
4177624
4177590
4198304
4175689
4192552
4192579
4201970
4200109
4200109
4306639
4306639
4307713
4307713
4306639
4294466
4292711
4293091
4291790
4291790
4290942
4291790
4291790
4289088
4289088
4289088
4289088
4290811
4290250
4290250
4290250
4290250
4287420
4283041
4286213
4286213
4285214
4285866
4285876
4285465
4285876
4285465
4284776
4285238
4284625
4283933
4283933
4283933
4279961
4279961
4279961
z
-65.8
-82.3
-57.9
-86.0
-91.4
-69.8
-69.8
-101.8
-68.3
-52.4
-57.3
-198.7
-233.0
-57.3
0.0
-4.3
-9.5
-8.8
-0.6
-6.9
-8.8
-9.0
4.0
-18.0
-36.0
-14.0
-9.1
-15.2
-15.2
-19.8
-17,4
-19.7
-22.6
-23.0
-24.1
0.8
-9.9
-16.8
-20.7
-16.8
-19.5
-8.5
-16.8
-3.8
-33.2
-10.7
-4.1
-5.7
-13.3
-17.9
P
999.2
999.7
999.8
999.2
999.2
999.2
1001.2
1002.1
1000.3
999.2
999.2
1000.0
1000.0
999.2
999.2
999.2
1008.0
999.3
1000.7
999.2
1000.1
999.5
999.3
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.3
999.2
999.3
999.3
999.2
999.2
1002.4
1007.6
1002.6
1000.5
63
-------
X
665369
665369
665369
665369
665369
665369
665369
665369
665369
665369
665369
665369
690651
690651
690651
690651
690651
690651
690651
690651
690651
690651
690651
508463
508463
508463
508463
508463
508463
508463
508463
532263
532263
532263
532263
532263
532263
532263
532263
532263
532263
532263
532263
532263
532263
532263
449343
462482
462482
462482
y
4325474
4325474
4325474
4325474
4325474
4325474
4325474
4325474
4325474
4325474
4325474
4325474
4315899
4315899
4315899
4315899
4315899
4315899
4315899
4315899
4315899
4315899
4315899
4266409
4266409
4266409
4266409
4266409
4266409
4266409
4266409
4127653
4127653
4127653
4127653
4127653
4127653
4127653
4127653
4127653
4127653
4127653
4127653
4127653
4127653
4127653
4318613
4307222
4307222
4307222
z
-84.1
-102.1
-121.0
-159.1
-178.3
-235.3
-253.6
-291.7
-310.6
-329.5
-348.7
-377.0
-305.1
-308.8
-313.3
-326.8
-364.5
-383.4
-431.0
-457.8
-498.0
-526.7
-555.0
-29.9
-39.0
-48.2
-57.3
-78.9
-92.7
-102.1
-111.6
-93.9
-111.0
-129.8
-139.3
-148.7
-158.2
-177.4
-196.6
-215.5
-233.5
-252.7
-271.6
-300.5
-310.3
-319.7
-137.2
-146.3
-195.1
-237.7
P
1023.5
1021.0
1018.6
1020.6
1033.0
1018.4
1017.8
1017.5
1017.9
1017.8
1018.7
1024.0
1023.5
1024.0
1024.6
1024.2
1037.5
1037.5
1024.3
1021.3
1024.0
1024.0
1024.3
1019.7
1031.5
1011.2
1005.5
999.5
1001.3
999.5
1000.4
1026.8
1021 .0
1017.6
1017.5
1018.4
1019.7
1020.1
1020.5
1021.9
1022.7
1023.8
1024.0
1023.4
1024.6
1036.2
999.2
999.3
1000.0
1000.0
X
493794
493794
446152
445164
445935
447103
446572
445440
445440
446266
474096
493342
493342
493342
493342
493342
493342
493342
493342
493833
493833
493833
493833
493833
493833
493833
494351
494351
494351
494351
494351
450877
478216
478216
478907
492257
493613
449497
480455
480455
420080
421216
406017
429953
437055
447781
454386
461489
468094
479672
y
4279961
4279961
4276339
4277572
4277613
4276882
4277674
4275266
4275266
4270594
4271638
4277288
4277288
4277288
4277288
4277288
4277288
4277288
4277288
4277324
4277324
4277324
4277324
4277324
4277324
4277324
4273698
4273698
4273698
4273698
4273698
4266554
4265970
4265970
4262483
4258722
4269682
4257141
4255757
4255757
4109399
4102296
4098318
4123746
4138449
4158265
4171476
4185682
4196905
4208128
z
-2.7
-2 7
-16.9
-19.5
-18.3
-19.8
-12.2
-19.8
-16.5
-23.6
-18.4
-8.8
-2.7
-4.2
-5.0
-8.8
-10.3
-11.8
-14.1
-4.2
-2.7
-3.4
-4.2
-5.0
-8.0
-9.5
-3.4
-5.7
-6.5
-8.8
-11.8
-1.7
5.5
5.5
-5.8
-17.8
-19.8
-21.3
-17.7
-65.7
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
P
1001.3
1008.2
999.2
999.3
999.2
999.2
999.2
999.2
999.2
999.2
999.2
1009.0
1000.1
1000.5
1003,4
1009.0
1008.9
1006.3
1000.2
999.5
1006.2
1002.7
1001.9
1002.9
1007.1
1010.3
1022.9
1017.6
1021.9
1018.6
1019.4
999.2
999.3
999.3
999.2
999.2
999.3
999.2
999.2
999.3
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
64
-------
X
445513
445513
445513
492063
493692
484073
495022
495022
424959
424959
424959
485824
482322
494761
494761
516177
452647
457483
447299
450352
440191
427199
433276
426371
423192
413230
413131
381495
386349
383287
474432
488153
474407
485926
490649
492203
492203
492203
492963
493622
493622
493622
493622
493622
494771
495139
495139
495139
495139
495546
y
4298316
4298316
4298316
4285982
4266053
4264973
4258646
4258646
4209049
4209049
4209049
4229800
4221801
4253464
4253464
4317494
4203684
4199978
4192577
4192929
4175551
4165917
4166160
4155489
4148222
4137301
4116059
4101481
4132317
4146958
4233178
4232525
4254373
4249917
4248877
4241833
4241833
4241833
4243841
4246664
4246664
4246664
4246664
4246664
4251403
4254651
4254651
4254651
4254651
4254637
z
-103.6
-170.7
-213.4
-75.6
-106.7
-91.4
-88.4
-137.2
-91.4
-152.4
-91.4
-42.7
-48.8
-128.0
-140.2
-106.7
-70.1
-61.0
-67.1
-54.9
-73.2
-54.9
-73.2
-73.2
-61.0
-39.6
-57.9
-42.7
-45.7
-30.5
-89.9
-91.1
-71.6
-90.2
-125.0
-117.0
-124.7
-134.9
-121.9
-119.8
-125.9
-132.0
-141.4
-150.1
-103.9
-105.8
-119.8
-120.4
-114.3
-125.6
P
999.2
999.4
999.6
999.3
999.3
999.3
999.2
999.9
999.2
999.7
999.3
999.3
999.3
999.5
999.6
999.3
999.3
999.3
999.2
999.2
999.3
999.2
999.2
1000.2
999.3
999.2
999.2
1000.5
999.3
999.3
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.3
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.3
999.2
999.2
X
492883
498494
505597
507585
507585
506094
506094
441 1 10
449436
460465
468280
474561
481353
493258
503629
512978
520793
523860
523860
523860
522327
519770
432784
420080
421216
406017
429953
437055
447781
454386
461489
468094
479672
492883
498494
505597
507585
507585
506094
506094
441110
449436
460465
468280
474561
481353
493258
503629
512978
520793
y
4219704
4229291
4243496
4256066
4280355
4291578
4302946
4117217
4132847
4149499
4163596
4177692
4189743
4202304
4212819
4225892
4241010
4255617
4270736
4281620
4292504
4305945
4108747
4109399
4102296
4098318
4123746
4138449
4158265
4171476
4185682
4196905
4208128
4219704
4229291
4243496
4256066
4280355
4291578
4302946
4117217
4132847
4149499
4163596
4177692
4189743
4202304
4212819
4225892
4241010
z
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
P
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
65
-------
X
495546
493571
495582
494158
494158
494158
494158
494158
494360
494962
494962
484226
474380
474380
483958
491043
489843
492246
492246
492246
489952
492978
492978
492978
493678
493678
493761
493761
493761
493761
493761
490663
495633
494202
494340
494989
494725
484163
479438
474489
493639
485941
474385
493540
495589
495589
494151
484139
479328
494596
y
4254637
4256188
4257574
4263600
4263600
4263600
4263600
4263600
4264920
4265185
4265185
4262435
4271053
4271053
4280699
4284865
4241777
4241853
4241853
4241853
4243705
4243920
4243920
4243920
4246382
4246382
4246617
4246617
4246617
4246617
4246617
4248925
4257666
4263601
4265003
4265170
4265417
4262483
4262824
4254337
4246598
4249891
4254362
4256130
4257595
4257595
4263613
4262404
4262832
4268212
z
-141.1
-106.7
-136.4
-105.9
-104.4
-114.0
-107.1
-105.6
-110.5
-107.9
-110.3
-90.8
-63.4
-61.3
-81.4
-65.2
-65.5
-78.8
-81.2
-68.1
-69.5
-78.3
-70.3
-69.7
0.0
0.0
-81.1
-80.8
-80.9
-86.6
-68.4
-78.6
-87.6
-86.3
-67.1
-58.4
0.0
-61.7
-50.8
-50.3
-57.2
-51.5
-22.3
-51.8
-54.1
-61.0
-60.7
-36.0
-25.9
0.0
P
999.2
999.3
999.8
999.3
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.3
999.2
999.2
999.2
999.2
999.2
999.3
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.6
999.2
999.2
999.3
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.3
999.2
999.2
999.2
999.3
999.2
999.2
999.2
999.2
999.2
X
523860
523860
523860
522327
519770
432784
420080
421216
406017
429953
437055
447781
454386
461489
468094
479672
492883
498494
505597
507585
507585
506094
506094
441110
449436
460465
468280
474561
481353
493258
503629
512978
520793
523860
523860
523860
522327
519770
432784
405022
397422
397422
398914
398914
398914
398914
414611
412622
406656
403531
y
4255617
4270736
4281620
4292504
4305945
4108747
4109399
4102296
4098318
4123746
4138449
4158265
4171476
4185682
4196905
4208128
4219704
4229291
4243496
4256066
4280355
4291578
4302946
4117217
4132847
4149499
4163596
4177692
4189743
4202304
4212819
4225892
4241010
4255617
4270736
4281620
4292504
4305945
4108747
4105633
4203579
4223395
4124099
4144413
4163094
4182769
4096541
4171333
4154144
4134825
z
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
P
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1020.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1025.0
1023.0
1018.0
1017.0
1022.0
1021.0
1020.0
1019.0
1025.0
1020.0
1021.0
1021.0
-------
X
493526
493835
483888
486730
486730
486730
486730
486730
486730
487544
488017
493785
485998
495676
484264
491081
488728
487637
488016
488016
474398
411298
414233
414214
414214
414214
414233
415888
412864
417576
419166
428115
426651
422266
432745
440259
435808
440286
412622
y
4270344
4273600
4280701
4289419
4289419
4289419
4289419
4289419
4289419
4290685
4291125
4246695
4249923
4257637
4262431
4284799
4285570
4290618
4291137
4291137
4271104
4121447
4112017
4110137
4110137
4110137
4112017
4127040
4125191
4145824
4155209
4151366
4155139
4166460
4168248
4173830
4177624
4177590
4104427
z
-61.0
-61.0
-50.6
-19.7
-28.7
-34.8
-41.6
-33.2
-28.7
-40.8
-39.9
-23.9
-17.2
-32.8
-23.8
-23.5
-22.0
-18.0
0.0
0.0
0.0
-9.1
-19.5
-15.2
-15.9
-16.8
-11.0
-1.8
-16.8
-6.1
-4.6
-5.2
-7.0
-10.7
-1.5
-0.9
-4.6
1.5
-100.0
p
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.3
999.2
1012,4
999.2
999.2
999.2
999.2
999.2
999.2
999.2
999.3
999.3
999.2
999.2
999.3
999.3
999.3
999.2
999.2
999.2
999.2
999.3
999.7
999.2
999.2
999.2
999.2
1024.0
X
406017
422708
418233
417736
409639
412622
405022
397422
397422
398914
398914
398914
398914
414611
412622
406656
403531
406017
422708
418233
417736
409639
412622
405022
397422
397422
398914
398914
398914
398914
414611
412622
406656
403531
406017
422708
418233
417736
409639
y
4116641
4187669
4196265
4212955
4223681
4104427
4105633
4203579
4223395
4124099
4144413
4163094
4182769
4096541
4171333
4154144
4134825
4116641
4187669
4196265
4212955
4223681
4104427
4105633
4203579
4223395
4124099
4144413
4163094
4182769
4096541
4171333
4154144
4134825
4116641
4187669
4196265
4212955
4223681
z
-25.0
-25.0
-25.0
-25.0
-25.0
-25.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-50.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
-100.0
p
1022.0
1019.0
1018.0
1017.0
1017.0
1024.0
1023.0
1018.0
1017.0
1022.0
1021.0
1020.0
1019.0
1025.0
1020.0
1021.0
1021.0
1022.0
1019.0
1018.0
1017.0
1017.0
1024.0
1023.0
1018.0
1017.0
1022.0
1021.0
1020.0
1019.0
1025.0
1020.0
1021.0
1021.0
1022.0
1019.0
1018.0
1017.0
1017.0
-------
Appendix A2
This appendix contains information on the area element mesh used in the MVAEM model in Chapter 1.
The area rncsh is shown in figure A.I, and the values used in the rnesli arc presented in the following pages.
Resistances are given in davs and heads in meters.
68
-------
Area
Figure A.I. MVAEM area element mesh
69
-------
No.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
utm-X]
433220
440764
435173
428131
432568
441553
441553
450018
449382
446537
432568
422342
425210
415363
427676
406859
414223
415552
416411
423494
385735
398197
409746
412222
415025
408882
413021
419907
423383
426031
390463
394374
421554
493426
488474
482081
490325
482720
420876
432516
420876
480984
464099
465791
465952
448841
449833
463972
447465
442584
2/1
4286003
4297782
4286006
42751 1 1
4266489
4269590
4269590
4280091
4290209
4268735
4266489
4247730
4246293
4234389
4235401
4234854
4244657
4252008
4222745
4207347
4293428
4282902
4273196
4277673
4282082
4287449
4291294
4297912
4297662
4305443
4298992
4308027
4263514
4282994
4275356
4272564
4256887
4251518
4201229
4197233
4201229
4287644
4323201
4320773
4324945
4345494
4342514
4266612
4312086
4330811
X-2
435173
447465
429761
432568
435973
451967
446537
451967
450598
447365
435973
425210
422342
414223
427671
420660
409693
411914
432937
400658
385735
398956
410677
416636
416579
410936
413997
421593
424966
428132
402007
392592
424029
487736
490122
482843
494895
487595
439812
442706
439812
488291
467777
464099
464191
465310
449350
470712
457771
448142
2/2
4286006
4312086
4274915
4266489
4265048
4279639
4268735
4279639
4290322
4270023
4265048
4246293
4247730
4244657
4229854
4233973
4244589
4253217
4224553
4212037
4276099
4272169
4268704
4273992
4282082
4285544
4288474
4299495
4296877
4305133
4293151
4327831
4265990
4277542
4270814
4271238
4256450
4249699
4205689
4199795
4205689
4291960
4326581
4323201
4323184
4345494
4340753
4273000
4314661
4338509
X3
440764
446067
428131
435344
441553
450018
447365
450598
450582
452230
428815
428815
415363
409693
418346
416411
411914
424043
431594
406365
398956
410677
416636
416579
413997
413997
423383
419360
428132
431988
404276
400711
423597
490122
483730
477098
494158
488372
442706
444106
431594
472855
469842
456030
448142
465310
465952
480634
460740
464191
2/3
4297782
4312409
4275111
4268341
4269590
4280091
4270023
4290322
4298617
4268919
4256138
4256138
4234389
4244589
4227663
4222745
4253217
4266090
4213072
4232987
4272169
4268704
4273992
4282082
4288474
4288474
4297662
4313699
4305133
4315315
4308027
4330306
4270681
4270814
4269772
4270741
4246224
4251707
4199795
4197909
4213072
4307687
4323878
4316628
4338509
4342844
4324945
4260083
4311384
4323184
X4
439241
439241
433220
429761
440610
440610
443040
449382
448931
451862
427133
427133
420660
406859
420660
406365
415552
425354
423494
416411
398197
409796
412222
415025
410936
413021
421593
416874
426031
429571
394374
404226
421061
495219
482081
476238
486188
483847
432516
438696
423494
470175
465791
457771
449350
449833
465310
473403
455090
456030
2/4
4297942
4297942
4286003
4274915
4271652
4271652
4271006
4290209
4296935
4267999
4256644
4256644
4233973
4234854
4233973
4232987
4252008
4264322
4207347
4222745
4282902
4273103
4277673
4282082
4285544
4291294
4299495
4313389
4305443
4315556
4308027
4308073
4271751
4271111
4272564
4271706
4245921
4253880
4197233
4192215
4207347
4300987
4320773
4314661
4340753
4342514
4342844
4255105
4303334
4316628
resist.
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
1.00E4
1.00E3
1.00E2
1 .OOE4
1.00E3
1.00E2
1.00E4
10
40
10
10
10
10
10
10
10
10
head
8
12
3
2
1
2
9
o
12
3
1
0
0
0
0
0
1
2
0
0
0
0
0
0
0
0
0
1.5
1
2
0
0
3
0
0
0
0
0
0
1
1
2
2.5
3
2
3
3
12
15
10
-------
PolylD
X4
2/4
resist. head
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
431988
423383
413997
416636
423597
428647
415363
441255
448931
441255
430805
450598
456743
435414
446919
447181
434239
420660
428815
434239
440383
446537
452230
463972
483498
477427
480634
483115
482843
483730
483730
492347
495541
488212
488212
494158
486188
478045
478045
431988
446067
429571
439236
400659
419360
411914
416636
438696
428517
411090
4315315
4297662
4288474
4273992
4270681
4259425
4234389
4288711
4296935
4288711
4273642
4290322
4290675
4240535
4245303
4254113
4244910
4233973
4256138
4244910
4258997
4268735
4268919
4266612
4284490
4279538
4260083
4270806
4271238
4269772
4269772
4263908
4265045
4262439
4262439
4246224
4245921
4243868
4243868
4315315
4312409
4315556
4338338
4330364
4313699
4253217
4273992
4192215
4174922
4140378
446067
436075
416579
432568
424055
425364
414223
437244
450018
440610
429761
450582
450598
434239
445262
447549
427676
427676
435973
440220
446537
440336
447549
456043
488474
482081
472174
480569
483115
488212
490122
490122
494895
483847
482720
485965
487595
470810
476553
442584
447465
437303
448651
410405
429571
410677
421061
441421
431054
413876
4312409
4291403
4282082
4266489
4266131
4264259
4244657
4290234
4280091
4271652
4274915
4298617
4290322
4244910
4247235
4255888
4235401
4235401
4265048
4258949
4268735
4258949
4255888
4254600
4275356
4272564
4271117
4266813
4270806
4262439
4270814
4270814
4256450
4253880
4251518
4231759
4249699
4245395
4240827
433081 1
4312086
4333126
4345402
4338000
4315556
4268704
4271751
4189956
4172799
4139660
436075
430143
428131
428131
428567
415535
415504
447465
440610
435409
437244
455090
451967
445262
447181
457240
427662
434239
440383
447549
441553
447549
457240
470389
482081
472174
476238
481431
480569
485238
492347
495219
490325
488372
473403
480777
480269
473403
480777
456030
457771
442584
449833
411811
423383
416636
421542
431054
422250
422250
4291403
4279426
4275111
42751 1 1
4259506
4252033
4251845
4312086
4271652
4268344
4290234
4303334
4279639
4247235
4254113
4256467
4229685
4244910
4259112
4255888
4269590
4255888
4256467
4244091
4272564
4271117
4271706
4265829
4266813
4261525
4263908
4271111
4256887
4251707
4255105
4234875
4252477
4255105
4234875
4316628
4314661
4330811
4342514
4335824
4297662
4273992
4263299
4172799
4157769
4157769
424966
413997
430143
416579
432568
422342
422342
455090
441255
430805
441255
456743
457437
446919
449258
456043
435414
428815
434239
445262
435973
451862
463972
473403
477427
470712
485238
483730
477098
481431
488212
495541
488212
490325
480634
486188
478045
480269
486188
446067
456030
431988
448142
401764
421593
416392
416300
428517
419220
419220
4296877
4288474
4279426
4282082
4266489
4247730
4247730
4303334
4288711
4273642
4288711
4290675
4285020
4245303
4254205
4254600
4240535
4256138
4244910
4247235
4265048
4267999
4266612
4255105
4279538
4273000
4261525
4269772
4270741
4265829
4262439
4265045
4262439
4256887
4260083
4245921
4243868
4252477
4245921
4312409
4316628
4315315
4338509
4325542
4299495
4257966
4258057
4174922
4158598
4158598
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
20
10
10
10
10
10
10
20
20
100
80
80
90
90
100
1.00E3
1 .OOE4
1 .OOE4
5.00E3
1.00E2
1.00E4
5.00E3
1.00E3
1.00E4
10
10
10
10
10
10
10
10
1.00E4
2.00E4
2.00E4
15
10
8
6
3
2
1
15
12
12
9
13
10
3
5
12
1
1.5
3
9
10
11
11
15
3
3
10
3
5
4
2
1
1
3
10
2
f
O
10
6
10
5
10
3
0
3
0
1
1
1
1
71
-------
PolylD
X4
2/4
resist. head
101
102
103
104
105
106
107
108
109
110
1 1 1
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
411090
409620
413278
424201
432877
440853
452474
445042
437869
427825
415757
413278
424257
432877
440853
449889
449919
448483
445390
446364
443081
452380
432937
429864
416411
439064
444910
443822
448483
444910
446919
446919
449919
470810
451269
457631
450679
457631
461588
467122
462487
462487
470389
472596
462907
472596
470130
464941
459924
438696
4140378
4121694
4121933
4142110
4158598
4171804
4186880
4170312
4157771
4141935
4115128
4121933
4142110
4158598
4171804
4218087
4227346
4231160
4216573
4214735
4211482
4218954
4224553
4224174
4222745
4225314
4225673
4230736
4231160
4225673
4245303
4245303
4231399
4245395
4227621
4230131
4237662
4230131
4225893
4232488
4237137
4237137
4244091
4229050
4225179
4229050
4217071
4206973
4197567
4192215
413876
413278
415757
427825
437702
445042
456167
448522
440806
431124
421124
413876
415979
422250
431054
452380
451269
449919
452380
445390
439596
461588
439064
427775
418346
444910
449919
441225
443822
441225
452823
452823
449373
478045
449919
461588
452823
462487
462907
468269
464573
464573
476510
476510
468269
476553
465545
461415
456312
444106
4139660
4121933
4115128
4141935
4157548
4170312
4184679
4168489
4155449
4141336
4113505
4139660
4144037
4157769
4172799
4218954
4227621
4231399
4218954
4216573
4212270
4225893
4225314
4229826
4227663
4225673
4227346
4243079
4230736
4243079
4242397
4242397
4237968
4243868
4231399
4225893
4242397
4237137
4225179
4231882
4248351
4248351
4240852
4240852
4231882
4225958
4219017
4208465
4198367
4197909
413278
415757
427825
437702
445042
452474
448522
440806
431124
421124
416714
415979
422250
431054
441421
451269
449919
449373
453457
439596
439812
462907
439812
435414
427773
445390
449889
446919
444910
435414
456043
450679
450679
476553
455165
452380
462487
467122
468269
470389
468864
456043
472596
480777
472596
470130
461415
456312
449279
445600
4121933
4115128
4141935
4157548
4170312
4186880
4168489
4155449
4141336
4113505
4106787
4144037
4157769
4172799
4189956
4227621
4231399
4237968
4217245
4212270
4205689
4225179
4205689
4240535
4229893
4216573
4218087
4245303
4225673
4240535
4254600
4237662
4237662
4240827
4232682
4218954
4237137
4232488
4231882
4244091
4245130
4254600
4229050
4234875
4229050
4217071
4208465
4198367
4188040
4195834
409620
413158
424201
432877
440853
449279
445042
437869
427825
415757
413278
424257
432877
440853
449279
449919
448483
447818
446364
443081
442706
453457
431594
439064
429864
439596
445390
448483
449919
439064
449258
447818
455165
470389
457631
451269
457631
461588
467122
468864
467122
452823
468269
476553
469785
465545
464941
459924
452474
441421
4121694
4108616
4142110
4158598
4171804
4188040
4170312
4157771
4141935
4115128
4108616
4142110
4158598
4171804
4188040
4227346
4231160
4238221
4214735
4211482
4199795
4217245
4213072
4225314
4224174
4212270
4216573
4231160
4227346
4225314
4254205
4238221
4232682
4244091
4230131
4227621
4230131
4225893
4232488
4245130
4232488
4242397
4231882
4225958
4225093
4219017
4206973
4197567
4186880
4189956
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
1.00E4
150
100
80
1.00E3
1.00E3
1.00E3
1.00E3
80
50
50
80
60
30
60
40
80
80
100
1.00E3
100
200
100
200
1.00E3
1.00E3
500
200
1.00E4
1.00E4
1.00E4
1.00E4
1 .OOE4
1.00E4
1.00E4
1.00E4
1
1
1
1
1
1
0
0
0
0
0
7
7
7
10
3
5
8
1.5
1
0
2
1
2
1
/I
10
10
10
4
/I
10
12
5
10
7
10
0
3
4
10
15
5
3
5
3
3
3
3
3
-------
PolylD
X4
2/4
resist. head
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
441421
448492
454623
453457
443081
443081
480777
476553
470130
464941
459924
488143
477747
467990
406859
409693
411914
400418
387155
387228
402007
402007
402007
404294
412625
437303
400658
420876
405851
400658
375045
425522
425522
416317
407999
405363
406509
416714
421124
431124
440806
448522
462629
468441
450889
485965
493280
494158
491815
450889
4189956
4199888
4209957
4217245
4211482
4211482
4234875
4225958
4217071
4206973
4197567
4291984
4302691
4322463
4234854
4244589
4253217
4265167
4268619
4268619
4293151
4293151
4293151
4308050
4335168
4333126
4212037
4201229
4231213
4212037
4241675
4177353
4177353
4159357
4140773
4121334
4107141
4106787
4113505
4141336
4155449
4168489
4196653
4205394
4167283
4231759
4233811
4246224
4294130
4167283
448492
454623
453457
462907
454623
454623
485965
481787
473829
468441
462629
491815
472855
474924
394994
402835
401852
383757
399588
383830
390371
415025
408882
419907
439236
439236
393307
432516
383622
393307
394994
428538
428420
419220
411090
406509
413158
420133
431124
434714
442841
450889
465437
470108
457000
493280
485965
500364
498442
457000
4199888
4209957
4217245
4225179
4209957
4209957
4231759
4224492
4215579
4205394
4196653
4294130
4307687
4312862
4247458
4250013
4256813
4262334
4268659
4262334
4298879
4282082
4287449
4297912
4338338
4338338
4193969
4197233
4238669
4193969
4247458
4174889
4174877
4158598
4140378
4107141
4108616
4104958
4141336
4140019
4153533
4167283
4195349
4204359
4164100
4233811
4231759
4247658
4296478
4164100
456312
461415
465545
469785
453457
445600
481787
473829
468441
462629
456167
481067
474924
477342
402835
401852
398956
389766
398956
400418
385438
409739
419907
412625
429571
448142
413418
425522
377151
371584
406859
438696
419220
411090
409620
413158
420133
423039
434714
442841
450889
459333
470108
481787
476293
500364
470108
502308
480663
440525
4198367
4208465
4219017
4225093
4217245
4195834
4224492
4215579
4205394
4196653
4184679
4306013
4312862
4314349
4250013
4256813
4272169
4246083
4272169
4265167
4293432
4273112
4297912
4335168
4315556
4338509
4185822
4177353
4219258
4203157
4234854
4192215
4158598
4140378
4121694
4108616
4104958
4112155
4140019
4153533
4167283
4183566
4204359
4224492
4199583
4247658
4204359
4268832
4316491
4137775
449279
456312
461415
465545
446364
442706
476553
470130
464941
459924
452474
477747
481067
469842
409693
411914
410677
402835
388507
399588
398197
398197
404294
401663
416874
442584
423494
413418
400658
377001
405851
432516
416317
407999
405363
409620
419107
421124
423039
440806
448522
456167
468441
473829
470108
494158
476293
495693
474924
434714
4188040
4198367
4208465
4219017
4214735
4199795
4225958
4217071
4206973
4197567
4186880
4302691
4306013
4323878
4244589
4253217
4268704
4250013
4275160
4268659
4282902
4282902
4308050
4325484
4313389
4330811
4207347
4185822
4212037
4218731
4231213
4197233
4159357
4140773
4121334
4121694
4102497
4113505
4112155
4155449
4168489
4184679
4205394
4215579
4204359
4246224
4199583
4267759
4312862
4140019
1.00E4
1.00E4
1.00E4
1 .OOE4
1.00E4
1.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
20
20
10
10
10
10
10
10
10
10
10
10
10
10
10
5.00E3
10
1.00E2
1.00E3
10
1 .OOE4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
2.00E4
5.00E1
5.00E4
7
7
5
3
5
3
0
0
0
0
0
0
0
0
1
2
2
1
1
0
1
3
4
4
4
6
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
-------
PolylD
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
utm-xi
434714
429710
440525
457000
479510
494989
500364
502308
519738
512928
512928
491975
465468
474431
491528
506398
529455
567906
612425
589800
573296
467777
465310
469380
469955
485914
489638
489638
393307
413418
416317
404122
383458
397587
395009
395757
406509
419107
429710
439432
474431
460740
467098
471457
476781
455090
456743
466853
471637
460398
2/1
4140019
4098933
4137775
4164100
4206078
4236965
4247658
4268832
4306994
4269879
4247443
4200754
4157473
4084078
4147773
4195296
4242059
4232252
4221629
4172334
4131674
4326581
4342844
4342566
4328638
4310536
4316792
4316792
4193969
4185822
4159357
4164744
4145046
4143113
4119457
4099111
4107141
4102497
4098933
4064416
4084078
4311384
4304398
4299376
4292792
4303334
4290675
4293872
4287288
4276664
a; 2
440525
447060
457704
465468
491976
506082
512928
512928
558322
536933
529455
506398
491528
523720
529495
546581
558322
618599
567906
546581
529495
469955
469380
474613
474613
489638
494830
506680
388705
404122
408114
397587
379270
408010
379270
395009
395757
412878
423784
447060
523720
458976
463741
466853
483498
458976
457437
457437
477907
452230
2/2
4137775
4090786
4129326
4157473
4200649
4232252
4247443
4269879
4313030
4261353
4242059
4195296
4147773
4085044
4141445
4187639
4313030
4301341
4232252
4187639
4141445
4328638
4342566
4342649
4332216
4316792
4333521
4301080
4167413
4164744
4140901
4143113
4119457
4141058
4119457
4119457
4099111
4088811
4072586
4090786
4085044
4308865
4299428
4293872
4284490
4308865
4285020
4285020
4279908
4268919
x?>
429710
457704
465468
491976
506082
512928
512928
519738
536933
529455
506398
491528
474431
529495
546581
567906
618599
659471
546581
529495
523720
480663
469955
474613
489638
506680
474613
528231
404122
416317
397587
383458
395009
405363
379224
405363
412878
423784
439432
474431
523720
463741
466853
471637
477907
463741
466853
460398
463972
443040
2/3
4098933
4129326
4157473
4200649
4232252
4247443
4269879
4306994
4261353
4242059
4195296
4147773
4084078
4141445
4187639
4232252
4301341
4287246
4187639
4141445
4085044
4316491
4328638
4332216
4316792
4301080
4342649
4314682
4164744
4159357
4143113
4145046
4119457
4121334
4099276
4121334
4089071
4072586
4064416
4084078
4043626
4299428
4293872
4287288
4279908
4299428
4293872
4276664
4266612
4271006
X4
419107
440525
457000
479510
494989
500364
502308
498442
5 1 2928
512928
491975
465468
447060
491528
506398
529455
567906
612425
589800
573296
557010
477342
465952
469955
485914
498442
474613
494830
413418
425522
404122
388705
397587
395009
395757
406509
419107
429710
447060
467044
467044
467098
471457
476781
471637
456743
463741
471637
460398
457437
2/4
4102497
4137775
4164100
4206078
4236965
4247658
4268832
4296478
4269879
4247443
4200754
4157473
4090786
4147773
4195296
4242059
4232252
4221629
4172334
4131674
4083364
4314349
4324945
4328638
4310536
4296478
4332216
4333521
4185822
4177353
4164744
4167413
4143113
4119457
4099111
4107141
4102497
4098933
4090786
4054705
4054705
4304398
4299376
4292792
4287288
4290675
4299428
4287288
4276664
4285020
resist.
5.00E4
1.00E5
1.00E5
1 .OOE5
1.00E5
5.00E4
5.00E4
1 .OOK5
1 .OOE5
2.00E5
2.00E5
2.00E5
2.00E5
5.00E5
5.00E5
5.00E5
5.00E5
1 .OOE6
1.00E6
1.00E6
1.00E6
5.00E1
5.00E1
5.00E1
5.00E1
5.00E3
1.00E5
1 .OOE6
1.00E4
2.00E4
2.00E4
1 .OOE4
1 .OOE4
2.00E4
1.00E4
2.00E4
5.00E4
5.00E4
1.00E5
1.00E5
5.00E5
10
10
10
10
10
10
10
10
10
head
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
10
0
0
0
0
0
0
0
0
0
0
0
0
0
8
9
6
9
15
16
17
15
15
74
-------
PolylD
251
252
253
254
255
256
257
258
259
utm-xi
452230
457771
480984
492440
491815
498625
493426
460740
456101
2/1
4268919
4314661
4287644
4290407
4294130
4295111
4282994
4311384
4184726
a; 2
460398
460740
483498
483498
492440
500282
500282
468742
459354
2/2
4276664
4311384
4284490
4284490
4290407
4283393
4283393
4321427
4183642
x?>
463972
468742
492440
487736
498625
493426
502308
474924
465437
2/3
4266612
4321427
4290407
4277542
4295111
4282994
4268832
4312862
4195349
X4
457183
467974
491815
493426
498442
492440
495693
470175
462629
2/4
4268048
4322469
4294130
4282994
4296478
4290407
4267759
4300987
4196653
resist.
10
10
5.00E2
1 .OOE3
1.00E3
1.00E4
2.00E4
10
2.00E4
head
16
4
3
0
0
0
0
2
0
-------
Appendix B
The FORTRAN code of the variable density module is included in this appendix. The variable density
module consists of two files: vardcns.cmn, which contains the common block of data used in the variable
density module, and vardens.for, which contains the functions and subroutines for the computation of variable
density flow. Input of density data makes use of the package Match, developed and written by P. A. Cundall,
which is explained in detail in Strack (1989). The cdir$ commands are vectorization directives for the
CrayC916.
Vardens.cmn
c NVDMX: Maximum number of density points
c NVDPTS: Number of density points
c DVDMUO: Additive constant nuO
c DVDBE: Scale factor beta
c DVDBSQ: Scale factor beta squared
c DVDX: Array with x-values of data points
c DVDY: Array with y-values of data points
c DVDZ: Array with z-values of data points
c DVDAL: Array with alpha values
c DVDDEL: Array with delta values
c DVDDELSQ: Array with delta square values
c DVDNUG: Array with given nu-values
c DVDLU: Matrix stored after LU decomposition
c INDX: Integer array used for LU decompostion
c DUM1,2,3: Dummy arrays for summation
PARAMETER (NVDMX=200)
PARAMETER (DZERCNO. DO, DOME= 1. DO, DTWO2. DO, DTHREE=3. DO )
PARAMETER (D1D6=1.DO/6.DO,DHALF=0.5DO,D1D3=1.DO/3.DO)
COMMON /CBVD/ NVDPTS,DVDNUO,DVDBE,DVDBSQ,
.DVDX(NVDMX),DVDY(NVDMX),DVDZ(NVDMX),
.DVDAL(NVDMX+1),DVDDEL(NVDMX),DVDDELSQ(NVDMX),DVDNUG(NVDMX),
.DVDLU(NVDMX+1,NVDMX+1),INDX(NVDMX+1),
.DUM(NVDMX),DUM1(NVDMX),DUM2(NVDMX),DUM3(NVDMX)
Vardens.for
BLOCK DATA BDVD
IMPLICIT REAL*8 (D)
INCLUDE 'vardens.cmn'
DATA NVDPTS,DVDBE /0,1/
END
SUBROUTINE VDINPUT
IMPLICIT REAL*8 (D)
IMPLICIT CHARACTER*! (A), LOGICAL (L)
SAVE
CHARACTER*32 AFILE
INCLUDE 'vardens.cmn'
INCLUDE 'MATCH.CMN'
DIMENSION AWORD(30)
DATA AWORD /'B','E','T','A',' ','S','0','L','V,' ',
'C' '0','N','T' ' ' ,'C' '0' 'I','N' ' ',
-------
'T'.'E'.'SVTV VQVU'.'IVT',' !'/
LERROR=.FALSE.
LMISS=.FALSE.
10 IF(LMISS.OR.LERROR) WRITE(9,9001)
LERROR=.FALSE.
LMISS=.FALSE.
MRITE(9,9005)
READ(8,9000) ALINE
CALL TIDY
IF(LMISS.OR.LERROR.OR.ILPNT(2).EQ.O) GOTO 10
IF (ILPNT(4).EQ.O) THEN
CALL MATCH(AWORD,!,JUMP)
IF (LERROR) GOTO 10
GOTO (100,200,300,400,500,600),JUMP
100 DVDBE=D¥AR(2)
DVDBSQ=D¥DBE**2
GOTO 10
200 CALL VDSOLVE
GOTO 10
300 CALL VDCONTROL
GOTO 10
400 DCOIN=DVAR(2)
CALL VDCOIN(DCOIN)
GOTO 10
500 GOTO 10
600 RETURN
ENDIF
NVDPTS=NVDPTS+1
DVDX(NVDPTS)=DVAR(1)
DVDY(N¥DPTS)=DVAR(2)
D¥DZ(NVDPTS)=DVAR(3)
DVDNUG(NVDPTS)=DVAR(4)
DVDDEL(NVDPTS)=DVAR(5)
DVDDELSQ(NVDPTS)=DVDDEL(NVDPTS)**2
IF(LMISS.OR.LERROR) NVDPTS=NVDPTS-1
GOTO 10
9000 FORMAT(80A1)
9001 FORMAT(' ERROR: ILLEGAL COMMAND OR MISSING PARAMETERS')
9005 FORMAT
.(' (x,y,z,nu,delta)..(b)...... '
. J(tol)..',/)
RETURN
END
SUBROUTINE VDCOIN(DCOIN)
IMPLICIT REAL*8 (D)
SAVE
INCLUDE 'vardens.cmn'
DTOL=DCOIN*DCOIN
DO 1=1,NVDPTS-1
DO J=I+1,NVDPTS
DIST=(DVDX(I)-DVDX(J))**2+
(DVDY(I)-DVDY(J))**2+
(DVDZ(I)-DVDZ(J))**2
IF (DIST.LT.DTOL) MRITE(7,9000)1,J,DSQRT(DIST)
ENDDO
-------
ENDDO
9000 FORMAT(' POINTS ',14,' AND ',14,' ARE A DISTANCE ',
1E13.5,' FROM EACH OTHER')
RETURN
END
SUBROUTINE VDCHECK
IMPLICIT REAL*8 (D)
SAVE
INCLUDE 'vardens.cmn'
MRITE(7,9000)DVDBE
WRITE(7,9001)
DO I=1,NVDPTS
WRITEC7,9002)1,DVDX(I),DVDY(I),DVDZ(I),DVDNUG(I),DVDDEL(I)
ENDDO
9000 FORMAT(' BETA: ',1P,E13.5)
9001 FORMAT(
.'IX Y Z NU DELTA')
9002 FORMAT(I5,1P,5E13.5)
RETURN
END
REAL*8 FUNCTION DFVDNUM(DX,DY,DZ,M)
IMPLICIT REAL*8 (D)
SAVE
INCLUDE 'vardens.cmn'
DNUM=DVDBSQ*((DX-DVDX(M))**2+(DY-DVDY(M))**2)+
(DZ-DVDZ(M))**2+DVDDELSQ(M)
DFVDNUM=DSQRT(DNUM)
RETURN
END
REAL*8 FUNCTION DFVDNU(DX,DY,DZ)
IMPLICIT REAL*8 (D)
SAVE
INCLUDE 'vardens.cmn'
DO M=1,NVDPTS
DNUM=DSQRT( DVDBSQ*((DX-DVDX(M))**2+(DY-DVDY(M))**2)+
(DZ-DVDZ(M))**2+DVDDELSQ(M) )
DUM(M)=DVDAL(M)*DNUM
ENDDO
DVDNU=O.DO
CDIR$ NOVECTOR
DO M=1,NVDPTS
DVDNU=DVDNU+DUM(M)
ENDDO
CDIR$ VECTOR
DFVDNU=DVDNU+DVDNUO
RETURN
END
SUBROUTINE VDCONTROL
IMPLICIT REAL*8 (D)
SAVE
INCLUDE 'vardens.cmn'
DALTOT=DZERO
78
-------
DO I=1,NVDPTS
DALTOT=DALTOT+DVDAL(I)
DNU=DFVDNU(DVDX(I),DVDY(I),DVDZ(I))
WRITEC7,9000)1,DVDAL(I),DNU,DVDNUG(I)
ENDDO
MRITE(7,9001)DVDNUO,DALTOT
9000 FORMAT(' I,ALPHA,NUCOMPUTED.NUGIVEN ',14,1P.3E13.5)
9001 FORMAT(' NUO,SUM OF ALPHA-S '.1P.2E13.5)
RETURN
END
REAL FUNCTION RFNUGRID(CZ)
IMPLICIT COMPLEX (C), REAL*8 (D)
SAVE
DX=REAL(CZ)
DY=AIMAG(CZ)
DZ=DFELEV()
RFNUGRID=DFVDNU(DX,DY,DZ)
RETURN
END
SUBROUTINE VDSOLVE
IMPLICIT REAL*8 (D)
SAVE
INCLUDE 'vardens.cmn'
NEQ=NVDPTS+1
DO I=1,NVDPTS
DO J=1,NVDPTS
DVDLU(I,J)=DFVDNUM(DVDX(I),DVDY(I),DVDZ(I),J)
ENDDO
DVDLU(I,NEQ)=DONE
ENDDO
DO J=1,NVDPTS
DVDLU(NEQ,J)=DONE
ENDDO
DVDLUCNEQ,NEQ)=DZERO
CALL LUDCMP(DVDLU,NEQ,NVDMX+1,INDX,DUM)
DO I=1,NVDPTS
DVDAL(I)=DVDNUG(I)
ENDDO
DVDAL(NEQ)=DZERO
CALL LUBKSB(DVDLU,NEQ,NVDMX+1,INDX,DVDAL)
DVDNUO=DVDAL(NEQ)
RETURN
END
SUBROUTINE VDSPECDIS2(DX,DY,DZ,DQX,DQY,DQZ)
IMPLICIT REAL*8 (D)
SAVE
INCLUDE 'vardens.cmn'
DK=DFAQK()
DH=DFAQH()
DZB=DFAQZB()
-------
DZT=DFAQZT()
DTOL=l.d-8
DO M=1,NVDPTS
DXMSq=DVDBSq*(DX-DVDX(M))**2
DYMSq=DVDBSq*(DY-DVDY(M))**2
DZM=DVDZ(M)
DZMSQ=(DZ-DZM)**2
DZTMSq=(DZT-DZM)**2
DZBMSQ=(DZB-DZM)**2
DELSQM=DVDDELSQ(M)
DRM=DSQRT(DXMSQ+DYMSQ+DZMSQ+DELSQM)
DRMT=DSqRT(DXMSq+DYMSq+DZTMSq+DELSqM)
DRMB=DSQRT(DXMSQ+DYMSQ+DZBMSQ+DELSQM)
DTESTSQ=DXMSQ+DYMSQ+DELSQM
IF (DTESTSq.LT.DTOL*DZMSq .AND. DZ.LT.DZM) THEN
DP1=0.5DO*DTESTSQ/(DZM-DZ)
ELSE
DP1=DZ-DZM+DRM
ENDIF
IF (DTESTSq.LT.DTOL*DZBMSQ .AND. DZB.LT.DZM) THEN
DP2=0.5DO*DTESTSQ/(DZM-DZB)
ELSE
DP2=DZB-DZM+DRMB
ENDIF
IF (DTESTSq.LT.DTOL*DZTMSq .AND. DZT.LT.DZM) THEN
DP3=0.5DO*DTESTSQ/(DZM-DZT)
ELSE
DP3=DZT-DZM+DRMT
endif
DLOG1=DLOG(dpl/dp3)
DLOG2=DLOG(dp2/dp3)
DDIS=DLOG1+((DZB-DZM)*DLOG2+(DRMT-DRMB))/DH
DUM1(M)=DVDAL(M)*(DX-DVDX(M))*DDIS
DUM2(M)=DVDAL(M)*(DY-DVDY(M))*DDIS
DPART1=DTHREE*( (DZT-DZ)/DH*(DRMT-DRMB)+DRM-DRMT )
DPART2=-DTWO*(DZ-DZM)*DLOG1+DTWO*(DZB-DZM)*(DZT-DZ)/DH*DLOG2
DPART3=DELSQM/DH*( -DH/dpl+(DZ-DZB)/dp3-(DZ-DZT)/dp2 )
DUM3(M)=DVDAL(M)*(DPART1+DPART2+DPART3)
ENDDO
DQX=O.DO
DqY=O.DO
DQZ=O.DO
CDIR$ NOVECTOR
DO M=1,NVDPTS
DQX=DQX+DUM1(M)
DqY=DqY+DUM2(M)
DQZ=DQZ+DUM3(M)
ENDDO
CDIR$ VECTOR
DQX=DK*D¥DBSq*DQX
DqY=DK*DVDBSq*DqY
Dqz=DK*DVDBsq*oqz
RETURN
END
80
-------
REAL*8 FUNCTION DFPOT2HEAD(DPOT,DX,DY,DZ)
IMPLICIT REAL*8 (D)
SAVE
INCLUDE 'vardens.cmn'
WRITE(*,*)'DPOT,DX,DY,DZ ',DPOT,DX,DY,DZ
DK=DFAQK()
DH=DFAQH()
DZB=DFAQZB()
DZT=DFAQZT()
DTOL=l.d-8
DO M=1,NVDPTS
DXMSq=DVDBSq*(DX-DVDX(M))**2
DYMSq=DVDBSq*(DY-DVDY(M))**2
DZM=DVDZ(M)
DZMZM=DZ~DZM
DZMSQ=(DZ-DZM)**2
DZTMSq=(DZT-DZM)**2
DZBMSQ=(DZB-DZM)**2
DELSQM=DVDDELSQ(M)
DRMSQ=(DXMSQ+DYMSq+DZMSQ+DELSQM)
DRMTSQ=(DXMSQ+DYMSQ+DZTMSQ+DELSQM)
DRMBSq=(DXMSq+DYMSq+DZBMSQ+DELSQM)
DRM=DSQRT(DRMSQ)
DRMT=DSqRT(DRMTSq)
DRMB=DSQRT(DRMBSQ)
DPART1=DHALF*(DRMSQ-DZMSQ)*DLOG(DZMZM+DRM)+DHALF*DZMZM*DRM
DPART2=DHALF*(DZT-DZM)*(DRMTSq-DZTMSq)*DLOG(DZT-DZM+DRMT)+
dhalf*DRMT*DZTMSQ-dld3*DRMTSQ*DRMT
DPART3=DHALF*(DZB-DZM)*(DRMBSq-DZBMSq)*DLOG(DZB-DZM+DRMB)+
dhalf*DRMB*DZBMSq-dld3*DRMBSq*DRMB
DUM1(M)=-DPART1+(DPART2-DPART3)/DH
ENDDO
DHEAD=DPOT/(DK*DH)-DVDNUO*(DZ-(DZT**2-DZB**2)/(2.DO*DH))
CDIR$ NOVECTOR
DO M=1,NVDPTS
DHEAD=DHEAD+DVDAL(M)*DUM1(M)
ENDDO
CDIR$ VECTOR
DFPOT2HEAD=DHEAD
RETURN
END
REAL*8 FUNCTION DFHEAD2POT(DHEAD,DX,DY,DZ)
IMPLICIT REAL*8 (D)
SAVE
INCLUDE 'vardens.cmn'
DK=DFAqK()
DH=DFAqH()
DZB=DFAqZB()
DZT=DFAqZT()
DTOL=l.d-8
DO M=1,NVDPTS
DXMSq=DVDBSq*(DX-DVDX(M))**2
DYMSq=DVDBSq*(DY-DVDY(M))**2
DZM=DVDZ(M)
DZMZM=DZ-DZM
81
-------
DZMSQ=(DZ-DZM)**2
DZTMSQ=(DZT-DZM)**2
DZBMSQ=(DZB-DZM)**2
DELSQM=DVDDELSQ(M)
DRMSQ=(DXMSQ+DYMSq+DZMSQ+DELSQM)
DRMTSQ=(DXMSQ+DYMSQ+DZTMSQ+DELSQM)
DRMBSQ=(DXMSQ+DYMSQ+DZBMSQ+DELSQM)
DRM=DSQRT(DRMSQ)
DRMT=DSQRT(DRMTSQ)
DRMB=DSQRT(DRMBSQ)
DPART1=DHALF*(DRMSQ-DZMSQ)*DLOG(DZMZM+DRM)+DHALF*DZMZM*DRM
DPART2=DHALF*(DZT-DZM)*(DRMTSQ-DZTMSQ)*DLOG(DZT-DZM+DRMT)+
dhalf*DRMT*DZTMSQ-dld3*DRMTSQ*DRMT
DPART3=DHALF*(DZB-DZM)=KDRMBSQ-DZBMSQ)*DLOG(DZB-DZM+DRMB)+
dhalf*DRMB*DZBMSQ--dld3*DRMBSQ*DRMB
DUM1(M)=-DPART1+(DPART2-DPART3)/DH
ENDDO
DPOT=DHEAD*DK*DH+DK*DH*DVDNUO*(DZ-(DZT**2-DZB**2)/(2.DO*DH))
CDIR$ NOVECTOR
DO M=1,NVDPTS
DPOT=DPOT-DK*DH*D¥DAL(M)*DUM1(M)
ENDDO
CDIR$ VECTOR
DFHEAD2POT=DPOT
RETURN
END
------- |