A study and evaluation of saltwater intrusion in the Floridan aquifer

Material Information

A study and evaluation of saltwater intrusion in the Floridan aquifer by means of a Hele-Shaw model
Evans, Andrew Joseph, 1940-
Publication Date:
Physical Description:
xiv, 213 leaves : ill., diagrs., maps ; 28cm.


Subjects / Keywords:
Aquifers ( jstor )
Cradles ( jstor )
Fresh water ( jstor )
Groundwater ( jstor )
Limestones ( jstor )
Modeling ( jstor )
Parametric models ( jstor )
Pumps ( jstor )
Saltwater ( jstor )
Velocity ( jstor )
Civil Engineering thesis Ph. D
Dissertations, Academic -- Civil Engineering -- UF
Saltwater encroachment -- Florida ( lcsh )
bibliography ( marcgt )
non-fiction ( marcgt )


Thesis--University of Florida.
Bibliography: leaves 210-212.
General Note:
General Note:
Statement of Responsibility:
by Andrew Joseph Evans, Jr.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Resource Identifier:
000158804 ( ALEPH )
02609315 ( OCLC )
AAS5126 ( NOTIS )

Full Text








To my supervisory committee, I wish to extend my

sincere appreciation for their efforts on my behalf. A

special thanks goes to Dr. B. A. Christensen and Dr. J. H.

Schaub for their continuing support and advice. It's not

an easy task to pull an oar in Hagar's longboat.

The efforts of Mr. C. L. White, Mr. William

Whitehead, and the staff of the Mechanical Engineering

machine shop in the construction of the model are grate-

fully acknowledged, as well as the assistance of Mr. Richard

Sweet and Mr. Tom Costello'in the operation of the model, and

the mathematical advice of Dr. Jonathan Lee. The participa-

tion of Mr. Floyd L. Combs in all aspects of the project was

especially valuable.

Without the financial support of the Office of Water

Resources Research, United States Department of the Interior,

and the administrative assistance of Dr. W. H. Morgan and

the staff of the Florida Water Resources Research Center, this

project would not have been launched.

Finally, to my mother and my many dear friends who

have watched the "ins and outs" of my academic and profes-

sional career, bless you for your interest and encouragement.

I want to leave the reader with this closing thought,

which I believe is paraphrased from Benjamin Franklin:

"Waste not want not, you never miss the water 'till the

well runs dry."




ACKNOWLEDGEMENTS ...................................... ii

LIST OF TABLES ........................................ vi

LIST OF FIGURES .... ........................... ........ vii

KEY TO SYMBOLS OR ABBREVIATIONS ........................ x

ABSTRACT .............................................. xiii


Topography ................................. 2
Western Highlands ....................... 3
Marianna Lowlands ........................ 6
Tallahassee Hills ....................... 7
Central Highlands ....................... 7
Coastal Lowlands ........................ 7
Climate .................................... 8
Geology .................................... 12

A: NUMERICAL METHODS ............................. 15
Method of Finite Differences ................ 15
Method of Finite Elements ................... 17
Relaxation Methods ......................... 18
B: PHYSICAL MODELS ............................... 19
Sandbox Model .............................. 21
Hele-Shaw Analog ............... ............. 22
Electric Analog .............................. 23
Continuous Electric Analog .............. 23
Discrete Electric Analog ................ 24
Ion Motion Analog ....................... 24
Membrane Analog .............................. 25
Summary .................................... 26


Chapter Page

III. THE HELE-SHAW MODEL ........................... 28
Viscous Flow Analog ........................ 29
Scaling .................................... 37
Time .................................... 39
Anisotropy .............................. 40
Leakage ................................. 50
Storativity ............................. 54
Discharge ............................... 55
Accretion ............................... 59
Volume .................................. 59

AND HYDROLOGY .............................. 61
Site Selection ............................. 61
Prototype Geology .......................... 67
Prototype Hydrology ......................... 69

OF MODEL ................................... 71
Design ..................................... 71
Prototype ............................... 71
Model ................................... 71
Construction ............................... 86
Frame ................................... 86
Plexiglas Plates and Manifolds .......... 105
Saltwater System ........................ 125
Freshwater System, General .............. 125
Freshwater System, Accretion ............ 132
Freshwater System, Wells ................ 132
Freshwater System, Flow Meters .......... 132
Operation ........ .......................... 156

Results .............. ....................... 158
Conclusions and Recommendations ............ 191


RECTANGULAR CROSS SECTION .................. 193

B. FLOW METER CALIBRATION ........................ 195

BIBLIOGRAPHY .......................................... 210

BIOGRAPHICAL SKETCH ................................... 213


Table Page


5.1 PROTOTYPE PARAMETERS ......................... 72

5.2 MODEL PARAMETERS ............................. 81

5.3 SIMILARITY RATIOS ............................ 82

6.1 DEPTH TO SALTWATER ........................... 186


Figure Page

1.1 Topographic divisions of Florida ................ 5

1.2 Mean annual precipitation ...................... 11

3.1 Free body flow diagram for Hele-Shaw model ..... 32

3.2 Section of anisotropic grooved zone in
Hele-Shaw model ............................. 44

3.3 Head loss in grooved anisotropic zone ........... 46

3.4 Section of leaky zone in Hele-Shaw model ....... 53

3.5 Storage manifold for Hele-Shaw model ........... 57

4.1 Regional area of prototype ..................... 64

5.1 Ghyben-Herzberg interface model ................ 75

5.2 Cradle and cradle dolly ........................ 88

5.3 Cradle rotation ................................ 91

5.4 Frame and model setup .......................... 93

5.5 Air hose and valve arrangement ................. 96

5.6 Stub shaft and pillow block arrangement ........ 98

5.7 Back up air supply .............................. 100

5.8 Internal support and sealing system ............ 102

5.9 Model mounting system .......................... 104

5.10 Model back and front plates .................... 107

5.11 Detail of the model front and back plate ....... 109

5.12 Detail of the model back plate ................. 111



Figure Page

5.13 Front plate with accretion manifolds ........... 114

5.14 Accretion manifolds ............................. 116

5.15 Detail of accretion manifolds .................. 118

5.16 Connections between model and fluid
supply system ............................... 120

5.17 Saltwater constant head tank ................... 122

5.18 Back and front plate clamp up .................. 124

5.19 Fluid supply network schematic ................. 127

5.20 Saltwater reservoir and pump ................... 129

5.21 Freshwater supply system ....................... 131

5.22 Freshwater reservoir and accretion pump ........ 134

5.23 Well supply manifold and pump .................. 136

5.24 Well supply manifold and accretion
supply manifold ............................. 138

5.25 Opposite view of Figure 5.24 ................... 141

5.26 Flow meter bank ................................ 143

5.27 Flow meter detail .............................. 145

5.28 Flow meter to model connections ................ 148

5.29 Flow meter pressure sensing lines .............. 150,

5.30 Flow meter switching device and
pressure transducers ........................ 152

5.31 Carrier demodulator and strip chart
recorder .................................... 155

6.1 Interface location, tm = 0 min. ................. 160

6.2 Interface location, tm = 32 min. ............... 162

6.3 Interface location, tm = 48 min. ................ 164
























tm = 62 min.

t = 79 min.

t = 92 min.
tm = 102 min

tm = 107 min

t = 117 min

t = 132 min
tm = 147 mir

tm = 162 mir

tm = 172 mir





6.14 Wedge location at time zero .....


............... 166

............... 168

............... 170

.............. 172

........ ..... 174

..... ......... 176

.............. 178

.............. 180

.............. 182

.............. 184

............... 188

6.15 Time overlay of wedge location ..














A = area, Fourier coefficient

b = width of model and prototype

B = body force, Fourier coefficient

f = subscript denoting freshwater, Floridan aquifer

g = acceleration due to gravity

h = height of water

j = summation limit

k = intrinsic permeability

.K = hydraulic conductivity

1 = distance between storativity tubes, subscript
denoting leaky layer

L = distance center to center of groove in anisotropic
zone, flow meter tube length

m = subscript denoting model

n = effective porosity

p = pressure, subscript denoting prototype

q = specific discharge

Q = total flow

r = subscript denoting ratio

R = accretion

R = Reynolds number



S = specific storage

t = time

T = transmissivity

U = volume

V = velocity

x = horizontal direction parallel to test section

y = horizontal direction perpendicular to test section

z = vertical direction

1,2 = subscripts denoting zone 1 and zone 2

a = width adjustment factor

y = unit weight

A = length adjustment factor

p = absolute viscosity

E = geometric parameter in anisotropic zone

p = mass density

v = kinematic viscosity

= potentiometric head

= velocity potential, potential


D. substantial derivative
V2 = Laplace operator

A = difference operator

C = Centigrade



cfs = cubic feet per second

*F = Fahrenheit

fps = feet per second

g/cm3 = grams per cubic centimeter

gpd = gallons per day

Hg = Mercury

I.D. = inner diameter

MGD = million gallons per day

msl = mean sea level

O.D. = outer diameter

pcf = pounds per cubic feet

psid = pounds per square inch differential

psig = pounds per square inch gage

RPM = revolutions per minute


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



Andrew Joseph Evans, Jr.

August, 1975

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

Continuing development of the coastline zone in the

middle Gulf area of Florida is increasing the demand for

