Analysis of freshwater lens formation in saline aquifers

MISSING IMAGE

Material Information

Title:
Analysis of freshwater lens formation in saline aquifers
Physical Description:
xvii, 172 leaves : ill. ; 28 cm.
Language:
English
Creator:
Glass, John Patrick, 1946-
Publication Date:

Subjects

Subjects / Keywords:
Aquifers   ( lcsh )
Genre:
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )

Notes

Thesis:
Thesis--University of Florida.
Bibliography:
Includes bibliographical references (leaves 168-171).
Statement of Responsibility:
by John Patrick Glass.
General Note:
Typescript.
General Note:
Vita.

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 000087927
notis - AAK3300
oclc - 05581693
System ID:
AA00003477:00001

Full Text











ANALYSIS OF FRESHWATER LENS FORMATION
IN SALINE AQUIFERS










BY

JOHN PATRICK GLASS


A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE
DEGREE OF DOCTOR OF PHILOSOPHY







UNIVERSITY OF FLORIDA


1979














ACKNOWLEDGEMENTS


This research project was initiated by Dr. B. A. Christensen and

Dr. Hillel Rubin. It was supported financially by the Office of Water

Research and Technology of the U.S. Department of the Interior.

The author extends his gratitude to Dr. B. A. Christensen for his

continued support and advice and for serving as Chairman of his

Supervisory Committee.

Special thanks are due to Dr. Joel Melville whose day-to-day

suggestions and assistance contributed substantially to the work.

Above all, I am thankful for the assistance of my wife, Claudette,

who participated in the experimental work, helped with the preparation

of the figures, typed the entire manuscript and supported me financially

throughout my graduate study.















TABLE OF CONTENTS


ACKNOWLEDGEMENT .


LIST OF TABLES


LIST OF FIGURES . . .

LIST OF SYMBOLS . . .

ABSTRACT . . .

CHAPTER

1 INTRODUCTION . .

2 LITERATURE SURVEY . .

General . .
Studies Dealing with Interface Location
Studies Dealing with Dispersion .
Studies Concerned Directly with Freshwater

3 BASIC EQUATIONS AND FLUID PROPERTIES .

Basic Equations . .
Darcy's Law . .
Continuity . .
Well Flow in Leaky Aquifers .
The Dispersion Equation .
Boundary Conditions at an Interface .
Fluid Properties . .
Density . .
Viscosity . .

4 THE HELE-SHAW MODEL . .

General . .
The Viscous Flow Analogy . .
Model Scaling . .
Model Construction . .
Properties of the Viscous Fluids .
Operation of the Model . .


. . 52


Page

ii


. . . . v


Lenses




. .
* .




* .


Analysis of Test Results












TABLE OF CONTENTS--continued


Page


5 SIMPLE CASES INVOLVING STEADY FLOW . .


General . .
Freshwater Injection in a Coastal Aquifer .
Saltwater Intrusion . .
Injection Into a Coastal Aquifer .
Confined coastal aquifer . .
Phreatic coastal aquifer . .
Axisymmetric Flow in a Leaky Aquifer .


6 A FINITE ELEMENT MODEL FOR THE STEADY FLOW CASE .

General . . .
Description of the Galerkin Finite Element Scheme .
The Fundamental Concept . .
The Element and its Shape Functions . .
Finite Element Formulation of the Governing Equation
Development of the Element Matrices . .
Computer Implementation . .
Program Verification . . .
Freshwater Lens in a Leaky Saline Aquifer .


68
69
69
70
74
80
82

94

94
94
94
95
99
105
117
118
122


7 NUMERICAL SIMULATION OF LENS GROWTH . .. 127


General . .
Description of the INTERCOMP Model . .
Model Verification . . .
Lens Formation in Aquifers With Weak Vertical Flow .
Lens Formation in Aquifers With Strong Vertical Flow
Dimensionless Representation of Lens Growth .
The Influence of Dispersion on Lens Shape .
Changes in Lens Shape During Recovery ...


8 CONCLUSIONS . . ... .... 149


APPENDIX . . . .


. 152


FINITE ELEMENT PROGRAM DOCUMENTATION . ... 152

REFERENCES . . ... . 168


BIOGRAPHICAL SKETCH . . .


. 172














LIST OF TABLES


Table Page

4.1 Summary of Similitude Requirements for Vertical
Hele-Shaw Model With Two Immiscible Fluids ........... 37

6.1 Shape Functions for the 20 Nodes Element ............... 101

6.2 Partial Derivatives of the Shape Functions ............. 108

7.1 Fluid and Aquifer Properties Used With the
INTERCOMP Model ...................................... 133














LIST OF FIGURES


Figure Page

3.1 Definition Diagram for Well Flow in a Leaky Aquifer ...... 15

3.2 Temperature vs. Density for Water with Various Salt
Concentrations ......................................... 22

3.3 Water Density vs. Salt Concentration at Various
Temperatures ......................................... 23

3.4 Viscosity vs. Temperature for Freshwater and Saltwater ... 25

4.1 Schematic of Gap Between Plates ........................ 28

4.2 Glass Plates with Infusion Reservoirs and Seal .......... 40

4.3 Assembled Model on Its Stand ............................ 41

4.4 Back Side of Model with Constant Head Tanks .............. 42

4.5 Infusion Pump Connected to Model ........................ 43

4.6 Temperature vs. Viscosity for the Two Viscous Fluids ..... 45

4.7 Test Run With Horizontal Model. Injection Rate =
20.44 cc/min ......................................... 48

4.8 Test Run With Horizontal Model. Injection Rate =
41.04 cc/min ......................................... 49

4.9 Instability Developing After 50 Seconds. Injection
Rate = 81.72 cc/min ....................................... 50

4.10 Unstable Interface After 3 Minutes. Injection Rate =
81.72 cc/min ....................................... 51

4.11 Interface After 1 Minute. Injection Rate =
20.44 cc/min ......................................... 53

4.12 Interface After 3 Minutes. Injection Rate =
20.44 cc/min ......................................... 54

4.13 Interface After 5 Minutes. Injection Rate =
20.44 cc/min ......................................... 55

4.14 Interface After 7 Minutes. Injection Rate =
20.44 cc/min ......................................... 56









LIST OF FIGURES--continued


Figure Page

4.15 Interface After 9 Minutes. Injection Rate
20.44 cc/min ......................................... 57

4.16 Interface After 11 Minutes. Injection Rate
20.44 cc/min ......................................... 58

4.17 Interface After 10 Seconds. Injection Rate =
8.16 cc/min Model Depth Reduced by half ............... 59

4.18 Interface After 1 Minute. Injection Rate =
8.16 cc/min Model Depth Reduced by Half ............... 60

4.19 Interface After 2 Minutes. Injection Rate =
8.16 cc/min Model Depth Reduced by Half ................ 61

4.20 Interface After 3 Minutes. Injection Rate
8.16 cc/min Model Depth Reduced by Half ............... 62

4.21 Interface After 4 Minutes. Injection Rate
8.16 cc/min Model Depth Reduced by Half ................ 63

4.22 Interface After 6 Minutes. Injection Rate =
8.16 cc/min Model Depth Reduced by Half ................ 64

4.23 Interface After 8 Minutes. Injection Rate =
8.16 cc/min Model Depth Reduced by Half ................ 65

4.24 Comparison of Experimental Results With Equations
for Interface Rotation ..................... ......... 67

5.1 Shape of Saltwater Wedge (Dupuit Parabola) ............... 71

5.2 Injection Into a Confined Coastal Aquifer ................ 72

5.3 Freshwater Lens Produced by Injection Into a Confined
Coastal Aquifer ......................................... 77

5.4 Dimensionless Representation of Interface Toe,
Q' = Q/ = 10.0 ......................................... 78

5.5 Dimensionless Representation of Interface Toe,
Q' = Q/R = 5.0 ........................................ 79

5.6 Injection Into a Phreatic Aquifer ...................... 81

5.7 Steady Freshwater Lens in a Leaky Saline Aquifer ......... 84

5.8 Differential Wedge of Leaky Aquifer ..................... 87









LIST OF FIGURES--continued


Figure Page

5.9 Integral Curves for Equation (5.47) ...................... 90
2
5.10 Q6/ (27Kb ) vs. rt/B for Three Values of a/b ............. 92

5.11 Section Through Steady Freshwater Lens for Two Values
of Leakage, B,: a = 0 ................................ 93

6.1 Three-Dimensional Element Showing Node Numbers ........... 96

6.2 Distribution of Two Quadratic Shape Functions Over a Two-
Dimensional Quadrilateral Element .................... 98

6.3 Local Coordinate System for the Three-Dimensional
Element ........................................... 100

6.4 Sample Points for Gaussian Quadrature of Equation (6.23).. 112

6.5 Sample Points and Weight Factors for Gaussian Quadrature
of Equation (6.25) ..................................... 114

6.6 Sample Points for Gaussian Quadrature of Equation (6.25)
a = 0.57735; Weight Factor = 1 ........................ 116

6.7 Division of Flow Domain Into Finite Elements ............. 120

6.8 Comparison of Solutions for Flow from an Injection
Well in a Bounded Leaky Aquifer (Constant Head
Boundary, s = 0 at r = 200m) ......................... 121

6.9 Flow Domain for Numerical Location of the Interface in
a Leaky Aquifer ........................................ 123

6.10 Comparison of Solutions for a Steady Lens in a Leaky
Aquifer ................................................ 125

6.11 Comparison of Solutions for a Steady Lens in a Leaky
Aquifer ....................................... .. .... 126

7.1 Finite Difference Grid Used With the INTERCOMP Model ..... 132

7.2 Comparison of the Concentration Distributions
Predicted by the INTERCOMP Model and the Approximate
Solution of Hoopes and Harleman ........................ 136

7.3 Lines of Constant Concentration After 40 Days of
Injection ........................................... 137

7.4 Lines of Constant Concentration After 40 Days of
Injection ............................................. 139


viii








LIST OF FIGURES--continued


Figure Page