groundwater supplies and, in turn, increasing the probability

of saltwater intrusion. Methods must be developed to make

long-range predictions on the effects of increased demands

on the Floridan aquifer.

A Hele-Shaw model is a physical model which fits the

requirements for long-range planning. It is well-suited for

handling anisotropic aquifers, difficult boundary conditions

and can simulate years of field conditions in minutes of

model time.


The section selected for study lies in a line from

the Gulf coast near Tarpon Springs to a point near Dade City

and passes through the Eldridge-Wilde well field. The

Eldridge-Wilde well field is the major water producer for

Pinellas County. This region has experienced several years

of dry weather, and pumping has lowered the water levels in

the aquifer by a significant amount. This loss of fresh-

water head is certain to induce saltwater intrusion.

A Hele-Shaw model has been built for this area, and

all pertinent geological and hydrological features of the

area are included. Steady state characteristics of the

aquifer system have been considered. In particular, the

long-term effects due to pumping and artificial recharge

were examined.



Florida, with the possible exception of California,

is the fastest growing state of the United States. The

rapid influx of people since World War II has greatly

increased the demands for land and water. In the past,

there has been an almost total lack of wide range planning

for the uses of these resources. Even fewer investigations

have been made into the consequences of their rapid and

unordered development. Recently, water supplies have had

to be rationed in South Florida. Overall, land and wet-

lands required for fish and wildlife have so diminished

that, in some instances, there has been a marked decrease

in their numbers. It seems reasonable to conclude that in

some areas of the state, land and water resources cannot

support much larger populations with current locally avail-

able supplies without almost irrecoverable damage to the

groundwater system in the form of saltwater encroachment.

With the growing affluence of the American people,

and the availability of economically priced air conditioning

units, it can be expected that even more people will leave

the colder northern climates for the southern and western

states. Florida can expect to receive more than its share

of the migration. Frequently now, environmental protection

groups are making forecasts of impending doom. At worst,

their predictions may come true and people are beginning to

look at all growth with a jaundiced eye.

It is doubtful, however, that growth and development

can be stopped. The history of man indicates a continual

effort to better his life-style, his private environment.

There is little doubt that this has sometimes caused a

degradation of other portions of his world. Unless the

cessation of all growth and development is acceptable, new

ways must be found of forecasting, or predicting, the results

of all growth so as to combat possible undesirable results.

Consequences of all growth must be known, even of those

resulting when pragmatic short-term solutions are used.

Hopefully, the remainder of this study will present

a modeling method which will be useful in forecasting the

results of pumpage and use of groundwater in our coastal

zones so that we may better plan their usage. But first,

a little background on Florida.


Florida lies between latitudes 24-40' and 31-00'

north, and longitudes 80-00' and 87*-40' west, and is the

most southerly unit of the continental United States. In

its southernmost extension it is less than 10 of latitude

north of the Tropic of Cancer.

Florida is bounded on the east by the Atlantic Ocean;

on the south by the Straits of Florida and the Gulf of Mexico;

on the west by the Gulf of Mexico and the state of Alabama;

and on the north by Alabama and Georgia. The shape of the

state in relation to the remainder of the United States

suggests two distinctive parts: the Floridan panhandle and

the peninsula of Florida. The panhandle is a strip roughly

225 miles long that stretches in an east-west direction.

The peninsula is a south-southeast extension at approximate

bearing S 17 east. From the northern boundary of the state

to the tip, not including the chain of keys, the peninsula

is approximately 415 miles long and includes two-thirds of

the land mass of the entire state. Its coastline, some

1350+ miles long, is the longest with the exception of

Alaska. No place in the interior of Florida is more than

60 miles from either the Gulf or the Atlantic coast.

Cooke (1939) divided the terrain of Florida into

five sections, Figure 1; the Western Highlands, the Marianna

Lowlands, the Tallahassee Hills, and a narrow band of Coastal

Lowlands, which comprise the panhandle, and the Central High-

lands and Coastal Lowlands, which comprise the peninsula.

The topography of each is described briefly.

Western Highlands

Extending eastward from the Perdido River (the

western boundary of Florida) to the Apalachicola River, the

northern part of this region near the Alabama state line


Topographic divisions of Florida.

Western Highlands

Marianna Lowlands
Tallahassee Hills
Central Highlands
Coastal Lowlands


0 50 100 150
Scale of Miles

After: Cooke (1939)

is not much higher than 300 feet. It is considered to be

hilly when compared to the broad gently rolling southern

parts of this region which drop to 100 feet elevation as

one approaches the coastal lowlands. The highest elevation

in the state, 345 feet, is found in this region in the north-

west corner of Walton County. The western highlands are

underlain with the sand of the Pliocene Citronelle Formation.

The steepness of the bankslopes at the headwaters of the

many streams is the most unique physiographic characteristic

of this section.

Marianna Lowlands

This roughly triangular-shaped region of Holmes,

Jackson, and Washington counties, with somewhat smaller

contributions from Bay and Calhoun counties, lies between

the Tallahassee Hills and the Western Highlands. It is

difficult to recognize this area of gently rolling hills

as lowlands. Cooke (1945) attributes the lower elevations

to the solubility and consequent degradation of the under-

lying limestone. This area is one of the two in the state

where the Ocala Formation is exposed to the surface and the

only area of the state where the Marianna limestone, the

soft white limestone of the Oligocene Age, is found exposed.

The region is dotted with sinks, sinkhole lakes, and springs.

Tallahassee Hills

From the Apalachicola River east to the Withlacoochee

River, the Tallahassee Hills extends along the Georgia-

Florida border and is only 25 miles in width. The western

section is a nearly level plateau some 300 feet above mean

sea level. The remainder consists of rolling hills carved

out of the Citronelle Formation. In addition to this, a

red clayey sand and Fuller's earth of the Hawthorn Formation

are found in this area. This is a fertile farming region.

Central Highlands

The Central Highlands forms the backbone of the

Floridan peninsula and extends from the Georgia line between

the Withlacoochee and St. Mary's rivers south-southeastward

some 250 miles into Glades County west of Lake Okeechobee.

This region is highly diversified. It includes high swampy

plains, hills, and innumerous lakes. Soils are sandy. Many

of them were derived from Pleistocene (Ice Age) marine ter-

races. However, a distinguishable amount comes from the

Miocene Hawthorn and Pliocene Citronelle Formations. The

lakes and sinks which dot the entire area indicate the

presence of limestone below the surface. Elevations of

this region average just slightly more than 150 feet; however,

they vary from less than 100 feet to approximately 300 feet.

Coastal Lowlands

The Coastal Lowlands, or Coastal Plains as it is

sometimes called, borders the entire 1350 mile Florida

coastline. Flanking on both sides of the Central Highlands,

the Coastal Lowlands is widest just south of Lake Placid and

narrowest between the western border and the Choctawhatchee

Bay just south of the Western Highlands. The elevations

everywhere within this region are less than 100 feet. The

soil for the most part is sandy except in the Everglades and

Big Cypress Swamp locales where Pliocene limestone, muck,

and peat prevail near the surface. The keys, which extend

some 100 miles into the Straits of Florida, are mostly sandy

oolitic limestone like that of the mainland; however, some

limestone with coral heads is found. The islands seldom

reach 15 feet elevation. The entire region is generally

flat, typical of recently deposited material with little or

no erosion.


The sea surface temperatures east and west of Florida

average, respectively, 780 and 77 Fahrenheit. Water tem-

peratures range from 740 to 83 Fahrenheit in the east and

700 to 84* Fahrenheit in the west. The coldest month in

both cases is February; the warmest month is likewise

August. The relatively homogeneous distribution of sea

temperature, the lack of high relief and the peninsula shape

of Florida contribute to its climate. The temperature is

everywhere subtropical. Mean annual average temperature in

the north is 68 Fahrenheit, and in the southern tip 750


The tradewinds, which shift from northern Florida to

southern Florida and back semiannually, bring a mildly mon-

soon effect to Florida. In November, the tradewinds are at

their southernmost extension and Florida's climate is con-

trolled by frontal, or cyclonic, activity moving in from the

continental United States. Rainfall during this period is

of low intensity and long duration.

Beginning in early May, the tradewinds move north,

again bringing with them the moist warm air of the Atlantic.

The cyclonic activity is greatly reduced over the state and

convectional instability begins to become established. June

through September is known as the rainy season in Florida.

The thunderstorms of this period are intense and very

specially varied. They usually occur during the hottest

part of the day and only on rare occasions last longer than

two hours. About 60 percent of the total average annual

rainfall occurs during this period, Figure 2. The mean

average rainfall of Florida is in the neighborhood of 53

inches. It varies from 38 to 40 inches in the lower keys

to over 65 inches in the southeast corner of the peninsula

and the western portion of the panhandle. Most of the

interior, that is the central highlands, receives approxi-

mately the mean annual average.


Mean annual precipitation (in inches).

I Over 66

S62 66

S58 62

I 54 58

50 54

S46 50

LZ Less than 46

0 50 100 150
Scale of Miles

^- <- -

Note: Precipitation normals
compiled from records
published in Climatological
Data: Florida Section,
December, 1971.

_I~_ _



The Floridan peninsula and the offshore submerged

lands above 50 fathoms, which Vaughan (1910) called the

Floridian Plateau, have existed for several million years.