7.5 Lines of Constant Concentration After 40 Days of
Injection ......................................... 140

7.6 Movement of the 50% Concentration Line With
Continuous Injection ................................. 142

7.7 A Dimensionless Representation of Gravitational
Segregation During Lens Growth ........................ 144

7.8 Comparison of Lens Boundaries After 30 Days of Injection
With Different Dispersivities ........................ 145

7.9 Movement of 50% Concentration Lines During Production
of Freshwater From a Lens ........................... 147

7.10 Movement of 50% Concentration Line During Production
of Freshwater From a Lens ........................... 148













LIST OF SYMBOLS
The following symbols are employed in this dissertation. Dimensions,

where applicable, are also given in terms of length L, mass M, time T

in square brackets.

A = Area [L2]

a = depth of the top of the confined aquifer below sea level, [L]

B = leakage coefficient, [L]

b = thickness of a confined aquifer [L]
b' = thickness of a leaky layer, [L]
bx = body force in the x-direction, [ML/T2
b = body force in the y-direction, [ML/T2]

bz = body force in the z-direction, [ML/T2]

c = fractional concentration of salinity

D = the dispersion tensor, [L2/T]
D = coefficient of longitudinal dispersion, [L2/T]

Dm = coefficient of molecular diffusion, [L2/T]

Dt = coefficient of transverse dispersion, [L2/T]

d = spacing between plates of the Hele-Shaw model, [L]

F = an arbitrary function

f. = the jth term of the forcing matrix, [L3/T]

g = acceleration due to gravity, [L/T ]
h = vertical thickness of the flow region, [L]

I0 = modified zeroth order Bessel function of the first kind

Il = modified first order Bessel function of the first kind








LIST OF SYMBOLS--continued


i = unit vector in the x-direction

[J] = the Jacobian matrix

Jx = x-direction component of head gradient

J = y-direction component of head gradient

Jz = z-direction component of head gradient

j = unit vector in the y-direction

K = hydraulic conductivity for isotropic porous media, [L/T]

K' = hydraulic conductivity of a leaky layer, [L]

KO = modified zeroth order Bessel function of the second kind
K0 = modified first order Bessel function of the second kind

K.. = components of the hydraulic conductivity tensor for
J (i, j = x, y, z), [L/T]
Ki. = terms of the global stiffness matrix (in Chapter 6), [L2/T]

k = intrinsic permeability, [L2]
A
k = unit vector in the z-direction

Kij = terms of the element stiffness matrix (in Chapter 6), [L2/T]

L = a linear operator
Nth
Ni = the it shape function
n = porosity of the porous medium

P = pressure [M/(LT2)]

Pf = pressure on the freshwater side of an interface, [M/(LT2]

Ps = pressure on the saltwater side of an interface, [M/(LT2)]

Q = volumetric rate of injection into a well, [L3/T]

Q' = dimensionless injection rate

Qx = uniform horizontal flow per unit width of aquifer, [L2/T]
q = discharge velocity vector, [L/T]








LIST OF SYMBOLS--continued


qt = longitudinal component of discharge velocity, [L/T]

qx = x-direction component of discharge velocity, [L/T]
qy = y-direction component of discharge velocity, [L/T]

qz = z-direction component of discharge velocity, [L/T]

r = the radial coordinate of a cylindrical coordinate system, [L]

re = radius of a constant head boundary, [L]

rt = radius of the toe of the interface, [L]

S = the boundary surface, [L2]

Ss = specific storage, [1/L]
s = drawdown or pushup, [L]

t = time, [L]

t' = dimensionless time

u = x-direction component of flow velocity, [L/T]

V = volume, [L3]

v = y-direction component of flow velocity, [L/T]

Wn = weighting factor for the nth sample point

w = z-direction component of flow velocity, [L/T]

ws = flow rate through the leaky layer, [L/T]

X = projection of the interface on the horizontal, [L]

Xi = x coordinate of node i, [L]

x = a horizontal coordinate of a cartesian coordinate system, [L]

x' = a dimensionless distance

xt = distance to the toe of the interface, [L]

Xw = distance from the injection well to the sea coast, [L]

x' = dimensionless position of the injection well
S= coordinate of node i, L]
Yi = y coordinate of node i, EL]








LIST OF SYMBOLS-continued


y = a horizontal coordinate of a cartesian coordinate system, [L]

y' = a dimensionless distance

Z = z coordinate of node i, [L]

z = elevation above a datum, or the vertical coordinate, [L]

z. = elevation of a point on the interface, [L]
1
a = coefficient of aquifer compressibility, [LT2/M]

a = a coefficient (in Chapter 5 and 6)

at = longitudinal dispersivity, [L]

at = transverse dispersivity, [L]

6 = coefficient of fluid compressibility, [LT2/M]

S= a constant (in Chapter 5), [L]

2 = a modified leakage coefficient (in Chapter 6) [T]

y = unit weight of a fluid, [M/(LT)2]

6= the Ghyben-Herzberg coefficient

C = the third dimensionless coordinate in an oblique
coordinate system

n = the second dimensionless coordinate in an oblique
coordinate system

0 = angular coordinate in a cylindrical coordinate system

S= dynamic viscosity, [M/(LT)]

v = kinematic viscosity, [L2/T]

S= the first dimensionless coordinate in an oblique
coordinate system


3
p = fluid density, [M/L3]
p = average density, [M/L ]

D = the Girinskii potential (in Chapter 5), [L3/T]

S= nodal values of for n = 1 .20, [L]


xiii









LIST OF SYMBOLS--continued


= piezometric head, [L]

= an approximation to the piezometric head, [L]

(0 = undisturbed piezometric head in a static aquifer, [L]

f = piezometric head of freshwater measured in terms of fresh-
water, [L]

ps = piezometric head of saltwater measured in terms of saltwater,
[L]


Subscripts

f -

i -

m -

p -

r -

s -

t -

x -

y -

z -


and Superscripts

freshwater or lighter fluid

on the interface

model

prototype

ratio

saltwater or heavier fluid

at the toe of the interface

x-direction

y-direction

z-direction














Abstract of Dissertation Presented to the Graduate Council
of the University of Florida in Partial Fulfillment
of the Requirements for the Degree of
Doctor of Philosophy



ANALYSIS OF FRESHWATER LENS FORMATION IN
SALINE AQUIFERS


By


John Patrick Glass

June, 1979

Chairman: B. A. Christensen
Major Department: Civil Engineering


The injection of freshwater into a saline aquifer can, under

favorable conditions, result in displacement of the saltwater and crea-

tion of a distinct body of freshwater called a freshwater lens. It may

be desirable to do this either as a method of storing the freshwater

for future use or in order to form a barrier against saltwater intrusion.

In either case, it is important to know beforehand whether the lens will

be stable or the injected water will be lost. The purpose of this

study is to develop methods for analyzing the hydrodynamic factors

involved in the formation and maintenance of such lenses.

A series of laboratory tests using vertical Hele-Shaw models was

performed to determine the shape and rate of movement of the initially

vertical interface between two fluids of different density. The

results of these tests were compared with the predictions of two formulas









proposed by previous investigators. It was found that for relatively

thick aquifers it is not correct to neglect the vertical components of

flow as has formerly been done.

An analytical method is presented for the location of the boundary

of a steady freshwater lens in a coastal aquifer. A vertically

integrated potential function is used to superimpose the flow from the

injected well on the pre-existing uniform seaward flow in the aquifer.

This results in a quasi-three dimensional solution for the position of

the coastal interface.

A pair of coupled ordinary differential equations is derived which

describe the steady flow of freshwater injected into a leaky saline

aquifer under the Dupuit assumption. The non-linear system of equations

is solved numerically by isolating the unique integral curve which

results in the only physically possible boundary conditions. A design

chart is presented by which the maximum possible size of freshwater

lens can be determined for various aquifer conditions.

A three dimensional finite element model is developed which can be

used to locate the boundaries of a steady freshwater lens under a wide

variety of conditions. The model uses isoparametric elements of the

quadrilateral family having twenty nodes each. This permits very close

approximation to the curved flow boundaries typical of freshwater lenses.

Inhomogeneous and anisotropic aquifers can be treated.

A finite difference model developed for the U.S. Geological Survey

is used to analyze the axisymmetric flow in transient freshwater lenses.

Simulation of lens formation in a variety of saline aquifer configura-

tions indicates that previously available methods of analysis which

neglect vertical components of flow may be reliable in very thin








aquifers. For thicker aquifers, however, vertical flow can predominate

resulting in a much more rapid rate of gravitational segregation of

the fluids.


xvii














CHAPTER 1

INTRODUCTION

Groundwater flow is an important part of the hydrologic cycle which

nature uses to even out the temporal and spatial variations inherent in

rainfall. It has long been man's practice to take water from groundwater

reservoirs with the assumption that the water removed will be replaced by

natural processes. In many cases, this arrangement works well, with

nature putting water into the ground and man taking it out. In some

areas, though, the increased use of groundwater supplies coupled with

decreased natural recharge due to man-made changes in the land surface

is either draining the freshwater aquifers or opening them to saltwater

intrusion. In recognition of this, efforts have been made to augment

natural recharge by such means as spray irrigation, water spreading,

seepage ponds and deep well injection. This signals a new conception

of aquifers as a form of underground storage which can be used in much

the same way as surface reservoirs. As an extension of this, it has been

proposed to make use of the storage capacity available in saline aquifers

by injecting freshwater into them in such a way that the resident salt-

water is displaced and a lens of freshwater is created.

The variety of possible flow patterns that may result from the

injection of freshwater into a saline aquifer can be viewed as a contin-

uous spectrum bounded by three idealized extremes. One possibility is

that the freshwater will disperse rapidly and lose its identity as a

separate flow, merely serving to dilute the resident saltwater. In





2


another case, a flow pattern may develop which transports the freshwater

away from the point of injection as a plume. A third alternative is

that the freshwater will tend to remain near the point of injection

creating an identifiable lens of freshwater which grows or diminishes

accordingly as fluid is added to or removed from it. If the freshwater

is injected in anticipation of its subsequent retrieval and re-use,

the formation and maintenance of such a lens adjacent to the point of

injection is essential.

The major factors influencing the formation of freshwater lenses

are: (1) buoyancy, (2) mode of injection, (3) the properties of the

porous medium, (4) dispersion, (5) flow of the resident saltwater.

The relative intensities of these'individual influences determine which

of the previously described flow patterns will predominate. Favorable

conditions for the development of freshwater lenses require that

buoyancy, dispersion, and preexisting flow be relatively weak influences.

In previous studies of freshwater lenses such ideal conditions have

been assumed. In reality, these effects are never entirely absent.

It is important, then, to determine at what point their effects become

substantial and to have techniques for including their influences in

the analysis of freshwater lenses.














CHAPTER 2

LITERATURE SURVEY


General

The analysis of freshwater lenses is one of a general class of

problems involving the stratified flow of two miscible liquids through

porous media. It is closely related to such other problems of that

class as saltwater intrusion, interface upcoming beneath wells, and

deep well disposal of liquid wastes. A wealth of published information

is available on each of these topics but much of it has no direct

application to the study of freshwater lenses. Many of the methods

developed for locating coastal interfaces, for instance, depend on

conformal mapping techniques or otherwise on the theory of complex

variables, which is limited to two dimensional problems. Likewise,

much of the work done with deep well disposal deals with the movement

and dispersion of plumes, emphasizing the rapid departure of the waste

fluid from the scene of injection. Nonetheless, some of the work in

each of these fields can be useful in lens analysis. It is only these

studies and those which are directly concerned with fresh or hot water

lenses that will be mentioned here.


Studies Dealing with Interface Location

The rotational movement of an initially vertical interface between

two fluids of different density was studied by-Gardner, Downie and

Kendall (1962) because of its applicability to recovery processes in

petroleum reservoirs. The initial motion of the interface and the









long term motion were found analytically using the assumption that the

porous medium was thin and vertical flow could be neglected. These two

solutions were then joined to give one equation which was assumed

accurate for intermediate times as well. It was assumed in these

derivations that the interface remained straight during its rotation

toward the horizontal. Laboratory tests using models of square and

circular cross-section packed with glass beads were performed which

seemed to agree with the equation.

Dagan and Bear (1967) treated the upward movement of a freshwater-

saltwater interface beneath a partially penetrating well using a pertur-

bation expansion to linearize the interface boundary condition. They

checked their results with a sand-packed model and found them to be

reasonable when the disturbance of the interface was small. The range

of applicability of this method, however, appears to be quite limited.

Shamir and Dagan (1971) developed a transient numerical model to

predict the landward movement of a saltwater wedge in a coastal aquifer.

The Dupuit assumption was used to formulate the governing equation which

was then solved by an implicit finite difference method. Dispersion was

neglected. The numerical results were checked against an analytical

solution similar to that of Gardner, Downie and Kendall and with a

Hele-Shaw model. The comparison with the analytical solution showed

good agreement that the interface moved as a straight line. Comparison

with the Hele-Shaw model, however, was not quite as good. The authors

attributed this to the fact that vertical velocities were neglected in

both the analytical and numerical solutions.

Strack (1976) studied the influence of inland production wells on

the steady state shape of the saltwater wedge in coastal aquifers. His









analysis makes use of a vertically integrated potential function which

can be defined so as to be single-valued, throughout the aquifer.

Lines of constant potential can then be found in the horizontal plane

which lead to a three-dimensional description of the interface based on

the Ghyben-Herzberg principle.

A similar problem was studied by Kishi and Fukuo (1977). They

considered the effect of inland pumping on the saltwater wedge in a fan

shaped coastal aquifer. The Laplace equation in radial coordinates is

integrated vertically to make it two dimensional. Using the Ghyben-

Herzberg principle the dependent variable is transformed so that it will

be single-valued in both the freshwater region and the interface region.

The resulting equation is solved by the method of Green's functions to

give a three-dimensional description of the interface.


Studies Dealing with Dispersion

Gardner, Downie and Wyllie (1962) considered the mixing across the

interface between methane gas and nitrogen used as a cushion gas in a

natural gas reservoir. A solution to the convective dispersion equation

in radial coordinates was combined with the gravity segregation formula

found by Gardner, Downie and Kendall.

Harleman and Rumer (1962) derived an equation for interface rota-

tion which is the same as that found by Gardner, Downie and Kendall for

the latter stages of movement. This equation also is based on the

assumption of a straight interface. The same report also deals with

longitudinal and lateral dispersion in both stationary and moving

coastal interfaces. Laboratory tests using sand-packed models to

simulate a coastal aquifer showed that in the absence of tidal fluctua-

tions the mixing zone was quite thin, especially near the toe of the









interface. An equation was presented which predicted the thickness of

the dispersion zone accurately in this region. When changing conditions

in either the saltwater or the freshwater caused a shift in the equi-

librium position of the interface the effects of longitudinal dispersion

became pronounced. It was concluded that such shifts are responsible

for most of the dispersion found in natural coastal aquifers.

Wooding (1963) obtained approximate solutions for the lateral

dispersion across a stationary interface by transforming the governing

equations into a natural coordinate system based on the shape of the

interface. Using order of magnitude reasoning these equations were

reduced to the form of Prandtl's boundary layer equations. For these,

the usual similarity solution was obtained which relates the thickness

of the mixing layer to the distance from the stagnation point and the

slope of the interface.

Hoopes and Harleman (1965 and 1967) analyzed the dispersion of

tracers injected into groundwater. Their work consisted of analytical

and numerical solutions to the convective dispersion equation and

laboratory experiments with a sand-packed model to verify these solu-

tions. The relative effects of convection, lateral and longitudinal

dispersion, and molecular diffusion were compared for the case of

axisymmetric radial flow from an injection well and for the flow between

an injection well and production well. It was shown that the effects

of lateral dispersion on the steady flow between two wells was almost

negligible except near the line joining the two wells.


Studies Concerned Directly with Freshwater Lenses

Glazunov (1967) proposed to use a cluster of three to five wells

operated in various patterns of injection and production in order to









form freshwater lenses of the desired shape. His analysis neglects

both dispersion and gravity segregation, treating the convection of an

interface in two horizontal directions by a method introduced by Muskat

(1946). In this method the interface is represented as a material

surface which is transported along the streamlines by the flow while

exerting no effect on it as a boundary.

Esmail and Kimbler (1967) studied the effects of radial flow

geometry on gravity segregation and longitudinal dispersion in labora-

tory models composed of glued sand grains. The solution to the convec-

tive dispersion equation given by Gardner, Downie and Wyllie was verified

for cyclic storage in freshwater lenses for several cycles of injection

and production in the model. A dimensionless parameter analysis of the

influence of radial geometry and interfacial mixing on the rate of

gravity segregation was conducted resulting in an empirical equation

for the rate of interface rotation. It was found that there is inter-

action between the processes of dispersion and interface rotation which

retards gravity segregation. The authors concluded that the horizontal

projection of the interface is proportional to the square root of the

density gradient across it and that gravity segregation in their model

had more effect on the recovery efficiency than did dispersion. Under

these conditions, interface dispersion may actually improve the recovery

efficiency by retarding gravity segregation.

Kumar and Kimbler (1970) continued the research with glued sand

models, finding that the recovery efficiency improves after several

cycles of injection and recovery.

Whitehead (1974) developed a computational method for predicting

the recovery efficiency of a single well cyclic storage system in thin









horizontal saline aquifers. His method was based on the empirical

equation of Esmail and Kimbler. Whitehead verified the results of

his computational scheme using the glued sand grain model and found

them to be accurate to within 10% for three storage cycles. The

recovery efficiencies achieved in the laboratory tests were in the

range of 70-100%. It was concluded that, given favorable subsurface

conditions, the use of saline aquifers for freshwater storage could be

economically feasible.

In a series of reports (Agrawal, 1975; D'Amico, 1975; Tate, 1976;

Kimbler and Whitehead, 1977) the results of Whitehead were extended to

include the effects of aquifer dip and viscosity difference between the

resident and injected fluids. These studies also used glued sand

models. The conclusions reached from these experiments were that the

best conditions for cyclic storage occur when the saline aquifer is

horizontal and the injected fluid has the same viscosity as the resident

fluid.

It is widely recognized that a pre-existing uniform flow in the

saline aquifer is detrimental to the recovery efficiency of a lens

storage system. Langhtee (1974) proposed to surround the storage area

with bounding wells which would be pumped in such a way that a stag-

nation zone would result inside the boundary. The required pumping

rates for these wells were found by superimposing the streamline and

equipotential pattern for the individual wells and solving for the

pumping rates which would give the best least-squares fit to the desired

piezometric head on the boundary.

Molz and Bell (1977) also considered the use of boundary wells for

head gradient control. They used a linear programming technique to find










the optimum set of bounding wells which would give the desired stagna-

tion in the storage area with the least expenditure of energy for

pumping.

The shape and movement of naturally occurring freshwater lenses on

Grand Cayman Island were studied by Childley and Lloyd (1977). They

used a modified Ghyben-Herzberg ratio based on field measurements of

salinity distribution in the mixing zone of the lens. The long term

transient effects of seasonal variations in recharge rate were modeled

numerically using the idea of a succession of steady states. Comparison

of the numerical results with records of actual groundwater conditions

on the island indicated that their model was useful in predicting long-

term changes due to water supply wells and changes in recharge condi-

tions. They concluded that, while it may take many years for a pumped

lens to reach equilibrium, an overpumped lens can collapse locally in

a matter of a few months.

Smith and Hanor (1975) report on a field test of the freshwater

lens storage concept. A 10 centimeter fully penetrating well was

drilled into a confined saline aquifer 36 meters thick to be used for

injection and retrieval of freshwater. After the construction of this

well it was discovered that a uniform flow of approximately 15 centi-

meters per day existed in the aquifer. Recognizing the unfavorable

conditions, the investigators proceeded with two tests, each using about

757,000 liters of freshwater. In the first test, the water was stored

for two hours and 41% was recovered. In the second test the storage

period was six days and the recovery efficiency was 24.5%.

Werner and Kley (1977) conducted a field test of the storage of

heated water in a water table aquifer 3.35 meters thick. A total of









430 cubic meters of water at a temperature of 45C was injected and the

movement of the heat was monitored with buried temperature sensors for

a period of 64 days. The heated water remained close to the injection

well throughout the test period but no attempt was made to recover it.

Molz, Warman and Jones (1978) reported a field test with a

partially penetrating well in a confined aquifer 21 meters thick. A

group of 14 observation wells were constructed to monitor the movement

of the hot water injected. In this test, storage of 7570 cubic meters

of hot water for 36 days resulted in an energy recovery factor of 0.69

which was considered promising.

Papadopulos and Larson (1978) applied a finite difference numerical

model developed for the U.S. Geological Survey (INTERCOMP, 1976) to

simulate the test results of Molz, Warman and Jones. The predicted

recovery efficiency was 75% versus a value of 69% measured in the field.

Considering the accuracy of the measured aquifer properties, this

prediction was considered quite good.













CHAPTER 3

BASIC EQUATIONS AND FLUID PROPERTIES


Basic Equations

Darcy's Law

The groundwater flow to be dealt with in this report is limited

to laminar flow in porous media with only primary porosity. The

equation of motion for this flow is Darcy's law which can be expressed,

in the cartesian coordinate system, as


qx = Kxx x + Kxy Jy + Kxz J


q = x J + Ky + K J


qz = Kzx J + Ky Jy + Kz J


q X9 qy, q


(3.la)


(3.1b)


(3.1c)


= components of the discharge velocity in the

x, y, and z directions, respectively

= components of the hydraulic conductivity tensor

for (i, j = x, y, z)


ax

ay


where









z z
zz

S = piezometric head

If it is assumed that the coordinate axes are parallel to the principal

directions of hydraulic conductivity then Darcy's law can be simplified

to


= K K j K l k (3.2)
xx ax yy ay zz a


where i, j, and k are unit vectors in the x, y, and z directions,

respectively. This assumption is not a major restriction to the use of

equation (3.2) because in practice it is usually difficult to distin-

guish any directional hydraulic conductivities except the vertical and

the horizontal. If the soil is isotropic or only one dimensional flow

is being considered equation (3.2) reduces to


q =-K V) (3.3)