The region has not been subject to violent earth movement

and, consequently, there has been a gentle doing resulting

in the formation of an oval arch above the basement rock.

The rock of the core underlying the plateau is hypothesized

to be pre-Cambrian; however, no drill has penetrated the

core. The oldest rocks penetrated, to date, are quartzites

found at about 4500 feet below the surface in Marion County.

The borehole encountered another 1680 feet of quartzite

before drilling was suspended. This metamorphized rock,

believed to be a continuation of the Piedmont region of

Georgia, was assigned by Cooke (1945) to the Pennsylvanian

period. The arch above the metamorphized basement, composed

of almost pure porous limestone, is known as the Lake City,

Avon Park, and Ocala Formations. Dated in the Eocene period,

the Ocala Formation has an estimated maximum thickness of

360 feet. It is found at, or above, mean sea level through-

out northeast and north central Florida and is this section's

principal aquifer. In southern Florida, in the vicinity of

the Everglades, the Ocala is found at depths approaching

1300 feet. The Lake City and Avon Park limestones found

below the Ocala are the principal aquifer used by

agricultural interests in central and south central Florida

and are known locally as the Floridan aquifer.

Above the Eocene series are the formations of the

Oligocene epoch. These are represented by the Marianna

limestone and the Byran limestones found and mined in the

Marianna Lowlands of the northern part of the state, and

the Suwanneelimestone found over the Ocala Formation as far

south as Hillsborough County.

The next higher formations are those of the Miocene

epoch. These are well represented by the Tampa limestones

of the early Miocene which are found above the Suwannee and

Ocala limestone in south Florida, the Chipola and Shoal

River Formations of the Alum Bluff group found in northwest

and north central Florida, the Hawthorn Formation, and the

Duplin marls. The latter three formations, Alum Bluffs,

Hawthorn, and Duplin, chiefly are sands, clays, and marls

that form a confining layer over the Eocene and Oligocene


The Hawthorn, with the possible exception of the

Ocala, is the most extensive formation within the state.

It occurs at, or near, the surface in most of north Florida.

It overlies the Tampa limestone formation in Hillsborough

County, and is, itself, overlain by the Duplin marls and

younger deposits in the south central and southern parts

of the state.


The surface material of most of the coastal lowlands

are of the Pliocene, Pleistocene, and Recent periods. The

most widely distributed are the sands formed along the old

shorelines of previous ocean levels. Cooke (1945) defines

seven of these marine terraces. Some small deposits of

coquina, oolite, coral reef limestone, and freshwater marls

are found among these deposits.



The purpose of this chapter is to enumerate some

of the more widely used modeling techniques in groundwater

flow, along with a brief description of each. The reader

is referred to Bear (1972) for additional information and



Numerical methods are used in many cases where the

partial differential equations governing flow through porous

media cannot be solved exactly. Various techniques have

been developed for obtaining numerical solutions.

Method of Finite Differences

The method of finite differences is one such

technique. The first step is to replace the differential

equations by algebraic finite difference equations. These

difference equations are relationships among values of the

dependent variable at neighboring points of the applicable

coordinate space.

The resulting series of simultaneous equations

is solved numerically and gives values of the dependent

variables at a predetermined number of discrete or "grid"

points throughout the region of investigation.

If the exact solution of the difference equations

is called D, the exact solution of the differential equation

is called S, and the numerical solution of the difference

equation is called N, two quantities of interest may be

defined. They are the truncation error, IS DI, and the

round-off error, ID NI. In order for the solution to

converge, it is necessary that IS DI 0 everywhere in

the solution domain. The stability requirement is such

that ID NJ ->- 0 everywhere in the solution domain. The

general problem is to find N so that IS NJ is smaller

than some predetermined error. Noting that (S N)

= (S D) + (D N), it is seen that the total error is

composed of the truncation error and the round-off error.

The arbitrary form selected for the finite difference

equation leads to the truncation error. This error is

frequently the major part of the total error.

The actual computation proceeds by one of two

schemes. They are the explicit, or forward-in-time,

scheme and the implicit, or back-in-time, scheme. The

explicit scheme is simpler but more time consuming than

the implicit scheme due to the stability constraint. The

implicit scheme is more efficient but requires a more

complicated program as compared to the explicit scheme.

Method of Finite Elements

The finite element technique employes a functional

associated with the partial differential equation, as

opposed to the finite difference method which is based on

a finite difference analog of the partial differential

equation. A correspondence which assigns a real number

to each function or curve belonging to some class is termed a


The calculus of variations is employed to minimize

the partial differential equation under consideration. This

is done by satisfying a set of associated equations called

the Euler equations. Thus, one seeks the functional for

which the governing equations are the Euler equations and

proceeds to solve the minimization problem directly, rather

than solving the differential equation.

The procedure is continued by partitioning the flow

field into elements, formulating the variational functional

within each element and taking derivatives with respect to

the dependent variables at all nodes of the elements. The

equations of all the elements are then collected. The

boundary condition is expressed in terms of nodal values

and incorporated into the equations. The equations are

then solved.

Relaxation Methods

This method may be applied to steady state problems

which are adequately described by the Laplace or Poisson

equations. The process involves obtaining steadily improved

approximations of the solution of simultaneous algebraic

difference equations.

The first step of the procedure is to replace the

continuous flow domain under investigation by a square or

rectangular grid system. The governing differential

equation is also replaced by corresponding difference

equations. Next, a residual, say R is defined corre-

sponding to the point o on the grid. Ro represents the

amount by which the equation is in error at that point.

If all values of the equation are correct, R0 will be zero

everywhere. In the initial step, values are assigned at

all grid points and, in general, the initial residuals will

not be zero everywhere. The process now consists in

adjusting values at each point so that eventually all

residuals approach zero, or to at least some required


The reduction of residuals is achieved by a
"relaxation pattern" which is repeated at different grid

points so as to gradually spread the residuals and reduce

their value.


As implied in section A, direct analytical solutions

are frequently inadequate or impractical for engineering

application. In many cases, the analytical solutions which

are found are difficult to interpret in a physical context.

In an attempt to circumvent some of the shortcomings of a

purely mathematical approach, model and analog methods are

frequently employed. The analog may be considered as a

single purpose computer which has been designed and built

for a given problem.

Modeling, then, is the technique of reproducing the

behavior of a phenomenon on a different and more convenient

scale. In modeling, two systems are considered: the proto-

type, or system under investigation, and the analog system.

These systems are analogus if the characteristic equations

describing their dynamic and kinematic behavior are similar

in form. This occurs only if there is a one-to-one corre-

spondence between elements of the two systems. A direct

analogy is a relationship between two systems in which

corresponding elements are related to each other in a

similar manner.

A model is an analog which has the same dimensions

as the prototype, and in which every prototype element is

reproduced, differing only in size. An analog is based

on the analogy between systems belonging to entirely

different physical categories. Similarity is recognized

in an analog by two characteristics: (1) for each dependent

variable and its derivatives in the equations describing

one system, there corresponds a variable with corresponding

derivatives in the second system's equations, and (2) inde-

pendent variables and associated derivatives are related to

each other in the same manner in the two sets of equations.

The analogy stems from the fact that the characteristic

equations in both systems represent the same principles of

conservation and transport that govern physical phenomena.

It is possible to develop analogs without referring to the

mathematical formulation; an approach which is particularly

advantageous when the mathematical expressions are exces-

sively complicated or are unknown.

Analogs may be classed as either discrete or

continuous with respect to space variables. In both cases,

time remains a continuous independent variable.

The need for complete information concerning the

flow field of a prototype system is obvious, and no method

of solution can bypass this requirement. However, in many

practical cases involving complicated geology and boundary

conditions, it is usually sufficient to base the initial

construction of the analog on available data and on rough

estimates of missing data. The analog is then calibrated by

reproducing in it the known past history of the prototype.

This is done by adjusting various analog components until a

satisfactory fit is obtained between the analog's response

and the response actually observed in the prototype. Once

the analog reproduces past history reliably, and within a

required range of accuracy, it may be used to predict the

prototype's response to planned future operations.

Sandbox Model

A reduced scale representation of a natural porous

medium domain is known as a sandbox model, or a seepage tank

model. Inasmuch as both prototype and model involve flow

through porous media, it is a true model.

A sandbox model is composed of a rigid, watertight

container, a porous matrix filler (sand, glass beads, or

crushed glass), one or several fluids, a fluid supply system

and measuring devices. The box geometry corresponds to that

of the investigated flow domain, the most common shapes being

rectangular, radial, and columnar. For one-dimensional flow

problems, the sand column is the most common experimental

tool. Transparent material is preferred for the box construc-

tion, especially when more than one liquid may be present and

a dye tracer is to be used. Porosity and permeability

variations in the prototype may be simulated by varying the

corresponding properties of the material used as a porous

matrix in the model according to the appropriate scaling

rules. The porous matrix may be anisotropic. In order to

measure piezometric heads and underpressures, piezometers

and tensiometers may be inserted into the flow domain of

the model.

Wall effects are often eliminated by gluing sand

grains to the walls of the box. This effect can also be

reduced by making the porous matrix sufficiently large in

the direction normal to the wall. Inlets and outlets in the

walls connected to fixed level reservoirs or to pumps are

used to simulate the proper boundary and initial conditions