The value of the hydraulic conductivity, K, or of its directional

components depends on the properties of both the soil and the fluid.

Since the properties of saltwater and freshwater are somewhat different

it may be necessary in the general case to consider the hydraulic

conductivity as a function of the salt concentration. The flow resis-

tance of the soil is then described in terms of its intrinsic permea-

bility, k, which is related to hydraulic conductivity by


K = kg/v (3.4)

where v is the kinematic viscosity of the fluid and g is the accelera-

tion due to gravity.









The piezometric head, q, appearing in equations (3.2) and (3.3) is

defined as


S= z + P/y (3.5)

which depends on the unit weight, y, of the fluid concerned. Since y is

also a function of salt concentration it is often more convenient to

deal directly with pressure, P, in the equation of motion.

Continuity

The continuity equation for groundwater flow is


V q =- (3.6)


in which Ss, the specific storage, is defined as


Ss = pg(a + nB) (3.7)

where

a = coefficient of soil compressibility

B = coefficient of compressibility for water

n = porosity

p = fluid density

The combination of equations (3.2) and (3.6) gives the governing

equation for groundwater flow


Kxx + K I + K S (3.8)
ax xx x ay 3yy y az zza z s t
.x j i /Z[ S 38








For steady flow in homogeneous and isotropic soil this reduces to the

Laplace equation



+= 0 (3.9)
ax 2 y az

Well Flow in Leaky Aquifers

In dealing with flow to and from wells, it is customary to denote

the head gradient in terms of drawdown, s, which is defined as


s = ( Qo (3.10)


where (0 is the head above the leaky layer (see Figure 3.1). The flow

from an injection well causes an increase of head, in which case s will

be called "push up" instead of drawdown. A leaky aquifer is one which

is confined by a layer of material which has low hydraulic conductivity,

but not so low that it can be considered impervious. The governing

equation for flow in a leaky aquifer is


V2s s/B2 s as (3.11)
Kb at


The leakage coefficient, B, is given by


B2 = Kbb'/K' (3.12)

where

K' = hydraulic conductivity of the leaky layer

b' = thickness of the leaky layer

b = thickness of the aquifer











The Dispersion Equation

The mixing of saltwater with freshwater at the periphery of a

freshwater lens is described by the convective dispersion equation



ac 1
+t- nq Vc = V D Vc (3.13)

in which c is the concentration of the fluid being followed and D is

the dispersion tensor. For isotropic porous media this tensor has only

two independent terms, D. and Dt, corresponding to longitudinal and

transverse dispersion respectively, relative to the direction of the

mean flow. The values of these coefficients depend not only on the

properties of the porous medium but also on the flow velocity and the

properties of the fluid involved. It is often assumed that these

coefficients can be represented by


= Dm + a q /n (3.14a)


Dt = Dm + t qq/n (3.14b)

where


Dm = coefficient of molecular diffusion in the porous medium

aC = longitudinal dispersivity

at = transverse dispersivity

qp = discharge velocity in the direction of mean flow


The dispersivity coefficients a, and at are properties of the

porous medium only. Their magnitudes are related to the intrinsic

permeability and the pore size distribution, but the actual values must








be found experimentally. The longitudinal dispersion coefficient, D,,

is usually 10 to 20 times greater than Dt but the ratio of the two

depends on the flow velocity, being close to unity when the flow is

very slow and the influence of molecular diffusion predominates.

For axisymmetric radial flow of a tracer-like substance from an

injection well, equation (3.13) can be simplified to


ac ac- nr
3c c = Q ) a c (3.15)
it 2nbnr ar X 2abnr @2 V

because there are no components of either velocity or concentration

gradient except in the radial direction. This equation has been

solved numerically for the case of steady flow from the well and con-

tinuous tracer injection beginning instantanously (Hoopes and Harleman,

1965; and Shamir and Harleman, 1966).

Freshwater injected into a saline aquifer does not behave as a

tracer because of the difference in density, so equation (3.15) does

not apply even though the flow may be asixymmetric. Because vertical

components of both velocity and concentration gradient are present,

equation (3.13) must be used. The numerical solution of this equation

will be the subject of Chapter 7 of this report.


Boundary Conditions at an Interface

When the mixing zone separating saltwater and freshwater is

relatively thin, it is often convenient to neglect dispersion and

assume that the two fluids are separated by a sharp interface. If

such a boundary between the two fluids exists, a separate governing

equation can be written for each zone and the two flows will be coupled

by their common boundary condition at the interface.








The interface is a material surface which is always composed of

the same fluid particles. The essential features of such a surface

are that both pressure and the velocity component normal to the

interface must be continuous across it. For a stationary interface

the normal velocity component is zero so the boundary condition is

concerned with pressure only. The pressure on the freshwater side of

the interface is


Pf = Yf (f zi) (3.16a)


and on the saltwater side


PS = YS (s zi) (3.16b)


where zi is the elevation of a point on the interface and the sub-

scripts f and s stand for freshwater and saltwater, respectively.

Requiring that these pressures be equal results in an expression for

the boundary condition at a stationary interface

Ys Yf
z. Y Y f (3.17)
i AY s AY f

where


A = Ys Yf


For the special case where the saltwater is at rest and a hydro-

static pressure distribution prevails in the freshwater zone this

reduces to the Ghyben-Herzberg relationship

Yf
Az =- AAf (3.18a)
Ay





19


or

Az = 6APf (3.18b)

where 6 = Yf/AY is the Ghyben-Herzberg ratio. This states that a change

in the piezometric head of the freshwater causes a change in the eleva-

tion of the interface that is opposite in sign and greater in magnitude

by a factor of 6.

If the interface is not stationary the condition of continuity of

the normal velocity component must also be used. The interfacial surface

can then be described by some function, F, such that


F (x, y, z, t) = 0 (3.19)

This surface is convected by the flow so that


F+ q vF = 0 (3.20)
at n

where

q = discharge velocity of interfacial particles

n = porosity of the aquifer

In terms of the velocities adjacent to the interface on both sides of it

this becomes


n | + q- VF = 0 (3.21a)

for the freshwater, and


aF 0
n -F + q VF = 0 (3.21b)
at s








for the saltwater. In terms of the elevation of a point on the inter-

face, zi, the function, F, can be written as


F(x, y, z, t) = z zi (x, y, t) = 0 (3.22)

Substituting this, together with Darcy's law, into equations (3.21) gives


az.
-n KVqf v(z z.) = 0 (3.23a)
at
3Z.
-n Kv v(z z ) = 0 (3.23b)
at 1


Performing the vector operations indicated in these equations results

in the boundary conditions for the freshwater zone and the saltwater

zone. The interface boundary condition for the freshwater zone, then,

is


Yf a f Ys 4s Yf 2 Y
nn K (V~ f) + K (f s)
AY at y at Ay f Ay f



K 0 (3.24a)
az


and for the saltwater zone,


Yf 4, Ys s 2 Yf
nf f s + K (s ) + K K (vqf )
AY-- at Ay at a- S Ay f


K s -= 0 (3.24b)
az









Fluid Properties

Density

Up to this point the fluids with which this report is concerned

have been referred to simply as freshwater and saltwater. These conven-

ient terms are used only in a relative sense and the actual properties

of the fluids deserve some further elaboration. In the ensuing analy-

sis the primary difference between the resident and injected fluids

will be their density. In general, this is a property that depends

on pressure, temperature and salt concentration. In this report,

however, the change in density due to compressibility will be neglected

because it is so small, and the thermal effects, while not necessarily

small, will be eliminated by assuming that the temperature of the

injected fluid is the same as that of the resident saltwater.

Figure 3.2 shows the relationship of density to temperature for

water with several different salt concentrations. With the exception

of the 3.5% concentration line, this figure is based on data for

solutions of pure sodium chloride in water. The 3.5% concentration

line shows the temperature versus density relation for seawater having

a dissolved solids concentration of 3.5%. Saline groundwater is

commonly associated with the intrusion of seawater into coastal aquifers.

However, it is also found in many locations which have no apparent

geological association with the oceans and where its concentration may

be much greater or less than that of seawater.

Figure 3.3 shows the relationship of density to sodium chloride

concentration. It is worthwhile to note that in this range the

relationship is apparently linear.















1.04





1.03




3. 5
1.02

E


S 1.01
I'--"

=7%


1.00 0%





0.99





0.98
0 10 20 30 40
TEMPERATURE (DEGREES C.)


FIGURE 3.2. TEMPERATURE VS. DENSITY FOR WATER
WITH VARIOUS SALT CONCENTRATIONS









































































































DP C?


23




































U-1




-r -,
c:a
ru cr












-___ ___ __----- ----- HL





LOU
~-- --- -r)I j














-
--- ---- i-'-

























UU
_,W/6)____SN30 a)