of the prototype.

Water is usually used in models which simulate ground-

water aquifers, although liquids of a higher viscosity may be

used to achieve a more suitable time scale.

The sandbox model is used extensively because of its

special features which permit studies of phenomena related

to the microscopic structure of the medium such as hydro-

dynamic dispersion, unsaturated flow, miscible and immiscible

displacement, simultaneous flow of two or more liquids at

different relative saturations, fingering, wettability, and

capillary pressure. The capillary fringe in a sandbox model

is disproportionately larger than the corresponding capillary

rise in the prototype and, for this reason, the sandbox model

is usually used to simulate flow under confined rather than

phreatic conditions.

Hele-Shaw Analog

The Hele-Shaw or viscous flow analog is based on the

similarity between the differential equations governing two-

dimensional, saturated flow in a porous medium and those

describing the flow of a viscous liquid in a narrow space

between two parallel planes. In practice, the planes are

transparent plates, and the plates are usually mounted in

a vertical or horizontal orientation.

The vertical Hele-Shaw analog was selected for this

study because it is more appropriate for the prototype system

under investigation. Also, it is not possible to model a

free goundwater table or percolation in a horizontal model.

A detailed description of this analog is presented

in Chapter III of this study.

Electric Analog

Three types of electric analogs are powerful tools

in the study of flow through porous media. They are the

continuous electric analog, the discrete electric analog

and the ion motion analog.

Continuous Electric Analog

This analog takes two forms: the electrolytic tank

and the conducting paper analogs. The analogy rests on the

similarity between the differential equations that govern

the flow of a homogeneous fluid through a porous medium, and

those governing the flow of electricity through conducting


In particular, Darcy's law for flow in a porous

medium and Ohm's law for the flow of an electric current

in a conductor may be compared. Also, the continuity

equation for an incompressible fluid flowing through a

rigid porous medium may be compared with the equation for

the steady flow of electricity in a conductor. One concludes

from this comparison that any problem of steady flow of an

incompressible fluid having a potential may be simulated by

the flow of electric current in an analog.

Discrete Electric Analog

This analog also takes two forms: the resistance

network analog for steady flow, and the resistance-capacitance

network for unsteady flow.

In this analog, electric circuit elements are

concentrated in the network's node points to simulate the

properties of portions of the continuous prototype field

around them. The unknown potentials are the solution of

the problem, and they can only be obtained for those points

which correspond to the nodes of the analog network. The

discrete electric analog is based on the finite-difference

approximation of the equations to be solved; therefore, the

errors involved in the discrete representation are the same

as those occurring in this approximation.

The electric resistor corresponds to the resistance

of soil to flow through it, and capacitors are used at the

nodes to simulate storage capacity of the prototype.

Ion Motion Analog

This analogy uses the fact that the velocity of ions

in an electrolytic solution under the action of a DC voltage

gradient is analogous to the average velocity of fluid

particles under imposed potential gradients in a porous

medium. In this case, both electric and elastic storativi-

ties are neglected. The primary advantage of the ion motion

analogy is that, in addition to the usual potential distri-

bution, it permits a direct visual observation of the movement

of an interface separating two immiscible fluids. In ground-

water interface problems where gravity is involved, this

analog cannot be used. Scaling for the analog is based on

the similarity between Darcy's law and Ohm's law governing

the ion motion in an electrolytic solution.

Physically the analog consists of an electrolytic

tank having the same geometry as the investigated flow

domain. Inflow and outflow boundaries are simulated by

positive and negative electrodes, and two- and three-

dimensional flow domains may be investigated.

Membrane Analog

The membrane analog consists of a thin rubber sheet,

stretched uniformly in all directions and clamped to a flat

plane frame. The achievement of.equilibrium of various

forces and stresses in the membrane (caused by distorting

the frame or transversal loads) leads to the Laplace equation

and to the Poisson equation. The analogy is based on the

similarity between these two equations and the corresponding

equations that describe the flow in the prototype.

This method is applicable mainly to cases of steady

two-dimensional flow involving complicated boundary geometry

and point sources and sinks within the flow field.


Following Bear (1972), Table 2.1 is presented as a

summary of the models and analogs discussed in section B of

this chapter. In section A, the numerical methods discussed

are most likely to be carried out on a digital computer. It

is important for the investigator to examine both the cost

and the applicability of these various numerical and physical

methods to his particular case. An analog is usually pre-

ferred to a digital solution when the accuracy and/or amount

of field data is small. In many simple cases, the analog is

likely to be less expensive than a digital computer; whereas,

for large regions or unsteady three-dimensional problems, the

computer may be less expensive.

The Hele-Shaw model also has definite advantages when

demonstration of the saltwater intrusion phenomenon to a

public body, or other laymen involved in political decision

making, is considered. This type of model allows for direct

observation of the phenomenon without the numerical interpre-

tations used in the computer models.




Dimensions of field
Steady or unsteady flow

Sandbox Model Hele-Shaw Analog
Vertical Horizontal


Simulation of phreatic surface yes'
Simulation of capillary yes
fringe and capillary pressure
Simulation of elastic yes,
storage di
Simulation of anisotropic yes
Simulation of medium yes
Simulation of leaky formations yes
Simulation of accretion yes

Flow of two liquids with appr
an abrupt interface
Simultaneous flow of two yes
immiscible fluids
Hydrodynamic dispersion yes
Observation of streamlines yes,

or three

two two
both both
yes* no
yes no

Electric Analogs
Electrolytic RC Network Ion Motion

two or

three two or three

for two yes yes yes, for two yes
mensions dimensions
yes yes yes yes
kx5kz kx ky
yes5 yes5 yes yes


yes yes yes' yes
yes yes yes, for two yes
yes yes no6 no6

no no

for two I
nsions, near
parent walls
three dimen-

no no
yes no



kxk y





I Subject to restrictions because of the presence of a capillary fringe.
z By trial and error for steady flow.
I By trial and error for steady flow, or, as an approximation, for relatively small phreatic surface fluctuations.
By scale distortion in all cases, except for the RC network and sometimes the Hele-Shaw analog where the hydraulic
conductivity of the analog can be made anisotropic.
s With certain constraints.
6 For a stationary interface by trial and error.




kx y







The viscous flow analog, more commonly referred to

as a parallel-plate or Hele-Shaw model, was first used by

H. S. Hele-Shaw (1897, 1898, 1899) to demonstrate two-

dimensional potential flow of fluid around a ship's hull

and other variously shaped objects. The analog is based on

the similarity of the differential equations which describe

two-dimensional laminar flow, or potential flow for that

matter, of a viscous fluid between two closely spaced parallel

plates; and those equations which describe the field of flow

below the phreatic surface of groundwater, namely Darcy's


qx = K (
x x ax

qz = K (3.1b)


qx, qz = Darcy velocity of specific discharge in the
x-direction and z-direction, respectively

K K = hydraulic conductivity in the x-direction
Kx z and z-direction, respectively

x = horizontal direction (major flow direction)

z = vertical direction

= potentiometric head

and, by use of the conservation of mass principle, the

Laplace equation

2 + 822 = 0 (3.2)
ax2 az2

Viscous Flow Analog

To demonstrate the analogy of model and prototype,
the equations of motion and continuity for laminar flow of

a viscous fluid between two closely spaced parallel plates

will be developed and then compared to equations 3.1a,

3.1b, and 3.2.

Consider a viscous incompressible fluid flowing

ever so slowly between two parallel plates which are spaced

such that the Reynolds' number, R based on the interspace

width is less than 500 (Aravin and Numerov, 1965). In

Cartesian coordinates, the general Navier-Stokes equations,

i.e., the equations of motion, are

DV B + 1 I- + V2Vx (3.3a)

DVy = B + + V2V (3.3b)
Dt y+ y ]

Vz = B + 1 3P + v2V (3 .3c)
Dt z p DZ


D = substantial derivative =- + V a
Dt at x ax
+ V a + V z -

V2 = Laplace operator = 2 + y2 +

Vx, V Vz = velocity in the x-, y-, and z-directions,

B B Bz = body forces in the x-, y-, and
xz z-directions, respectively

p = pressure

p = density of the fluid

p = absolute viscosity of the fluid