LU-
--`- ~ --- ~ ~ --LU





(~-w/6.-l) __
AIISN-O D LU-

---- --------I--~-----IL.-


NT-~-j--


" `








Viscosity

It was shown in equation (3.4) that the hydraulic conductivity

of an aquifer is partially determined by the kinematic viscosity of its

water. Viscosity is a function of both temperature and salt concentra-

tion but within the range of present interest temperature is the more

important factor of the two. The temperature dependence of kinematic

viscosity for both pure water and seawater is shown in Figure 3.4. It

is interesting to note that the kinematic viscosity of pure sodium

chloride solutions, as given by Kaufmann (1960) for concentrations up

to 4%, differ so little from the kinematic viscosity of pure water that

they could not be distinguished from it if plotted at the scale of

Figure 3.4. Since differences in temperature are not being considered

in this study and the influence of salt concentration on viscosity is

relatively weak, viscosity differences will be neglected also.












1.4





1.3





1.2

in
o
I--




>-









FRESHWATER






0.8





0.7
S10 20 30 4

TEMPERATURE (DEGREES C.)

FIGURE 3.4. VISCOSITY VS. TEMPERATURE FOR
0 gg FRESHWATER






0.8





0.7
0 10 20 30 40

TEMPERATURE (DEGREES C.)

FIGURE 3.4. VISCOSITY VS. TEMPERATURE FOR


FRESHWATER AND SALTWATER














CHAPTER 4

THE HELE-SHAW MODEL


General

Valuable insights into the flow phenomena associated with the

injection of fresh water into saline groundwater can be obtained by

studying certain simplified aspects of the problem. This chapter will

describe the use of a physical model to analyze the movement of the

interface between the two fluids in the absence of dispersion. The

Hele-Shaw analog was chosen for this work because it is simple and

inexpensive to construct and the test results are easily observable

without elaborate instrumentation. The Hele-Shaw analog consists of

one or more viscous liquids flowing in the gap between two closely

spaced flat plates. Its usefulness in groundwater studies is based

on the fact that the Navier-Stokes equations governing the viscous

flow in this gap can be reduced to a form similar to Darcy's law in

two dimensions.


The Viscous Flow Analogy

If u, v, and w are the velocity components in the x, y, and z

directions, respectively then the Navier-Stokes equations for incom-

pressible flow in a cartesian coordinate system are


u + u + v u w b + vV 2u (4.1a)
at x + y+ az x p ax

+ u @v + v + w b + vv (4.1b)
at ax ay az y P y









aw aw cw aw ) aP 2
w+ + v + w W b P + v2w (4.1c)
t ax 3y -z z p Tz

where

bx, b and bz = components of acceleration due to body forces in

the x, y, and z directions, respectively

p = density of the fluid

v = kinematic viscosity of the fluid

P = pressure

Applying these equations to the flow between the plates shown in

Figure 4.1, it can be assumed that because of the low Reynolds number

(Re < 1 for this model) the flow is laminar and there is no velocity

component, v, in the direction normal to the plates. Note also that,

since the velocity at the surface of the plates is zero, the velocity

gradients normal to the plates are very large compared to those in the

x and z directions so that the latter can be neglected. Finally, the

only body force acting on the fluid is gravity, so bz = g and

bx = b = 0. Incorporating these observations into the Navier-Stokes
x y
equations and multiplying by p results in


aP a2u
=x 2 (4.2a)
ay

DP
0 (4.2b)


.P 82u
Y + z 2U (4.2c)
Y 3y 2
"y

where y = pg is the unit weight of the fluid and p = pv is its dynamic

viscosity. Recalling that the piezometric head, 4, is defined as





28







x













U-1
ui


i-


u-i
1=
a-
SI-













0 /01
Iz:
LU
U-
'=4




N -__




/ 3
>1w








{ = z + P/y (4.3)

it follows that equations (4.2) when expressed in terms of q and divided

by y become


ac p a2u (4.4a)
ax y 2


S= 0 (4.4b)
ay


= a2u (4.4c)
az y 2

From the second of these equations it is seen that at any point (x, z)

in the plane of the model, the piezometric head is not a function of y.

Therefore, the first and third of these equations can be integrated

with respect to y giving


y au + C (4.5a)
ax Y ay


y = 2w + C (4.5b)
az y ay

The constant of integration, C, is found to be zero by observing that

at the center of the flow region, where y = 0 both and are also

zero, so that


y = u (4.6a)
ax y ay


y = w (4.6b)
az y Dy


Integrating these equations again with respect to y gives









Y- i- U + C1
2 Dx Y 1


2
2 az
2 az


(4.7a)


(4.7b)


p w + C1
Y 1


Here, the constant of integration, C1, is found by requiring that, at

the surface of the plates, where y = d/2, u = w = 0. With these

boundary conditions the expressions for the velocity distribution

across the gap are found to be


u = (y2 d2) 2

S(y 2
w = (y2 d2 )
2p 4 az


where d is the thickness of the

to find a velocity analogous to

through porous media, equations

thickness of the gap as follows


q= d/2

-d/2
-d/2

d/2

-d/2


gap between the plates. Now, in order

Darcy's discharge velocity for flow

(4.8) must be averaged across the


dy d Y 4
y 12 P ax




dy = 12
12 P az


(4.9a)




(4.9b)


If the constant coefficients in equations (4.9) are used to define an

effective hydraulic conductivity, Km, for the model


Km 12


(4.8a)


(4.8b)








then the governing equations for flow between the plates become


q Km a
S m- ax

qz = Km -
qz m Dz


in which the resemblance to the

obvious. In further analogy to

ability for the Hele-Shaw model,


so that


two dimensional form of Darcy's law is

groundwater flow, an intrinsic permea-

km, can be defined as


km = d2/12




Km = km


(4.12)




(4.13)


The continuity equation for the flow between the plates is


au + w -_
ax az


(4.14)


since there is no velocity in the y direction. Averaging these veloc-

ities across the gap, as before, results in


(4.15)


ax (Kmax + (Km )= 0


which has the same form as the governing equation for two dimensional

flow in an incompressible porous medium. The model hydraulic conductiv-

ity, Km, is a constant depending only on the spacing, d, of the plates

and the unit weight and viscosity of the fluid. Consequently, for

parallel plates and homogeneous fluids the flow is homogeneous and

isotropic, so that equation (4.15) reduces to the Laplace equation in


(4.11a)


(4.11b)








two dimensions



S + 2 = 0 (4.16)
ax az2


Since Km is a function of physical parameters which are relatively

easy to measure and control, the Hele-Shaw model is a very convenient

tool for the scale modelling of two dimensional groundwater flow.


Model Scaling

In the preceding section it was shown that two dimensional ground-

water flow can be modeled by the Hele-Shaw analog because the governing

equations for the two cases are of the same form. If quantitative

results are to be obtained from the model, however, the scaling ratios

between model and prototype must be chosen so that the governing

equations are actually identical. In order to achieve this, the model

laws which fit the Hele-Shaw model to the prototype case of groundwater

flow must be developed.

Since the present case involves the modeling of two immiscible

liquids separated by an interface the equations on which the model

laws depend are Darcy's law, the continuity equation and the boundary

condition for interface flow. Letting the subscript m denote model

parameters and p denote those in the prototype, Darcy's law for the

model is

i
q m= -Km (4.17a)


1 -m a7b)
qzm -Km (4.17b)









where the superscript


: I :


for the lighter fluid

for the heavier fluid


The corresponding equations for the prototype are


q =
xp


ai1
q = -K -
zp p azp


Let the scaling ratios between

be denoted by the subscript r,


model parameters and prototype parameters

so that


Kr = Km/Kp


(4.19a)


(4.19b)


xr = x/xp


(4.19c)


etc.

Then, substituting these relationships into equations (4.18) gives


K X41
S i Km r m
q xm/qxr K ri X


qm m m
zm zr Kr 3azr
r r r


If these are to be identical to equations (4.17) then it must be

ture that Ki
qi r
xr Xr


(4.20a)



(4.20b)


(4.21a)


-K -
p p


(4.18a)



(4.18b)









K 1
qr rr (4.21b)
zr zr


Equations (4.21) provide one of the model laws which allow results

from the model to be transferred to the prototype.

The continuity equations for isotropic flow in the model and

prototype are

2i 2i
+ =0 (4.22)
Xm m

and

2 i 2i
S+ = 0 (4.23)
ax 2z
p p

Substituting the appropriate parameter ratios into equation (4.23) and

requiring that it be identical to equation (4.22) results in


xr = zr (4.24)


The third similitude requirement is that the boundary conditions

at the moving interface be given by identical equations in both model

and prototype. These interface equations were given in Chapter 3. For

the model they are


Sf s s f 2 s
n m mm n m m K (V ) + K (Vm .Vs) -
mm Y t m y t m Ay m m Am m mn
am Ym m m m



m m = 0 (4.25a)
K zm








and

f f s s s 2 f
m m Ym mm Ym m f s
n m n ( K Vm) + K (Vm Vs) -
Aym atm m ay t m Ay m

s
Kmiz 1 =0 (4.25b)


The corresponding equations for the prototype are the same as equations

(4.25) except that the m subscripts are replaced by p. When the param-

eter ratios are substituted into the prototype equations and equality

is required with the model equations, 16 groupingsof parameter ratios

are generated, which must all be equal to unity. Many of these group-

ings are repetitive and, since they provide no new information, they

will not be listed here. The groupings which are useful are


1 AY, t 1 AY t 1 Ay xr2 1 z 1 z
r r r r r r r r
f--2 s- (4.26)
nr Yr r nr Yr r K Yr r Kr r Kr r

The equality of the first and third terms of this equation gives the

time scale for the model


nr xr
tr r (4.27)
K 4

Equality of the first, second, and fourth terms of equation (4.26)

together with equation (4.27) results in


ff s A
yr = r r = AYr zr (4.28)








Finally, the equality between the last two groupings in equation (4.26)

requires that


f s
=s r (4.29)

which, in turn, means that


f s
Y = Y, (4.30)


The complete list of similitude requirements for matching the

prototype to the model is given in Table (4.1).


Model Construction

The design of the Hele-Shaw model involved the following con-

sideration:

1. The plates must be transparent, relatively inflexible and

demonstrably flat.

2. Provision must be made for the controlled introduction of the

viscous fluids and the release of entrapped air.

3. The model must be easily disassembled for cleaning between

tests.

4. It must be possible to rotate the model from horizontal to

vertical so that the initial shape of the interface can be

controlled.

The availability of materials and the requirement for transparency

dictated that the material to be used for the flat plates be either

1/4 inch plate glass or 1/2 plexiglass. Since 1/2 plexiglass is more

than 2.5 times more flexible than 1/4 inch plate glass, it was decided

to use the latter.














TABLE 4.1

SUMMARY OF SIMILITUDE REQUIREMENTS FOR VERTICAL
HELE-SHAW MODEL WITH TWO IMMISCIBLE FLUIDS


Velocity Ratios:


i r ri
q -
xr Xr



fh
where i =
S,


and


1
Kz r
q
zr r


for the lighter fluid

for the heavier fluid


Geometric Ratios:


x =z
r r


Time Ratio:
2
nr Xr
t
r f
Kr r


Head and Density Ratios:


ff s s
Y, r r r dr

and


s
Y= r








It was originally thought that all plate glass was very flat and

that its flatness must be controlled to very close tolerances during

its manufacture. Conversations with production engineers at several

glass factories revealed, however, that most emphasis is placed on the

uniformity of thickness of the glass and not on its flatness. Essen-

tially, all standard plate glass is now produced by the flotation

process. This gives it a very smooth surface and uniform thickness and

it is simply assumed to be flat enough for most purposes. Numerous

pieces of commercially available plate glass were optically tested for

flatness, and it was found that some pieces are much flatter than

others. The tests were run by placing the carefully cleaned surfaces

of two glass plates together and illuminating them with monochromatic

light from a sodium vapor lamp. The uniformity of the interference

fringes produced gives a qualitative indication of the flatness of the

two glass panes. Three pieces of special plate glass were obtained

which were manufactured by the traditional process of mechanical grind-

ing and polishing. The manufacturer claimed that this glass was flat

to within two interference fringes per inch, where one interference

fringe represents a deviation from the mean surface of approximately

2.75 x 10-4 millimeters. When these pieces of plate glass were observed

under the sodium vapor lamp, their superior flatness was evident.

Selected pieces of the commercially available glass, however, were

nearly as good. These were used in the construction of the Hele-Shaw

model.

In order to provide for the controlled infusion of the viscous

liquids into the gap between the plates, it was necessary to glue a

section of machined plexiglass to each end of both glass plates. In








this way, an infusion reservoir was created on either end of the model

so that the liquids could be injected uniformly into the test section.

This is shown in Figure 4.2.

The gap between the two glass plates was maintained at a uniform

thickness by a series of brass spacers 1/16 of an inch (1.59 mm) thick

placed at intervals around the edges of the model. The plates were

held in position by C-clamps which pressed them tightly against the

spacers. The model was sealed by a length of surgical rubber tubing

which was placed around the edges, just inside of the spacers, and

crushed flat by the action of the C-clamps. The wall thickness of this

tubing was approximately 1/32 of an inch so that it could be crushed

between the plates without exerting much force against the glass.

The assembled model, mounted on its stand is shown in Figure 4.3.

The length of the test section is 76 centimeters and its height is 38

centimeters. With the rubber seal in place, however, the actual height

available for flow is 35.6 centimeters.

The viscous liquids are supplied to the infusion reservoirs on

either end of the model through plastic tubes. The plexiglass end

pieces on the back plate are each supplied with four fittings for 1/2

inch (1.27 cm) plastic tubing. The heavier fluid is conducted to the

model through this tubing from constant head tanks mounted on a ring

stand. This is shown in Figure 4.4. On the front plate, the end

pieces are provided with air vents on one end, and four nipples for

the connection of 1/4 inch (0.64 cm) plastic tubing on the other.

The lighter viscous fluid is injected through this tubing into one of

the infusion reservoirs by a positive displacement syringe pump. This

setup is illustrated in Figure 4.5.







40














V)u



V)
v-J

uL







-j Lu
V)A
w -









C C,
Q- ca
co









LU






VA
uLu
CC
w LLU


crc











'-4
(A
a:D













.-f








(A
(A
F-


VD




Cj


V)u
LL c LJ

























< ::D
__j



































V) CL
< I C-)r c






U-LL)Iv





























-)
0 <
>


LLII





































I~


I.







V







0



C
I-













CC
oc






LL
t)
I-
I-I














UJ
C
-I






f-i
LL







42









































LU
1' -4













in


ov U-

Lii

V)



Io
mr



LCC

U-












j A w









Properties of the Viscous Fluids

Dow Corning Corporation Series 200 silicone oil was used as the

heavier fluid in the model. This is a clear dimethyl siloxane which

is characterized by oxidation resistance and low vapor pressure. Its

specific gravity is 0.977 and its kinematic viscosity is nominally 500

centistokes at 250C. For the lighter fluid in the model, a mixture of

90 SAE and 30 SAE non-detergent motor oils was used. This mixture was

adjusted until its kinematic viscosity was the same as that of the

silicone oil at a temperature of 24.50C. The temperature versus

viscosity relationships for the two fluids in the neighborhood of 250C

is shown in Figure 4.6. This data was obtained with a Cannon-Zietfuchs

viscometer with calibration traceable to the National Bureau of Stan-

dards. The specific gravity of the light oil mixture was 0.890 at

250C.

These two fluids are not readily miscible, yet the surface tension

effects between the two are not so great that capillary action at the

interface is noticeable.


Operation of the Model

The Hele-Shaw model is a versatile analytical tool and its mode

of operation depends somewhat on the type of test being run. In the

present case the object of the tests was to observe the shape and rate

of rotation of the initially vertical interface between the two fluids

of different density. The idea was that the behavior of such an

interface would be similar to that of the interface formed by the

injection of freshwater from a well into a confined saline aquifer,

neglecting dispersion. Of course, the two processes are not the same

because the flow away from a well is radial while the Hele-Shaw model




























MOTOR OIL MIXTURE


DOW 200 FLUID


2 23 24 25 2
TEMPERATURE (DEGREES C.)


FIGURE 4.6.


TEMPERATURE VS. VISCOSITY FOR THE
TWO VISCOUS FLUIDS


540


530


520





510





500








is limited to two-dimensional flow in cartesian coordinates. Direct

application of the Hele-Shaw results to a prototype well situation is,

therefore, not possible. With increasing distance from the well,

however, the two flows become more and more similar so that the model

tests are useful for the qualitative description of interface movement.

The results of the model tests are also useful for comparison with

numerical models which can be made to work in either radial or

cartesian coordinate systems.

In setting up the model for a test, it is necessary first to be

sure that all parts of it are clean and dry. This can only be done by

completely disassembling it and washing each part. A commercial pre-

paration of trichloroethane cleaning solution was found useful in

washing the oils from the glass and plastic parts. After the model had

been cleaned and reassembled, with the plates rotated to the horizontal

position, the heavier oil was introduced at the downstream end; the

upstream end being the one into which the lighter oil would be injected.

During this process the upstream end was raised a few centimeters so

that all of the air in the test section would be expelled. Air

entrapped in the downstream infusion reservoir was released by loosen-

ing the brass plugs provided for this purpose. The introduction of

silicone oil was continued until the oil reached the end of the glass

and was just ready to run into the upstream infusion reservoir. At

this point, the flow was stopped and the oil front was held at the

edge of the glass by proper adjustment of the constant head supply

tank. The 1/2 inch (1.27 cm) plastic tubes on the bottom side of the

upstream end were then clamped off tightly with hose clamps. Three

of the 1/4 inch (0.635 cm) tubes from the syringe pump were then









hooked up and the upstream infusion reservoir was filled with motor

oil. The fourth nipple was left open for the escape of air. After

all of the air had escaped and the model was entirely filled with the

two viscous fluids, the fourth tube was attached and the interface

between the two oils was driven into the model to the starting point

for the test. It was found helpful to let the interface sit in this

position for several hours so that any irregularities in it would be

straightened out by gravity. The elevated end of the model was then

lowered to a level position and the test begun.

The test was started by simultaneously rotating the model to the

vertical position, turning on the syringe pump and starting the timer.

The movement of the interface could then be photographed until the

downstream end of the model was reached.

Several variations of this procedure were used. The first test

was run with the model left in its horizontal position so that the

straightness of the advancing interface could be checked. It was found

that, with higher rates of injection of the lighter fluid, the interface

moved faster in the center of the model. This indicated that the

higher pressure required to drive the viscous fluids across the model

could cause the glass plates to bend, changing the thickness of the

gap. Figures 4.7 and 4.8 show the shape of the interface when

advancing at rates of 3.65 cm/min and 7.28 cm/min, respectively. Since