y = horizontal direction (minor flow

t = time

Referring to the free-body diagram of the idealized
flow regime shown in Figure 3.1, if no slip conditions

(adherence to the walls) of the fluid particles are assumed

in the molecules closest to the walls of the parallel plates,

it is easily seen that the velocity gradient in the

y-direction is much larger than the velocity gradient in

either the x- or z-directions. Thus., the first and second

order partial derivatives taken with respect to both x and

z may be neglected when compared to those taken in the

y-direction. Secondly, because of the very low velocities

("creeping" motion) the intertia terms, that is the terms


Free body flow diagram for Hele-Shaw model.



_ _II~



on the left side of equations 3.3, are very small when

compared to the viscous terms, those on the right side of

equations 3.3, and may be neglected. Thirdly, because of

the restriction to two dimensions, the velocity in the

y-direction is taken to be zero; consequently, all rates

of change of velocity in the y-direction.must be zero.

Finally, the only noncancelable body force acting on the

fluid is gravity which acts only in the vertical. Mathe-

matically, Bx -x (gz) = 0; By = (gz) = 0; and
Dz -~ a (gz) = -g = -32.17 ft./sec.2. Incorporating all

of the above arguments and values into equations 3.3, the

equation of motion becomes:

3+ 0 (3.4a)

2= 0 (3.4b)

pg a- = 0 (3.4c)

Defining the potentiometric head, or potential, i = z + p/y,

where y equals the unit weight of the fluid, and taking the

partial derivative with respect to x, y, and z, the following

results are obtained after multiplying through by the unit

weight of the fluid:

Y M =_ ap (3.5a)
Y ax ax

Y -- = DP (3.5b)
Tay ay

Y (y y 'I

Introducing these relationships into equations 3.4 and

dividing through by the unit weight yields:

3t = P2 x (3.6a)
3x y ay2

D| 0 (3.6b)

S_ 1 2 z (3.6c)
2z y ay2

It is evident from equation 3.6b that the potentiometric head

is constant in the y-direction. It is possible then to inte-

grate the first and third equations of equation 3.6 with

respect to y. After separating variables and integrating
once, using the boundary condition y = 0, -x = 0, and
3V ay
z = 0, the following equations are obtained:

y a 1 x (3.7a)
a x y 8y

y = z (3.7b)
Yz y ay

Integrating, once again, using the second boundary condition

y = b/2, Vx = 0, Vz = 0 (no slip) the above becomes, after

solving for the respective velocities:

Vx y2 ) (3.8a)

Vz 2 (2 ] (3.8b)

Note, where b is the spacing between the plates, that if a
potential = y- b-4 is defined, equations 3.8 can
be written:

V (3.9a)
x ax

V (3.9b)

P, the velocity potential, is dependent only on y. Inte-
grating the velocity profiles established by equations 3.8
between the limits of b/2, and dividing by b, the
directional specific discharges are obtained:

q +b/ 2 y= b/2
qx b _/2 x y b = x N -b/2

=- b2 Y (3.10a)

qz b j-b/2 vz dy 2-i z (3 j-b/2

b2 Y (3.10b)
TZ T az

It is obvious that for a model of constant spacing b, the

quantity b2 does not vary in either the x- or z-direction.

Defining the model hydraulic conductivity as K = K
= b* ~ equations 3.10a and 3.10b become:

qx = Kxm x (3.11a)

qz = Km (3.11b)

which, of course, is analogous to equations 3.1.

Consider, now, the two-dimensional continuity

equation for flow between parallel plates:

aV 3V
_x + z 0 (3.12)
ax az

The specific discharge, or Darcy velocity, is

related to the velocity by the vector equation n q = V,

where ne is the effective porosity of the flow media. In

the model, ne equals 1. From the analogy Vx = q x; Vz = qz'

and substituting the relationship obtained from equation

3.11 into equation 3.12:

x- m a-x x- Kzm = z (3.13)

or dividing by -Km and recalling that for a model Km = Kxm

= Kzm

9 + 2i-- (3.14)
ax2 az2 (3.14)

which is clearly analogous to equation 3.2.

The similarity of equations 3.1 and 3.11 and
equations 3.2 and 3.14 establish the analogy.


The two-dimensional equation along the free surface,
or water table, of an anisotropic porous media given by Bear

(1972) is:

K p .2 + K 2 = n (3.15)
xp xp zp z az ep at
Sj p p

where the subscript p denotes the prototype. For a Hele-

Shaw model, using the subscript m, the equation can be

written as:

2 Mjt 2
K Im + K =- n mm (3.16)
xm jx-m- zm azm a m em atm

Introducing the similitude ratios, denoted by the

subscript r, of the corresponding parameter of model and


Kxr = K (3.17a)
Kzr = (3.17b)

Xr -x (3.17c)
zr = z (3.17d)

r -I (3.17e)
ner = n (3.17f)
tr = t-- (3.17g)
and substituting these relationships into equation 3.15, the

following is obtained:

Kxm [( m ) 2 Kzm [*a(m/r) 2 (m r)
T Tv K+ _zz
Kxr m/x Kzr Zm/Zr) (Zm/Z

nem < tm/tr) (3.18)

The ratios of model to prototype quantities are constant and
can be removed from behind the differential; therefore,
equation 3.18 can be rearranged in the following way:

Kxr r2 xm 2 Z r m 2
K r 2 xm ax m zm K ^ 2 3z

Zr m'I tr m
Kzr r -zmJ nerzr nem tm (3.19)

Comparing the equations 3.16 and 3.19, it is evident that, if

the equations are identical, the following must be true:

x 2 z 2 z t
1 K r r r r (3.20)
xr r zr r zr'r er r

Solving the third equality for zr, the following important

relationship is found:

zr =r (3.21)

The second equality, after cross-multiplying, yields:

S xr K (3.22)
zr Kzr

Recalling the definitions of Kxr and Kzr (equations 3.17a

and 3.17b) and remembering that K = K in an isotropic
xm zm
model, the above equation can be rewritten:

x r 2 K/Kxp K
r xm xp zp (3.23)
r K zm zp Kxp

The ratio of zK is called the ratio or degree of anisotropy
of the prototype.


Using the fourth equality of equation 3.20, the time

ratio of the model and prototype is established:

n z
t er r (3.24a)
r K

n x 2
t er r (3.24b)
r K rr

Substituting the vertical ratio, zr, for the potentiometric

head established by equation 3.20 and the similitude ratios

of time, hydraulic conductivity and porosity into equation

3.24b results in

tm nem X r (3.25)
tp nep Kxm z

The effective porosity of an isotropic model, nem, is unity.

The hydraulic conductivity of the model was defined pre-

viously as p?, thus the time scale for the model is finally

written as:

K x 2
t 12 vK xx r t (3.26)
m g ep zr P


The Hele-Shaw model is normally isotropic. This is

because of the nonvariance of the. spacing of the parallel

plates. There are, however, two methods for simulating

anisotropy in a model. Equation 3.23 gives a clue as to

the first possibility of simulating anisotropy:

r xm x_ zp (3.23)
r zm /Kzp xp

Since Km = Kzm, the x or z ratio can be adjusted so that the

model's hydraulic conductivities are kept equal. This is

usually done by choosing a suitable horizontal ratio. Know-

ing the prototype parameters, a vertical scale for the model

is computed so that the aforementioned conductivities are

kept equal, demonstrating:

zr zm X m (3.27)
p zp p

Solving for zm:

K x
zzm p xm z (3.28)
zp p

Unfortunately, the geometric distortion method is adequate

for modeling only one ratio of anisotropy. If there is a

second aquifer, within the prototype which has a different

vertical or horizontal hydraulic conductivity, the second

aquifer cannot be correctly simulated unless, of course,

the second aquifer's ratio of anisotropy is the same as the

ratio of the first. This restriction would severely limit

the use of the Hele-Shaw analog in modeling of regional

groundwater problems unless another method were available

to correct the ratio of anisotropy.

Polubarinova-Kochina (1962) suggests using a grooved

plate within the model to correct the ratio of anisotropy of

the second flow zone. The plate may be grooved in any one

of several methods. It matters little whether a grooved

plate is sandwiched between the parallel plates, or if

rectangular bars are attached to the front or back plate.

The degree of anisotropy of the second aquifer and the amount

of geometric distortion used to model the first flow zone

determines the directions in which the grooves, or bars,

are placed; however, the grooves are normally placed

horizontally or vertically. Collins and Gelhar (1970) have

developed the conductivity equations for the flow zone in

which Polubarinova-Kochina's grooved plate is used. The

analysis assumes one-dimensional flow and can be used

equally well with either vertical or horizontal orientation

of the grooves.

Following Collins and Gelhar (1970), consider flow
in a grooved portion of a model. For simplicity, assume

Figure 3.2 is a section of the grooved zone. Assuming such,

the horizontal direction then corresponds to the x-direction

and the grooves, which are vertical, lie in the z-direction.

Area 1 is associated with the wider spacing of length ab.

Area 2 is associated with the narrower spacing of length b.

Since flow area 1 is the much larger of the two areas, most

of the frictional head loss occurring through the total

length L is developed in flow area 2 which has length

(1 M)L. Lambda, X, is a length correction factor..

Referring to Figure 3.3, the potentiometric gradient

ax across area 2 is:


Section of anisotropic grooved zone in Hele-Shaw model.

Back Plate


AL Front Plate
L (1 A)L

Plan View

Perspective View


Head loss in grooved anisotropic zone.

A 1-)

x or z



--- --


ax (1 A)L (3.29)

For high values of a:

___ A = 2C(3 .30)
ax L L

but, from equation 3.29, 2 = (1 A) 3 so that:

a = (1 ) ~2 (3.31)
ax ax

Applying Darcy's law to area 2:

q= K (3.32)
qx x2 ax

and, substituting the previous expression for :2-
S- x2 1 (3.33)
x (1 A) ax

The effective hydraulic conductivity in the x-direction

then is:

K Kx b2g (3.34)
xm (1 X) 12v 1 ) (334)

Consider vertical flow through the grooved zone

illustrated in Figure 3.2. In particular, consider flow

downward through areas 1 and 2. The total discharge

through these areas can be written as the sum of the dis-

charge through each, that is, Qzm = Q1 + Q2. Applying

Darcy's law for the total discharge, Q:


Q, = K (abAL) (3.35)

Q2 = K ( ) Lb (3.36)
z2 az

Adding Qi and Qz:

Qm= K (aX) + K (1 )]bL (3.37)
zm z1 z2z z

For flow area 2 it is not unreasonable to assume that the

frictional forces in the fluid boundary to either side of

area 2 are negligible. Therefore, the vertical conduc-

tivity in this area is the same as defined by the earlier

analysis, that is:

K b2 = b2 (3.38)
Z2 = l- 12

where v is the kinematic viscosity. Furthermore, if b/ab

<< 1, the flow in area 1 can be assumed roughly equivalent

to flow through a rectangular hole. According to Rouse

(1959), the equation of motion through a rectangular

cross section of length XL and width ab is given by:

v x(x XL) + X sin JTx
z 2i 3z jl AL

x A. cosh 1J + B. sinh (3.39)



A. 2y(AL) 2 (cos j I) (3.40)


cosh jab 1(
B. = A. (3.41)
sinh jab3

By integrating Vz over the area and dividing by the total
area ALab, the mean velocity is given by:

V= (a b)2 (3.42)
z az


S= 192 b j (cos (j ) 1)2 _
= 2 -L j( 4j5 2cab


Since the terms of the infinite series decrease as j5s, only
the first term of the series need be considered and retained,
so that:

1 19 tanh fXL (3.44)

Additional information may be found in Appendix A. From
equation 3.42, the equivalent hydraulic conductivity in
area 1 is given by:

K Y (ab)2 (3.45)
Kz P 12

Finally, introducing the values found for K and K into
zE Z2
equation 3.37:

Qzm =-] (b) (aX) + b (1 U) bL (3.46)


Qzm = (3 + 1 A) bL (3.47)

from which it is seen that the effective vertical hydraulic

conductivity is given by:

Kzm = (a 3X1 + 1 X) (3.48)

after defining Qzm/bL = Vz, where Vz = qz is the effective

vertical specific discharge. Equations 3.34 and 3.48 give

the second method available to correct the hydraulic con-

ductivity of a model so that it can simulate the true

ratios of anisotropy found in the prototype.


An aquiclude can be defined as a soil stratification

in which the hydraulic conductivities are zero. In certain

geohydrologic problems, it is convenient to assume such

conditions. However, in reality few soil masses are truly

impervious. The degree of perviousness in a stratum is

referred to as leakance and it is generally assumed that the

direction of flow is only vertical. There is no horizontal

flow, that is,

Kxp = 0 (3.49)

Bear et al. (1968) suggest the use of vertical slots

to model such a semipervious layer. To accomplish this, the

spacing between the parallel plates of the Hele-Shaw analog

is filled with a slotted middle plate. See Figure 3.4.

The analysis to determine the effective vertical

hydraulic conductivity of a model's leaky layer closely

parallels that for an anisotropic grooved zone. Again,

following Collins (1970), Darcy's law for flow through a

vertical slot is:

QZ = ALab (3.50)
z = z 3z

The effective specific discharge through the slot found by

integrating the Rouse equation (equation 3.39) is the same

as equation 3.42 from which is found the hydraulic conduc-


Kz (b)2 (3.51)

and, introducing the above into equation 3.50:

Qz = (b)2 Lab (3.52)
z U -2-- aL z


Section of leaky zone in Hele-Shaw model.


Back Plate
Solid ab
x l b

F-I v
PL (1-X)L t
Front Plate L

Plan View

Perspective View

Again, the effective specific discharge, or mean velocity,

is equal to:

z 3 b (3.53)
Vz bL 12 az 353)

so that the effective hydraulic conductivity of a leaky layer

in the model is:

Kzm a b2 (3.54)


While the problem of storage has not been completely

solved, it has, in general, been neglected by most researchers.

Bear (1960) suggests that discrete tubes attached to either

the front or back plate and connected to the aquifer be used

to model the specific storage of a confined aquifer. For a

nonisotropic aquifer, the right hand side of equation 3.13

is not zero, but, in fact, equals the specific storage, S ,

times the rate of change of the potentiometric head. Rewrit-

ing the two-dimensional equation 3.13 for both model and

prototype to include the above gives the following:

K -2 + K --- = S -- (3.55)
xp ax 2 zp aZ 2 op at
p p p

K 2m + K m m (3.56)
xm Dxm2 zm 3Z%2 om atm
m m m

Defining a ratio of storativity:

Sor = (3.57)

it follows from inspection that,

z K z
K Sr- -- (3.58)
xr x2 zr or t

or that,

K t S
zr r om (359)
or z S (3.59)
r op

Referring to Figure 3.5, the storage represented by

the model in the discrete length 1m is equal to:

S om_ m (3.60)
mo m b 1 z

where A is the cross-sectional area of the storativity tube.

Introducing the above into equation 3.59 and solving for Am

A bl z Sop K r (3.61)
m mmm op zr z 2

The discharge scales are.obtained from Darcy's law.

Written for both prototype and model with the usual sub-

scripting, these are in the x-direction:

Q K -- b z (3.62)
xp = xp ax pp


Storage manifold for Hele-Shaw model.

Area, A
bm m



End Section





xm -xm 3x mm

Dividing equation 3.63 by equation 3.62 and recalling the

definitions for the various parameters' ratios, it follows


z z 2
S K r b r=K b r
xr xr r r x xr rx

Similarly, in the z-direction,

Q = K r b r K b x
zr zr r r zr zr r r



Solving equation 3.22 for the hydraulic conductivity in the


Kx K
xr zr




and, substituting


Qxr = Kzr

this result into equation 3.64, it follows

r b r K b x
r r


Qxr = Qzr = Q
*xr -zr -r




Accretion, R, is the rate at which a net quantity

(precipitation and surface inflow minus evapotranspiration,

runoff, etc.) of liquid is taken into the flow system at

the phreatic surface. It is measured as a volume per unit

horizontal area per unit time, that is:

Rr (3.69)

From equations 3.64 or 3.65, it follows that:

K z 2
R x Kzr (3.70)


On occasion, volume, U, is of some importance. The

volume scale follows directly from continuity, that is:

Ur = Q tr (3.71)

Substituting the values found from equations 3.65 and 3.24a

for Qr and tr, respectively, the above equation becomes:

r zr r r er K r er r r (3.72)

As inferred earlier in this section's opening sentence, the

volume scale is usually neglected; however, in the case of

free surface water bodies, lakes, rivers, etc., if the volume

exchange of liquid is of interest and has to be modeled, the


volume scale requires an additional restriction. In the

following analysis, the bar above the width dimension

indicates the free water surface of a river, lake, ocean,

or such.

In the portion of the model simulating the body of

water, the spacing of the model is increased to maintain

hydrostatic pressure distributions within the model. The

narrower spacing of the model is, of course, a measure of

the hydraulic conductivity of the aquifer. In the proto-

type, however, the width of the open water and the aquifer

are equal and this leads to the following (Bear, 1960) for

the model and prototype, respectively:

U = n b x z (3.73a)
r er r r r

U = n b x z (3.73b)
r er r r r

The same volume ratio must be applicable to both the narrow

and the enlarged interspace; therefore, Ur Ur. It follows


n er b = ner b (3.74)


n em 1 (3.75)

br = ner br (3.76)

Note that for an anisotropic media, nem does not necessarily

equal one.


Site Selection

The site selected for this study is the middle Gulf

area of Florida. This region has a rapidly expanding popu-

lation with a corresponding growth in water demand. The

increased pumping to satisfy this demand also increases

the likelihood of saltwater intrusion, and, in fact, a

number of municipal supply wells in the coastal zone have

been shut down in recent and past years due to chloride


Black et al. (1953) list eight factors responsible

for saltwater intrusion. They are:

1. Increased water demands by municipalities.

2. Increased water demands by agriculture.

3. Increased water demands by industry.

4. Excessive drainage.

5. Lack of protective works against tidewater in
bayous, canals, and rivers.

6. Improper location of wells.

7. Highly variable annual rainfall with insufficient
surface storage during droughts.

8. Uncapped wells and leakage.

Of these eight factors, numbers 1, 2, 3, 6, and 7 would apply

to this area. The city of St. Petersburg is an outstanding

example of these problems. Their original water supply was

from local artesian wells, but increased demands caused salt-

water intrusion and forced the closing of these wells. In

1929 the present Cosme-Odessa field was located farther

inland to escape this problem.

One of the major water supply systems in this region

is the Pinellas County Water System, and this study is

centered around the Eldridge-Wilde well field of this system.

The location of Eldridge-Wilde in relation to several of the

population centers of this region is shown in Figure 4.1.

It is about 8 miles east of the Gulf of Mexico and encom-

passes an area in the northeast corner of Pinellas County,

at the intersection of the boundaries of Pinellas, Hills-

borough, and Pasco counties.

This system was instituted in 1937 to supply the

towns along the Gulf coast from Belleair Beach to Pass-a-

Grille. Its original form was that of a raw water reservoir,

and the first wells were not drilled until 1946 in the

McKay Creek area. These wells were soon contaminated with

salt water, and investigations were begun in 1951 to locate

the well field at its present site. This well field has

grown over the years, and in 1970 (Black, Crow and Eidsness,

Inc., 1970), the waterworks facilities at Eldridge-Wilde

included: sixty-one water wells, over 11 miles of raw water