the straightness of the interface at the lower speed was considered

acceptable most of the remaining tests were run at this rate of

injection.

Figures 4.9 and 4.10 show a test which was run at a higher injec-

tion rate with the model in the vertical position. The bowing out of

the glass plates caused the center section of the interface to advance






48




























I-t








hi.
0







i .i












(J
N





i t





L..




o.







49
































~~CD


uj




~' tLw
CD








I C)







F-





LLJ
iy


::Do

LL






















S
C-

0
f ~C' ---- I('4


Co
II










H




C-,
LU
r)





C-,


L~,LU

tc LA
rr~j !gta ~ --1 "~1~ 1""a::
lrLii
H
r, IILL





F a-


C-


-J
'--
~~cc
aI
I: U)
S,
'--

I i- i ~i~b- LsParC ~ I ~r~02~ '




k' ; LL






51
















I-I -'








I-
co

CD

A- 17)
7-









mmm.L--' ~i
7p,)

L-)


L-4







Lit!i






C

Li-
I t
-r- .- 'wE f









more rapidly than the ends. Rotation of the upper part of the interface

past the vertical results in instability, which is clearly evident here.

The sequence of Figures 4.11 through 4.16 shows a test that was

run with an injection rate of 20.44 cubic centimeters per minute.

This is the same injection rate that resulted in the relatively straight

interface shown in Figure 4.7 when the model was horizontal.

The sequence of Figures 4.17 through 4.23 shows a test in which

only the lower half of the model was used. This was accomplished by

inserting a barrier composed of rubber splicing compound and silicon

rubber sealant between the glass plates to reduce the width of the test

section by half. The injection rate used here was 8.16 cubic centi-

meters per minute.


Analysis of Test Results

The rotation of an initially vertical interface between two fluids

of different density has been studied in several previous investigations.

Gardner, Downie and Kendall (1962) proposed an equation to predict the

rate of interface rotation which, for an isotropic aquifer and fluids of

equal viscosity, has the form


()2 16 (t/t')2 (4.31)
3 1 + t/t'


where


X = projection of the interface on the horizontal

b = thickness of the aquifer

t = time

t' = 4/3 nbp
KAp
p = average density of the two fluids.













































-*-1-~t- !
- I -~- ~
4 -
I '



1


~~1


jI,

ji


4-. 7
.i13
r4


II



C

















0
I-'
LL_
I:













I-
U-


I-
I-


LU

NE


(I i ~'

i '."







I~ sl





























* 1 2





It
* 1- i .


hi 1 1

it


I,


I
* 1~


* 4


, i


I' "


*.~ 7.
' 1'L~lr


I


4 4


z



C)


d-
cx

0
F-


LU











-<
a:
I-
0


z:



I-



cC
LL



I-i
'-


IItr jr


t


4























tz
S-t---, .. I'' -, -t ,- C --k ...: .
"t~fi -e I 4 j / ;l', ir
! i [ i '






I ~ ;~ rP-4
-t .- --t .. '- .-LU
r 4 j ,r ,

F- 1 1 V ,
4. I I .C
4..., II -


*1 *. -J
I ; l I :, I '









i-iin
---i- --V -- ...... L --- -
4 If.. i '- ii



. : .:

1. 1; I
4- 1 ". 1 .-.-



,. I ..- 1 .





,,,"4 m*-
I I
a:J

'~ LI Cr
I I
j74.i ~~dls 1ik 4iL a:
IYC'-LFlr LU
H
r t Iv -
I~~* >40c

)fi
~t< C'
1.2 3
1*- *- IL






























-J .
I I




'1,


t t


1.1.


-"I

a'.i


:.1


I-


#sJS i
b0y
jL~

i i,,


I'


^ ,





S



Cj
*1 C)'





ItI
-1-~ C-i -









I-
-t I :Fi 44U41 t -rS ,



C-P

















LL
=D
. .
,t -





p







58


















4'



c-1 L U1

iN- --Vc


S- I 41- -



Cf)





LUJ


V)


LUJ
jo-





LUJ

IL-








C-0
F;~j i-- -i-LL.







59




























00


C-,
C)















I-
(L, <
m


14










IV z







LU
7o







LL
!C)
CL













LLUj
i' I-







4.-




LUJ
D
il- -i

















Cl,
UL-

.L LJ






60






















co

,T A rII


O~LL-
~ Cf. P I-cI
F- <


LLu



Li u
Lu J
CK

U- ui

u L
2--
E


i~ Lu







Lu

Li'..







61






















00


--1 II




ILU





2t7T9.IILL.
V; u
Lii


` `-i. -
~-lLLJ







Q-I











(!J
Li)
i' -r-- Fi' -cc
I I'


I c W

i I !I i ~ Uo

1: i~i L-







62




















I-I




W.




co





II
rc



C)2
LU -


'-4

C)


LU Z










LU
L:









IL
rll cOtf


F 4r






Ar,,









63



























C6)




C-)




LU
co




Li..

-f-i


*D LUJ
rl





LL- LUI








Li..

LUJ









LUI



LL-







64
















z






06






*LUJ

CA)
LU
-liJ-J
'-t-i~l~~1'-4
t LLJ







-'-4
oJ







C'3



LU
t =D L J













a:

LLJ
L J
.j r L Li o
LL
r,
s i r'll~ =









1.4







65































-i;--L.) c tI
'" "! liLl


L~LU




tii
4,i F t- -J
LL- uj
CD
g ii-~: z >
r .;7:~ cLl
Z;S 4'

~ r~tY~dig1138 ~ s i

LL
tl~~ If~g~~ i I vCY
ZD
't) ICD
LL.








Equation (4.31) was constructed intuitively by joining the solutions to

the extreme cases of nearly vertical and nearly horizontal interfaces.

The solutions to these extreme cases were based on the assumptions of

no vertical flow and no curvature of the interface.

Another equation for the movement of an initially vertical inter-

face was given by Harleman and Rumer (1962). Once again, assuming a

straight interface and no vertical flow, they obtained


X 2t
S = 8/3 (4.32)


This equation is of the same form as the one for nearly horizontal

interfaces from which equation (4.31) was derived.

Figure (4.24) shows the comparison of equations (4.31) and (4.32)

with data taken from three runs of the vertical Hele-Shaw model. Runs

number two and three are shown in Figures 4.11 through 4.16 and 4.17

through 4.23, respectively. Run number one was performed in the same

manner as run number two except that the injection rate was twice as

high. The third run gave the best fit to equation (4.31). Runs number

one and two, while fairly consistent between themselves, do not fit

either equation very well. Since the main difference between these two

runs and the third run is in the depth of the flow region, it is most

logical to suspect that it is caused by the existence of vertical

velocity components which were ignored in the derivation of the

equations. The greater curvature of the interface which is apparent in

the first two runs tends to support this view.






67







3-


I I
2





S.8
I "-Y -- T 1 -- i -t .- ....- .. / "

-' .. ..


. HARLEMAN & RUMER


.4


GARDNER, DOWNIE & KENDALL

.2


S----RUN 1
o ----RUN 2
.I ----RUN 3
.04 .06 .08 .1 .2 .4 .6 .8 1

t/t'

FIGURE 4.24. COMPARISON OF EXPERIMENTAL RESULTS WITH
EQUATIONS FOR INTERFACE ROTATION













CHAPTER 5

SIMPLE CASES INVOLVING STEADY FLOW


General
It was shown in Chapter 3 that the general problem of locating a

moving interface between two fluids of different density is very complex

even when the standard simplifying assumptions involving homogeneity,

isotropy and immiscibility are used. It requires the solution of two

coupled flow problems, one in the freshwater and one in the saltwater,

joined by a nonlinear boundary condition on the interface, of which the

position is unknown. The only methods presently available for getting

a solution to such a problem are by physical and numerical modeling.

However, if the problem is limited to the special case of a stationary

freshwater lens in a static saline aquifer, then it becomes much more

tractable. In this special case, it is only necessary to solve the

continuity and motion equations in the freshwater region and the location

of the interface, while still not initially known, can be treated by the

Ghyben-Herzberg principle. This makes it possible to formulate some

problems which can be solved analytically.

Steady flow in a freshwater lens is not possible unless the aquifer

possesses a constant head boundary which permits freshwater to flow

out of the lens at the same rate that it is being injected. One such

constant head boundary is found at the sea coast where injection wells

may be used to create freshwater lenses in aquifers affected by saltwater

intrusion. Another case is that of injection into a leaky saline aquifer









where the leaky layer serves as a boundary which permits efflux of

freshwater from the lens. Even in these cases, steady flow may only

occur during a small part of the history of the lens. It represents a

condition of equilibrium in which buoyant forces are counteracted by the

loss of energy due to the flow resistance of the porous medium. It is

useful to study this condition because it indicates the limiting size

of lens that can be achieved under the given flow conditions.

It should also be noted that, in the operation of a storage lens,

certain advantages can be gained by maintaining the lens in a steady

state through continuous freshwater injection. The first is that

continuous injection is necessary to preserve the potential energy of

the lens so that the maximum amount of freshwater can be withdrawn when

it is needed. Note that as soon as injection is stopped, the lens

begins to decay into a thin layer of freshwater at the top of the aqui-

fer which could not be selectively retrieved even if fresh and saltwater

did not mix. The second advantage is that the continuous injection of

freshwater to the lens tends to wash out the saltwater which had re-

mained in the pores and adhered to the particles of the porous medium

after the passage of the mixing zone as the lens was being formed. In

this way, the salinity of the freshwater and the thickness of the

mixing zone are reduced.


Freshwater Injection in a Coastal Aquifer

Saltwater Intrusion

It is commonly observed in coastal aquifers which discharge to the

sea that a wedge of nearly static saltwater extends inland from the

coast underneath the flowing freshwater. Actually, it was the study of








this phenomenon which led W. Badon Ghyben and A. Herzberg to develop

their hydrostatic theory for interfaces. Numerous subsequent investi-

gators have brought forth more refined methods for analyzing saltwater