Regional area of prototype.

il MEXIUI0GQ;:1;:.: ,

iiii -- / ko
V -,^ / 10'

Dade City

60' /


, \



^ ./ I


N "- 50'

S 30'

Eldridge-Wilde Well Field

Lake Tarpon
0O Tampa

:: :1 :i:: Peter .sbug : ::? :::: ::::::l

::.. : : ... .. : :::::: ::::: :

...................i 9

collection piping, water treatment facilities consisting of

aeration and chemical treatment, including chlorination and

fluoridation, and high service pumping units.

All wells are open hole and penetrate the Floridan

aquifer at depths from 140 to 809 feet below ground surface,

averaging 354 feet. The design capacity, of the field at the

present time is 69 million gallons per day, although the

maximum allowable pumpage has been set by the Southwest

Florida Water Management.District at 28 million gallons per

day on the average with a maximum day of 44 million gallons

per day.

In selecting the prototype location within the site

area, two characteristics of the vertical Hele-Shaw analog

must be considered. The first characteristic is that there

can be no general flow normal to the parallel walls of the

model. This means that the flow from one end of the model

to the other is streamline flow. The second characteristic

is that the ends of the model are finite. Therefore, the

prototype must be along a streamline in the flow domain and

have boundary conditions which are "infinite" reservoirs or

water divides.

The prototype selected meets the above requirements

and includes the point of interest, i.e., Eldridge-Wilde

well field. The center line of the prototype is shown in

Figure 4.1 as the unbroken line passing through Eldridge-

Wilde in a southwest to northeast direction. The dotted

contours in the figure define the potentiometric surface of

the Floridan aquifer in feet above mean sea level as of May,

1971. They were obtained from a map publication entitled

"Potentiometric Surface of Floridan Aquifer Southwest Florida

Water Management District, May, 1971" prepared by the U. S.

Geological Survey in cooperation with the Southwest Florida

Water Management District and the Bureau of Geology, Florida

Department of Natural Resources. Now, in a flow field,

streamlines are perpendicular to potentiometric lines. As

can be seen from the figure, the prototype orientation

reasonably satisfies the streamline requirement. The proto-

type is terminated on the southwestern end at the 15 feet

depth contour in the Gulf of Mexico, and it is assumed that

this satisfies the infinite reservoir boundary condition.

The northeastern terminus is located in the center of the

80 feet contour, southwest of Dade City. This location

satisfies the water divide boundary condition. The area in

the vicinity of the 80 feet contour is known as the Pasco

High. The overall length of the prototype is 36 statute

miles. The width of the prototype is taken to be 3.5

statute miles. This dimension is sufficient to include the

cone of depression caused by pumping in Eldridge-Wilde well

field, and is based on the results of a study by Mr. Evans

(employing a numerical model) for Black, Crow and Eidsness,

Inc. The land surface contours of the prototype were

obtained from U. S. Geological Survey topographic maps.

The bottom boundary of the prototype is taken to be the base

of the Lake City Formation, with depths being determined

from available well logs of wells in the prototype vicinity.

The maximum depth from highest land surface to deepest point

is 1340 feet.

Prototype Geology

Stewart (1968) identifies eight formationsas being

of interest in terms of water production in the prototype

area. They are in descending order, the Undifferentiated

Deposits, Tampa Limestone, Suwannee Limestone, Crystal

River Formation, Williston Formation, Inglis Formation,

Avon Park Limestone, and Lake City Limestone. Underlying

the Lake City Limestone is the Oldsmar Limestone which is

not used as a source of water at present.

The Undifferentiated Deposits are interbedded sand,

silt, and clay of Post-Miocene age and range in thickness

from zero near the Pasco High to 60 feet in the Eldridge-

Wilde well field. The thickest deposits are in northeast

Pinellas County around the north end of Lake Tarpon where

sand dunes, as much as 40 feet high, overlie alternating

layers of clay, thin limestone beds, and sand greater

than 70 feet. thick.

The Tampa Limestone is a hard, dense, sandy, white

to light tan, or yellowish-tan fossiliferous limestone of

Miocene age. This limestone is near the surface in the

area of the Pasco High and about 80 feet below land surface

at the Eldridge-Wilde well field. At Eldridge-Wilde, the

thickness varies erractically from about 20 to 240 feet.

The Tampa Limestone is a poor to fair producer of water.

The Suwannee Limestone is a soft to hard, nodular

or grandular, fossiliferous white to tan limestone of

Oligocene age and is about 200 feet thick. The Suwannee

and Tampa Limestones are the major water producers for

wells in the area.

The Crystal River, Williston, and Inglis Formations

comprise the Ocala Group of late Eocene age. The Crystal

River and Williston Formation are lithologically similar

units of white to cream, porous, soft, coquinoid limestone

and are generally poor producers of water. The Inglis

Formation is a hard, cream to brown to gray fossiliferous

limestone and is generally a good producer of water.

The Avon Park and Lake City Limestones are litho-

logically similar units of soft to hard, cream to brown,

fossiliferous limestone with beds of dolomitic limestone

and some gypsum. Both formations are good producers of

poor quality water.

The Oldsmar Limestone is a fragmental dolomitic

limestone with lenses of chert, thin shale beds, and

some gypsum.

In this study, two formations are considered,

the Undifferentiated Deposits and the Floridan aquifer.

The Floridan aquifer is considered to contain all formations

from the Tampa to, and including, the Lake City Limestone.

The transmissivity of the Floridan aquifer ranges

from about 165,000 to 550,000 gallons per day per foot,

and the coefficient of storage ranges from about 0.0005

to 0.0015. The coefficient of leakage is approximately

0.0015 gallons per day per cubic foot.

Based on groundwater discharge and water levels,

the estimated recharge (leakance) to the Floridan aquifer

was computed (Stewart, 1968) to be about 103 million

gallons per day. Based on aquifer test data, the estimated

recharge for a 250 square mile area was 90 million gallons

per day.

The Undifferentiated Deposits act as a confining

layer, and the Floridan aquifer is thus under artesian


Prototype Hydrology

The surface waters of the area consist of many

lakes and few streams. Because of the flat topography,

little water runs off into streams, and swampy wetlands

are numerous. Most rainfall evaporates or is transpired

by plants.

The Floridan aquifer is recharged through the

Undifferentiated Deposits by surface and groundwater

derived from local rainfall. Many millions of gallons

of water are also admitted to the aquifer by numerous sink-

holes in the region. Water levels in the Floridan aquifer

respond to rainfall since this is the recharge source.

This response is not immediate, but usually fluctuates with

the wet and dry seasons. Water levels in wells which are

not directly affected by local pumping show yearly lows in

the dry season, April and May, and yearly highs during the

wet season, late summer or early fall (Black, Crow and

Eidsness, Inc., 1970).

The aquifer recharge has been estimated (Black,

Crow and Eidsness, Inc., 1970) from available data and the

use of the following formula:

Aquifer Recharge = P + SWI + GWI ET R GWO

The basin area is 575 square miles, changes in

storage are assumed zero, and evapotranspiration is assumed

to be 75 percent of the precipitation. The applicable

values are listed below in million gallons per day:

P = Precipitation

SWI = Surface Water Inflow

GWI = Groundwater Inflow

ET = Evapotranspiration

R = Runoff

GWO = Groundwater Outflow
Aquifer Recharge

This value is in reasonable

reported values (Stewart, 1968).

= +1492

= + 0

= + 0

= -1119

= 218

= 37
= + 118

agreement with previous





The selection of the prototype area was discussed

in Chapter IV. Table 5.1 is a summary of the prototype

characteristics. The leaky layer is synonymous with the

undifferentiated deposits. The top and bottom of the

Floridan aquifer were determined by straight-line extrapo-

lation from available well logs.

The hydraulic data are within the reported range

of values and are the result of a trial and error process

to stay within the range and still-produce a reasonable



The purpose of the Hele-Shaw analog in this study

is to model saltwater intrusion." Before discussing the

model design, it seems appropriate at this point to pro-

vide some background about saltwater intrusion. Water,

in general, whether it be surface water or groundwater,

is continually migrating towards the sea, where an

equilibrium, or moving freshwater/saltwater interface,



Parameters Aquifer Leaky Layer


x (ft.)

z (ft.)

yp = bp (ft.)


Txp (gpd/ft.)

Tzp (gpd/ft.)

Kxp /Kzp

Leakance (gpd/ft.3)

Kzp (gpd/ft.2)

Kxp (gpd/ft.2)


v @ 77* F (ft.2/sec.)

gp = 32.2 ft./sec.2










0.965 (10-5)








0.965 (10-5)

is established. The two fluids are miscible, but because

of the difference in densities and the very low velocities,

the interface is formed. Across the interface, the salinity

varies from that of the fresh groundwater to that of the

ocean. The transition zone, as it is called, is due to

hydrodynamic dispersion and, although it is anything but

abrupt, it is usually assumed to be. The interface then

is generally selected to occur at some measured electric

conductivity or salt (chloride) concentration.

The earliest investigations of saltwater encroch-

ment were made by Badon-Ghyben (1888) in Holland and

Herzberg (1901) in Germany. Working independently, both

investigated the equilibrium relationships between the

shape and position of the freshwater/saltwater interface.

Figure 5.1 shows a coastal phreatic aquifer and the

Ghyben-Herzberg interface model. Badon-Ghyben and Herz-

berg assumed static equilibrium and a hydrostatic pressure

distribution in the fresh groundwater and stationary saline

groundwater near the interface.

Considering a point P on the interface, and choosing

mean sea level as the datum, the pressure at point P is:

pp = h 5 (5.1)


h = vertical distance from mean sea level to
s point P

s = unit weight of sea water


Ghyben-Herzberg interface model.


Water Table


h h

pdcA Interface

Y^ \pdA Y

This pressure may, also, be expressed by:

p = {hf + hs] f (5.2)

in which hf equals the vertical distance from mean sea level

to the phreatic line at the location of p and yf equals the

unit weight of freshwater. Equating the, preceding two


hSys = hfyf + hsf (5.3)

and, rearranging,

hsYs hf = hff (5.4)

Solving for hs :

hs Y= Yf hf (5.5)

Introducing y = pg, where p is the density, factoring and

canceling out the gravity term, the Ghyben-Herzberg

relation is found:

h = -f hf (5.6)
S P-S f

for a saltwater density of 1.981 lb sec.2/ft.4 and a fresh-

water density of 1.933 lb sec.2/ft.4, @ 250 C, the quantity

pf/ps pf) = 40. The implications of equation 5.6 are
rather dramatic. For instance for every foot of freshwater

above the datum, there is 40 feet of freshwater below

the datum. More importantly however, consider the effects

of lowering the phreatic surface. For every one foot drop

of the water table, the interface raises 40 feet. It must

be remembered that the above analysis assumes static

conditions. This, in fact, is not always the case. The

position of the interface is a function of dynamic conditions

rather than static. Even so, in cases where flow is quasi-

horizontal, i.e., the equipotential lines are nearly

vertical, equation 5.6 is valid.

Many investigators have incorporated dynamic forces

into the analysis of the stationary interface. Hubbert

(1940) was able to ascertain a more accurate determination

of the shape of the interface near the coast line. He

assumed that at the interface the tangential velocity was

zero in the saltwater, but increases with horizontal dis-

tance in the freshwater as the coast line is approached.

This then is the cause for the interface to tilt upwards

as the sea is approached and the greater depths found than

those estimated by the Ghyben-Herzberg relationship.

Hubbert showed that the Ghyben-Herzberg equation holds

between points on the water table-and the interface along

an equipotential line, rather than along a vertical plane.

R. E. Glover (1959) modeled an infinitely deep

coastal aquifer by assuming no flow in the saltwater

region, a horizontal water table and a horizontal seepage

face located seaward of the coast line. He found an exact

solution for the shape of the wedge, giving the following


z2 2qx + Pf 2 (5.7)
K s f K2 S f-
Pf Pf

where x and z are the horizontal and vertical directions,

respectively, q is the seepage rate per unit width and K

is the hydraulic conductivity. De Wiest (1962) using

complex variables and a velocity potential of D = Kx .

p Pf)/pf derived the same equation.

Bear and Dagan (1964b), using the Dupuit assumptions

and the Ghyben-Herzberg equation, developed the approximate

shape of the interface for a shallow aquifer of constant


All of the above investigated the equilibrium

position of the saltwater/freshwater interface.

If there is a change in the freshwater flow regime,

a transition period is caused during which the interface

moves to a new point of equilibrium. The nonlinear

boundary conditions along the interface make the solution

for the shape and position of the transient interface all

but impossible except for the simplest geometries. Bear

and Dagan (1964a), as well as other investigators, have

used the Dupuit assumptions to approximate the rate of

movement of an interface in a confined aquifer. Following

Polubarinova-Kochina's (1962) suggestion, they assumed

quasi-steady flow and were able to approximate the inter-

face shape and position for both a receding motion and

landward motion of the interface.

Characteristically, the solutions obtained by

investigators to date have all had simple,geometries and

involved simplifying assumptions, some of which have had

little resemblance to actual conditions. Therein lies the

advantages of a Hele-Shaw analog, complex geometries and

boundary conditions can be modeled with relative ease.

In order to satisfy additional similitude require-

ments for the flow of two liquids with an abrupt interface,

as in this study, and to provide a suitable time ratio,

two liquid silicone fluids were chosen to be used in the

model. Dow Corning Corporation Series 200 silicone fluid

was used to model fresh water. Series 510 silicone fluid

from the same company was used to simulate salt water.

The 200 Series fluid and the 510 Series fluid have densi-

ties of 0.977 gm/cm3 and 1.00 gm/cm3, respectively, at

250 C. Both fluids have kinematic viscosities of 500

centistokes at 250 C. Dow Corning 200 fluid is a clear

dimethyl siloxane which is characterized by oxidation

resistance, a relatively flat viscosity-temperature slope

and low vapor pressure. Dow Corning 510 fluid is a clear

phenylmethyl polysiloxane which also has a relatively flat

viscosity-temperature slope. In order to locate and follow

the interface movement, the denser fluid was dyed blue and

the lighter fluid dyed orange.

The dimensions of the model were selected so that a

unit of reasonable size would be produced. These dimensions

are 11.75 feet long, 0.5 inch wide inside, and 2 feet deep.

The parameters xr, Zr, and br are therefore set, and the

result is a distorted model. This distortion requires the

use of a slatted inner zone as discussed in Chapter III.

Tables 5.2 and 5.3 list the applicable parameters.

The following analysis is shown for the design of

the model Floridan aquifer, leaky layer and storage


A. Slatted Anisotropic Zone for the Floridan Aquifer.

1. Compute k xm/kzm from equation 3.22

x 2 Kxr
r _= xr (3.22)
Sr Kzr

Noting that k = K v (5.8)
Then 3.22 becomes,

x '2 k k
zr XM Z kz (5.9)
r ixp zm

and finally,
p [ -6.182 10-0|2
xMm = (1.22) 0 [1 = 0.00209
zm Kzp zr1.493 x 1



Parameters Aquifer Leaky Layer


x (ft.)

z (ft.)

ym (ft.)



b (in.)

L (in.)

ab (in.)

XL (in.)

(1 A)L (in.)


kxm (in.2)
kzm (in.2)

xm /kzm
S (ft.-1)

vm @ 77* F (ft.2/sec.)

gm = 32.2 ft./sec.2












3.94. x 10-4




5.382 x 10-3













1.05 x 10-4


5.382 x 10-3




This is the value which the slatted zone must


2. Select dimensions for the slatted zone as shown

in Figure 3.2, as follows:

hold ab = 0.5 inch

select b

find a

select (1 X)L

select XL

find L

find x

3. Compute kxm b2 ( 1 (3.34)

4. Compute kzm = ((3AE + 1 A) (3.48)

5. Compute kxm/kzm and compare to the results of

step (1).

6. Repeat the process until the result of step (5)

equals the result of step (1).

B. Slotted Leaky Layer. Refer to Figure 3.4 for this


kxlm = 0 (5.10)

b2 (3.54)
kzlm = i-Z a (3.54)


Also, continuity of flow between the leaky layer and

the Floridan aquifer requires that:


Qlr = Qfr


Q = Qzr = Qxr = K bzrb x
"r zr "xr zr rr


K b x =K b x
zlr rr zfr rr


Kzlr = Kzfr


kzlr = kzfr

(3.67) & (3.68)





kzlm = kzlp kzfr

Now rewrite equation 3.54 as,

a3 k 12
zlm U_

1. Set b = 0.5 inch.

2. a3A = constant, since b is set and kzlm can be

computed from previous information.


3. Select ab and find a.

4. Select XL and L and find X.

5. Compute E.

6. Compute a3A.

7. Compare the result of step (6) to the result

of step (2).

8. Repeat the process until the result of step (6)

equals the result of step (2).

9. Make sure the physical size of this layer is

compatible to the slatted layer, especially

in terms of slot spacing, i.e., blockage.

C. Storage Coefficient Manifolds. Refer to Figure 3.5

for this section (bm = 0.5").

1. Take average storage coefficient and average

depth in prototype to compute S

2. Select convenient time ratios, in this case

1 minute = 1 year, and compute Sor from

equation 3.59.

3. Compute S om

4. Equation 3.60 is now used to compute A for
various 1 's with z averaged over 1 .
m m m

In this case the model was apportioned into five

zones with one manifold per zone.


Because the Hele-Shaw analog is capable of modeling

complex geometries and boundary conditions, it is desirable

that it be as adaptable to as many different prototype

geometries and hydraulic parameters as possible. This

would facilitate model construction and investigations of

many different areas in the state of Florida where saltwater

intrusion is, or in the future might be, a problem. A

reduction in cost of investigation would also be achieved

if many of the parts were reusable.

A list of general specifications would then be as


1. The Hele-Shaw model should be housed in a frame
in which it can be easily installed and removed.

2. The front and back plates with interior model
parts should not be permanently sealed together.

3. The front and back plates should be as adaptable
as possible to different situations.

4. The model should have as few opaque parts as

5. The model should be mobile.


As shown in Figure 5.2, the frame is composed of

two assemblies; the cradle and the cradle dolly. The

function of the cradle is to support and orient the

Plexiglas plates. It also contains the inflatable neoprene

hose which seals the plates. It is fabricated of 2-1/2"

x 2" x 3/8" steel angles which are welded into a channel