wedges (Glover, 1959; Henry, 1959; Bear and Dagan, 1964; Van Der Veer,

1977) and have found the Ghyben-Herzberg results to be quite accurate

except for that part of the interface closest to the sea.

The application of this principle, combined with the Dupuit assump-

tion, to a confined aquifer of constant thickness gives the following

relationship between h and x.


x = Kh2/(26Qx) (5.1)

in which the symbols are as shown in Figure 5.1. The position of the

toe of the interface is given by

-Kb2
x = -2 (5.2)
t 26Q

Injection Into a Coastal Aquifer

If an injection well pumps freshwater at a constant rate into the

coastal aquifer pictured in Figure 5.2 then the equilibrium position of

the saline wedge will be modified and a steady freshwater lens will be

formed in what was formerly a saline part of the aquifer. This is a

case of injection into an aquifer in which only the freshwater is flow-

ing, so that the Ghyben-Herzberg principles can still be used. The flow

pattern around a well in a uniformly flowing aquifer is customarily

studied by the methods of potential theory in which the velocity poten-

tials for source flow and uniform flow are superimposed. However, in

the presence of the saltwater wedge, even though the vertical velocity

components are neglected, flow is not two dimensional because of the







71





























o
N I


-


L. re"

I LA-
'--





\LUl
i tn /--
u I \- I

Cz
I

u ,






O I I-




L I
< Lu
Lu I



SICD



-. i
(1 I -y





I-- X
- I I



C I II ----
N



I ~U
; I L .
u I e8







72













Sll*







LL




--




LX /



-I--
LiiL
N 0j


LL




0 i
^ /,- ^ -- '"- --^ -----L
.I.I I L.C)
Lu "J I \





L I \
--------- U
u L\






/ I
- |
O ^ '









changing thickness of the flow region. Strack (1976) has found a

modified form of the vertically integrated potentials originally pro-

posed by Girinskii (see Bear, 1972; p. 157) which can account for the

presence of the interface. If the thickness of the flow region is a

linear function of the piezometric head, 4, i.e.,


h = #i + 8 (5.3)

then Strack gives the potential for the interface zone as


Di = 1/2aK(4 + B/a)2 (5.4)

For the inland part of the aquifer where there is no interface, the

potential is defined by


P = Kb) + c (5.5)

where the constant, c, is used to make the potential function continuous

at the toe of the interface. It is easily verified from equations (5.4)

and (5.5) that the discharge components in the x and y directions are

obtained from 0 in the usual way.


Qx = -a3/ax and Qy = -D/ay (5.6)

in which Q and Q represent the discharge through a unit width and the
x y
total thickness of the aquifer in the x and y directions, respectively.

From equations (5.6) it follows that the continuity condition for the

confined coastal aquifer is given by the Laplace equation


a20/ax2 + a2o/ay2 = 0 (5.7)

and 0 is a harmonic function of x and y. The effect of this is that








through use of Girinskii's potential, variations in the vertical direc-

tion have been integrated out of the problem. The resulting potential

function can be manipulated, using the theory of complex variables to

describe a two-dimensional flow field with the desired boundary condi-

tions. At the same time, the numerical value of the potential function

at any point gives the depth to the interface at that point through

the Ghyben-Herzberg relation and thereby supplies three-dimensional

information about the freshwater lens. Before this can be done, however,

it is necessary to determine the constants a and B in equation (5.4) and

c in equation (5.5).


Confined coastal aquifer. From the geometry of the saltwater wedge,

as shown in Figure 5.2, it is seen that

s = b a (5.8)

The thickness of the freshwater region, h, is related to the "push up",

s, through

h = 6s a (5.9)

which is a consequence of the Ghyben-Herzberg principle.

The combination of equations (5.8) and (5.9) results in


h = 64 6b (1 + 6)a (5.10)

But it has already been assumed that the relationship between h and <

in the interface zone is given by equation (5.3). Consequently the

values of a and B must be:


a = 6 and 0 = 6b (1 + 6)a (5.11)

Substituting these values into equation (5.4) gives the definition of








the potential function in the interface zone for a confined coastal

aquifer.

i 6K[p b (1 + 1/6)a]2 (5.12)


Now, in order that the potential function may be continuous throughout

the aquifer, it is necessary to determine c in equation (5.5). Note

that at the toe of the saltwater edge


s = (a + b)/6 (5.13)

which along with equation (5.8) gives


: = (a + b) (1 + 1/6) (5.14)

Using this value of p in both equations (5.5) and (5.12) and requiring

that they be equal at the toe of the wedge results in


c = 2- Kb(a + b) (1 + 1/6) (5.15)

The definition of the potential function in the inland zone of a

confined coastal aquifer then becomes

P = Kb[b + 1/2 b/6- (a + b) (1 + 1/6)] (5.16)

At the toe of the saltwater wedge this reduces to


(t = 1/2 Kb2/6 (5.17)

In order to construct the complex potential function for flow in

the aquifer the boundary condition at the sea coast must be known.

From Figure 5.2 it is seen that s = a/6 at the coast. Substituting this

into equation (5.12) gives the boundary condition








( = 0 at x = 0 (5.18)

This condition can be satisfied by locating an image well, with the

same strength as the injection well but opposite sign, on the other side

of the boundary at the same distance from it, x,, as the injection well.

The superposition of the velocity potential for the two wells and the

uniform seaward flow results in

(x Xw) y22
= -xQ + In ( )2 y2 (5.19)
x 47r (x + x) 2+ y2


When particular values of 0 are computed from equation (5.12) and subs-

tituted into equation (5.19) which is then solved for the particular

equipotential line, the shape of the freshwater lens produced by the

injection well is outlined. Such a case is shown in Figure 5.3.

From a practical point of view it is most important to locate the

toe of the saltwater wedge around the periphery of the lens because this

is the outer boundary of the usable portion of the freshwater. The

equipotential line which delineates this outer boundary has the value

given by equation (5.17). Solutions for this line can be made dimen-

sionless by defining the following dimensionless parameters:


x' x (5.20a)
t

y, = Qx y (5.20b)
It
x Q x (5.20c)
x w t w

Q. = Q (5.20d)
(t

A series of such dimensionless solutions for various conditions of

injection is shown in Figure 5.4. and Figure 5.5.








































-.)


z cr


S V


SI V



'I














0

C\J









r.-

II




Lfl






SII
i"
C)



LU
(-5








I-
LLI





Z










I-m
LL








Lu
C)

Lu
-j
cr:
VL



I-4





2-:
LJ
LI-I

LO

C,







L-
u-


o Cn


x($/Ix-)= ,x-






























Lfl

















LO

II








U)


X(,/Xb-)= ,x-








Phreatic coastal aquifer. The Girinskii potentials defined by

Strack for interface flow can also be used in aquifers which have no

upper confining layer because the thickness of the freshwater region is

still a linear function of the piezometric head. The linear function

is not the same as for a confined aquifer, however, so the potentials

must be redefined to conform to the aquifer geometry shown in Figure

5.6. In the inland zone the specific discharge, Qx is given by

K2
Qx -K K2 (5.21)
ax

The potential for this region, then, should be defined as


D = 1/2 KD2 + c (5.22)

In the interface zone, note that


s = ( b (5.23)

and from the Ghyben-Herzberg principle the thickness of the flow region

is

h = s(l + 6) (5.24)

which leads to


h = (1 + 6)( (1 + 6)b (5.25)

If this is to be in the form of equation (5.3) then the constants a and

6 must be


a = (1 + 6) and B = -(1 + 6)b (5.26)

Substituting these values into equation (5.4) results in a definition






81














UJ



rw

LU
C,
0



--









LLJ < /--c
SI LJ l
LU I_

LU <







r _j I u.
SLLU

I I I-
ULU




LU 0 II-





F-



------4- -------





















r \
0 D<



LU I------~
-4




0Z








for the velocity potential in the interface zone


Di = 1/2 (1 + 6) K(q b)2 (5.27)

At the toe of the saltwater wedge the piezometric head is


t = b(l + 1/6) (5.28)

Substituting this into equations (5.22) and (5.27) and requiring that
be single valued at this point leads to the value of the constant c in
equation (5.22).


c = -1/2 Kb2 (1 + 1/6) (5.29)

The definition of the velocity potential in the inland zone of the
phreatic aquifer .then becomes


S= 1/2 K[P2 (1 + 1/6)b2] (5.30)

At the toe of the saltwater wedge this reduces to



Dt = 1/2 Kb2(l + 6)/62 (5.31)

which is greater than the corresponding value for a confined aquifer by

a factor of (1 + 1/6).

Axisymmetric Flow in a Leaky Aquifer
Consider a homogeneous saline aquifer with hydraulic conductivity,
K, of uniform thickness, b, and unlimited area extent. It is confined








at the bottom by a horizontal impervious layer and at the top by a

horizontal aquitard of thickness, b', and hydraulic conductivity, K'.

A fully penetrating well injects freshwater at a constant rate, Q,

forming a freshwater lens, which remains centered on the well since

there is no pre-existing flow in the surrounding salt water. Injection

continues until a steady state is eventually reached in which the injec-

tion rate is matched by the rate of leakage through the aquitard. The

mixing zone at the boundary of the lens then becomes relatively thin

and can be approximated by a sharp interface. This situation is

depicted in Figure 5.7.

The flow region is divided into two zones, each of which has its

own governing equation. In zone 1, which is occupied entirely by

freshwater, the customary analysis for well flow in leaky aquifers can

be used. This includes the assumption of 90 degree streamline refrac-

tion at the aquitard so that flow in the aquifer is horizontal while in

the aquitard it is vertical. The governing equation is (see De Weist,

1965, p. 264)

d2s 1 ds 1
d2 + d- )2 s = 0 (5.32)
dr2 r dr
where

B2 = Kbb'/K'

This is a modified Bessel equation of order zero, for which the

general solution is given by


s = C1I0(r/B) + C2KO(r/B) (5.33)

where 10 and K0 are the modified zeroth order Bessel functions of the

first and second kinds, respectively. The boundary condition at the