WATER POLLUTION CONTROL RESEARCH SERIES
16060 ELJ 03/72
DENSITY INDUCED MIXING
IN COUCINED AQUIFERS
U.S. ENVIRONMENTAL PROTECTION AGENCY
-------
WATER POLLUTION CONTROL RESEARCH SERIES
The Water Pollution Control Research Series describes the
results and progress in the control and abatement of pollution
in our Nation's waters. They provide a central source of
information on the research, development, and demonstration
activities in the water research program of the Environmental
Protection Agency, through in-house research and grants and
contracts with Federal, state, and local agencies, research
institutions, and industrial organizations.
Inquiries pertaining to Water Pollution Control Research Reports
should be directed to the Chief, Publications Branch (Water),
Research Information Division, R&M, Environmental Protection
Agency, Washington, D. C. 20460
-------
DENSITY INDUCED MIXING IN CONFINED AQUIFERS
by
L. W. Gelhar
J. L. Wilson
J. S. Miller
J. M. Hamrick
Ralph M. Parsons Laboratory for Water Resources and Hydrodynamics
Department of Civil Engineering, Massachusetts Institute of Technology
Cambridge, Massachusetts 02139
for the
Office of Research and Monitoring
ENVIRONMENTAL PROTECTION AGENCY
Project #16060 ELJ
March 1972
-------
EPA Review Notice
This report has been reviewed by the Environmental
Protection Agency and approved for publication.
Approval does not signify that the contents neces-
sarily reflect the views and policies of the Environ-
mental Protection Agency, nor does mention of trade
names or commercial products constitute endorsement
or recommendation for use.
11
-------
ABSTRACT
Analytical techniques are developed to describe the mixing between two
fluids of different density in a confined aquifer, in which one fluid
is introduced to the aquifer by well recharge. The immiscible displace-
ment process in both linear and radial flows is analyzed and the effects
of longitudinal and lateral dispersion are included using a boundary
layer approximation. The theoretical results demonstrate the effect of
hydrodynamic dispersion in retarding gravity segregation due to density
differences.
The theoretical results are compared with observations of aquifer mixing
in linear and radial flow laboratory models. During recharge excellent
agreement between the theoretical predictions and experimental results
is found and the predicted retarding effects of longitudinal dispersion
are verified. During withdrawal some systematic differences between the
theory and observation are noted. Theoretical predictions of recovery
efficiency during a recharge-storage-withdrawal sequence show trends
similar to those observed but are typically 5 to 10% lower than those
observed.
Direct theoretical predictions of recovery efficiency during single or
multiple sequences of recharge-storage-withdrawal are developed for an
immiscible system, and similar developments outlined for miscible dis-
placement . The results are suggested for use in developing preliminary
estimates of operating conditions in the field.
This report was submitted in fulfillment of Project Number 16060ELJ
under the sponsorship of the Office of Research and Monitoring, Environ-
mental Protection Agency.
iii
-------
CONTENTS
Section Page
Conclusions 1
Recommendations 3
1 Introduction 5
2 Analysis of Density Induced Mixing in Linear Flow 9
3 Analysis of Density Induced Mixing in Radial Flow 31
4 Linear Flow Experiments 71
5 Radial Flow Experiments 79
Acknoxjledgments 111
References 113
Definition of Symbols 117
Appendices 121
-------
FIGURES
Figure Title Page
2.1 Definition Sketch for Linear Flow 10
2.2 Comparison of One-Dimensional Concentration Model 10
2.3 Theoretical Horizontal Interface Projection for
Linear Flow 23
2.4 6/L as a Function of Time for the Case of No
Dynamic Effects of Grain Dispersion 25
2.5 6/L as a Function of Time for Interaction between
Grain Dispersion and Density Induced Mixing 26
2.6 g(L/6) for One-Dimensional Linear Flow 29
2.7 Total Dispersion Coefficient, Et, as a Function
of Time for Linear Flow 29
3.1 Definition Sketch of a Radial Confined Aquifer 32
3.2 Comparison of Interface Shapes for Various Values
of e During a Recharge Operation 37
3.3 Inverse Interface Slope, 0, as a Function of Time
During Recharge " 56
3.4 Inverse Rate of Change of Relative Concentration,
At, as a Function of Time During Recharge 57
3.5 A/BH as a Function of Time During Recharge for the
Case of no Dynamic Effects of Grain Dispersion 59
3.6 A/3H as a Function of Time During Recharge for the
Case of Interaction betx^een Grain Dispersion and
Density Induced Mixing 62
3.7 Sketch of the Interface at the End of Withdrawal 66
3.8 Recovery Ratio, R, as a Function of e and e
(ts=0) r w ^
4.1 View of the Linear Sand Box Model 72
4.2 Horizontal Projection of Interface, L, as a
Function of Time 75
vii
-------
Figure Title Page
4.3 A Typical Linear Flow Experimental Breakthrough
Curve Compared with Theory 77
5.1 View of the Radial Model 80
5.2 Photographs of the Radial Model in Operation 81
5.3 Location of Conductivity Probes 83
5.4 Detail of Well Construction 85
5.5 Calibration of Averaging Conductivity Probes 88
5.6 Media Dispersivity Calculated from the Averaging
Probes 93
5.7 Piezometric Head Distribution in Radial Model 94
5.8 Comparison of Visual and Electronic Measurement
of the Interface with the Immiscible Theory 97
5.9 Comparison of Theory and Experimental Time
Dependent Concentration Slopes During Recharge 99
5.10 Volumetric Length of the Interface, gH = AV, as a
Function of Time 100
5.11 Comparison of the Coupled Theory and Experimental
Time Dependent Concentration Slopes During
Recharge 101
5.12 Experimental Breakthrough Curves at the Well 103
5.13 Experimental Breakthrough Curves During Recharge
and Withdrawal 104
5.14 Comparison of Experimental and Theoretical
Breakthrough Curves During Recharge 105
5.15 Comparison of Experimental and Theoretical
Breakthrough Curves During Recharge 105
5.16 Comparison of Experimental and Theoretical
Breakthrough Curves During Withdrawal 105
5.17 At as a Function of Time During a Recharge-Storage-
Withdrawal Sequence, Run A-2 106
5.18 At as a Function of Time During a Recharge-Storage-
Withdrawal Sequence, Run A-7 106
viii
-------
TABLES
Table Title Page
3.1 Summary of Equations Applicable to a Recharge,
Storage, Withdrawal Operation 69
4.1 Range of Test Results for the Linear Model 76
5.1 Radial Model Horizontal Hydraulic Conductivity 89
5.2 Determination of Dispersivity from Averaging
Probes 92
5.3 Experimental Conditions for the Radial Model 96
5.4 Observed Values of -QAt/V during Withdrawal 108
5.5 Comparison of Theoretical Recovery Ratios with
Experiments 109
IX
-------
CONCLUSIONS
1. Analytical techniques have been developed to describe the combined
effects of density difference and hydrodynamic dispersion on mixing in
linear and radial flows in confined aquifer systems. The theoretical
results demonstrate the retarding effects that dispersion may have on
the rate of tilting of the interface. The roles of both lateral and
longitudinal dispersion in the mixing are established.
2. Laboratory experiments with linear and radial flow aquifer models
confirm several aspects of the theoretical developments. During re-
charge operations the predictions are found to be in excellent agree-
ment with the observed interface tilting. In most of the recharge
experiments, the interface motion is dominated by immiscible dynamics
but some experiments confirm the predicted retarding effects of hydro-
dynamic dispersion. However, during withdrawal there are significant
differences between the observed and predicted interface motion. Pre-
dictions of recovery ratios showed trends similar to those observed
although the predicted values were typically 5 to 10 percent lower.
3. Direct application of the analytical results in the evaluation of
field operations is feasible. The results of the immiscible theory
(Eq. 3.119) provide conservative estimates of the recovery ratio for a
recharge-storage-withdrawal sequence in terms of aquifer and fluid pro-
perties and flow rates. Multiple recharge-storage-withdrawal sequences
can be treated and improved recovery efficiency is found after a few
sequences. Conditions under which dispersive effects become important
can be evaluated (see Section 3.4) and these effects can be included
using the computational procedure developed in Section 3.5.
-------
RECOMMENDATIONS
1. Several additional factors of importance In field situations should
be considered in future analytical and experimental studies; these in-
clude partially penetrating wells, phreatic aquifers, sloping and multi-
layered aquifers, and ambient flow in the aquifer.
2. Analytical or numerical solution techniques, should be developed to
evaluate the implicit limitations of the proposed analytical models;
these should include evaluation of the effects of interface curvature
during the storage and withdrawal cycles and influence of no-flux boun-
dary conditions at the aquifer boundaries.
3. Field experiments are suggested for further evaluation of the pro-
posed analytical techniques. This involves the development of a series
of experiments in which breakthrough curves are observed at several
monitoring wells and at the recharge-pumping well during a recharge-
storage-withdrawal sequence. The hydraulic conductivity and dispersivity
of the aquifer should be determined as part of the experimental program.
-------
1. INTRODUCTION
1.1 General Discussion
The pressures of population and economic growth lead inevitably to
problems in water resource management. One method of trying to optimize
present water resources to keep up with growth involves the utilization
of ground water aquifers for storage or disposal of surface or waste
waters. During periods of the year when there is excess surface water
available, it can be pumped into aquifers through wells for storage and
use during dry periods. This use of underground aquifers for water sto-
rage is attractive because of the scarcity of natural reservoir sites,
the low initial costs compared to ordinary surface reservoir construc-
tion, the elimination of evaporative losses, and the possibility of par-
tial purification of the stored water during its passage through the
aquifer.
These storage schemes can be used to either mix and dilute contaminated
surface waters with naturally occuring ground water, or for storage of
surface water of high quality in an aquifer initially containing water
of lower quality. This latter case in which fresh water is stored, for
example, in a saline aquifer could find wide application in the United
States where one third of the land area has saline ground water less
than 500 feet below the land surface (see Kohout, 1970). Moulder (1970)
has reported field studies on the recharge of saline aquifers with fresh
water and has suggested specific geographic areas where this scheme may
be feasible.
Ground water aquifers are also finding wide use for disposal of liquid
wastes and the trend is toward increased use of aquifers for this pur-
pose. The number of industrial waste injection wells, for example,
increased from close to 100 in 1967 to nearly 150 in 1969 (see Kohout,
1970). In addition to the disposal of normal wastes, ground water
aquifers can also be used to dispose of thermally polluted waters.
When surface water is introduced into an aquifer, mixing occurs with the
native ground water. Usually the surface water and native ground water
differ in physical, chemical, and biological characteristics. In de-
signing and operating recharge or disposal systems, it becomes necessary
to predict and control, as much as possible, the amount of mixing between
the injected and native waters.
The mixing in the aquifer is caused by the large scale convective motion
of the native ground water as well as the flow from wells during recharge,
disposal, or recovery. Further mixing is due to hydrodynamic dispersion
which includes the effects of molecular diffusion as well as what may be
termed grain dispersion. Grain dispersion is caused by the velocity vari-
ation of the fluid flowing along different paths within the pore structure
of the medium. A final mixing process comes about because of differences
in the physical properties of the native and injected waters. Differences
in density may be a significant cause of mixing due to gravity currents
-------
that develop. The resulting convective motion increases mixing in the
aquifer. The storage of fresh surface water in saline aquifers is an
example where density induced mixing is important.
1.2 Review of Previous Work
Considerable work has been done on hydrodynamic dispersion and its
effects on the mixing of miscible fluids flowing in porous media.
Schiedegger (1961) and Bear and Bachmat (1964) theorized a linear re-
lationship between the dispersion tensor and the velocity of flow
through porous media. Harleman and Rumer (1963) and List and Brooks
(1967) gave experimental evidence that some nonlinearity between the
velocity and dispersion coefficient may exist. Studies by Li and Yen
(1968) indicated that differences in viscosity and density have little
or no effect on the local mixing due to dispersion at the interface.
More recently Gelhar and Collins (1971) developed an analytical solu-
tion that describes the concentration of a substance at a particular
radius for flow from a well. Numerical solutions of the dispersion
equation have been developed by Shamir and Harleman (1967), Reddell and
Sunada (1971), and Guymon, Scott, and Hermann (1970). Rose and Passioura
(1971) showed that current methods for calculating the coefficient of
hydrodynamic dispersion during lab tests may be profoundly affected by
even seemingly negligible density differences between the original fluid
and the fluid used as a tracer during tests. Furthermore Theis (1967)
and Mercado (1967), demonstrated that field determinations of dispersion
coefficients may be several orders of magnitude greater than those in-
dicated by laboratory tests because of the large scale heterogeneity of
natural deposits. These factors should be considered when studying and
applying results from tests on hydrodynamic dispersion in the laboratory.
There has also been extensive work done on salt water intrusion in
coastal aquifers. Dispersion and density induced mixing are both im-
portant in the study of salt water intrusion. Henry (1960) studied the
combined effects of gravity and dispersive mixing for a steady saline
wedge in a confined aquifer. Rumer and Harleman (1963) studied the ef-
fects of tidal motion on interfacial dispersion. Bear and Dagan (1964)
developed an approximate analysis of the unsteady movement of a salt
water wedge due to a change in the fresh water flow. They compared their
results with model experiments. Shamir and Dagan (1971) developed numer-
ical solutions for the motion of the seawater interface in coastal
aquifers.
There have been a few investigations of dispersion and density induced
mixing more directly related to this paper. Esmail and Kimbler (1967)
and Kumer and Kimbler (1970) studied the possibility of storing fresh
water in saline aquifers. Their work was derived from an "ad hoc" as-
sumption of extending results based on the rate of interfacial tilting
in a linear model with no net flow to radial flow from a well using a
numerical procedure which was not explicitly described in the papers.
Although their procedure involves this extrapolation from linear to
radial flow, they did report good agreement between their numerical
-------
simulations and data collected in a three dimensional radial flow model.
Perrine and Gay (1964) developed numerical solutions of the problem of
density induced motion during radial flow from a well with density and
viscosity difference between the two fluids. Their solutions did not
take into account hydrodynamic dispersion at the interface, and no com-
parison was made between their theory and experimental or field data.
Green and Cox (1968) developed numerical solutions of the equations of
motion and mass conservation based on the method of characteristics.
They analyzed miscible displacement of fluids with different density and
equal viscosity in a linear model using constant dispersion coefficients
in a numerical scheme. Results of the numerical analysis were compared
with those from a laboratory model. Keulegan (1954) and Gardner, Downie,
and Kendall (1962) reported experiments on the tilting of the interface
between fluids of different density in linear laboratory models with no
net flow. De Josselin de Jong (1960) developed the potential flow solu-
tion for the initial motion of a vertical interface in a linear model
with no discharge. Gardner, Downie, and Kendall (1962) have also used
the immiscible approximation to develop analytical descriptions of inter-
face tilting for small and large times under conditions of no net flow.
1.3 Objectives and Scope
The primary objective of this study is to investigate the effects of hydro-
dynamic dispersion and density stratification on mixing in ground water
systems. Of particular interest is the mixing process which occurs when
a fluid whose density differs from that of the native fluid is injected
into a confined aquifer. These mixing phenomena can be important in ap-
plications involving the storage of fresh water in saline aquifers or
deep-well waste disposal. Methods of evaluating and predicting these
mixing phenomena under field conditions are sought.
The general approach includes the development of analytical models which
describe the effects of density differences and dispersion for two flow
configurations. A series of laboratory experiments in also developed to
evaluate the analytical predictions.
The analytical methods are first developed for mixing in a linear flow
with a constant net convective velocity, and later extended to mixing in
a radial flow well system. Initially, the theoretical analyses are based
on immiscible displacement for the case of fluids with different density
and equal viscosity. The effect of dispersion on the density induced
tilting of the interface and the overall mixing is then analyzed theore-
tically using boundary layer type approximations. From this analysis it
is possible to determine when dispersion will influence the interface tilt-
ing and the overall mixing. The analyses also reveal the conditions under
which longitudinal or lateral dispersion will dominate the dispersive mix-
ing and influence the interface tilting. The relative importance of grain
dispersion and density induced mixing is characterized by a dimensionless
parameter for both flow regimes.
-------
Based on the results of the linear flow analysis., a concept of gross dis-
persion due to density difference is introduced and a total dispersion co-
efficient representing both density and grain dispersion is also derived.
The radial flow system is analyzed for the three processes of recharge,
storage and withdrawal. The results are used to simulate multiple se-
quences of recharge, storage and withdrawal.
The results of laboratory experiments designed to evaluate the theoretical
developments are presented. The linear flow case is simulated in a hori-
zontal linear sand model. In the radial flow case a 15° radial sector
model of consolidated plastic beads was constructed. Interfacial tilting
and gross mixing were measured by visual and electronic means.
The application of the theoretical results to field operations is consid-
ered and predictions of recovery efficiency in a recharge-storage-withdrawal
operation are developed.
-------
2. ANALYSIS OF DENSITY INDUCED MIXING IN LINEAR FLOW
The objective in this section is to develop methods of analysis which
can be used to describe the combined effects of dispersion and gravity
segregation on aquifer mixing. The flow configuration considered here
is a one-dimensional displacement in a confined system with horizontal
boundaries, i.e., a linear flow. A theoretical analysis of immiscible
displacement in a model with a net flow is developed for the case of
fluids with different density but equal viscosity. Based on this
analysis, a new concept of gross dispersion due to density differences
is introduced. Using boundary layer approximations, the effects of
grain dispersion on the tilting of the interface and the overall mixing
are analyzed. The effects of both longitudinal and lateral dispersion
are included in the analysis. In Section 3, these same concepts and
methods of analysis are extended to the case of the radial flow system.
2.1 Immiscible Displacement in a Confined Aquifer of Constant Thickness
The problem of mixing of two liquids of different density in a homogen-
eous, confined aquifer may be simplified by neglecting grain dispersion
and molecular diffusion. In such a case, the mixing process is des-
cribed by an immiscible displacement model. Consider the displacement
of a heavier liquid with density p.,, by a lighter liquid with density
p,, in a horizontal confined aquifer of constant thickness H (see
Figure 2.1). If the two liquids have approximately the same kinematic
viscosity, v, (i.e., fresh and salt water), the hydraulic conductivity,
K, may be taken as a constant throughout the flow.
The following analysis is based on the Dupuit approximation, which
assumes hydrostatic pressure distributions and thus essentially hori-
zontal flow. The analysis is strictly correct only when the interface
slope is small and hence the vertical velocities are small.
The continuity equations for the two liquids are:
IF (2-2>
where £ is the distance to the interface, in liquid 1, measured from the
top of the aquifer, and r\ is the distance to the interface, in liquid 2,
measured from the bottom of the aquifer. The seepage velocities of the
two liquids, u.. and u«, are given by the Darcy equations for the two
liquids,
K 3*1
Ul = - --jr-i (2.3)
1 n 3x
-------
1
«f»1-H/2
u
<}>0-H/2
H
//ANN //A\\//A.\s //A\\//A\V/A \V/ANV /A\\/ /' NV/ /\\ V//Qv
f-i L *]
Figure 2.1 Definition Sketch for Linear Flow
1.0
C/C .6 -
o
.4 -
.2 -
1 T
Eq. 2.33
Eq. 2.24
I
r L
-.6 -.4 -.2 _0 .2 .4 .6
x/L
Figure 2.2 Comparison of One-Dimensional Concentration Models
10
-------
K
Vertical velocities are neglected due to the Dupuit approximation. The
porosity of the aquifer is n, and the piezometric heads of the two liq-
uids are (j>. and~. Along the interface, the fluid pressure must be
continuous. The pressure continuity equation is
*2 = *l - (C-H/2) (2.5)
Ap = P2 - PI
If the density difference, Ap , is small compared to p9, and the rela-
tive density difference is identified by Ap/p_ = Ap/p, Eq . (2.5) becomes
(2-6)
The continuity equation for the entire flow is obtained by adding Eqs,
(2.1) and (2.2) and integrating
U;LC + u2n = UH (2.7)
The convection of the interface through the aquifer is due to the net
seepage velocity, U.
The solution of the problem consists of determining the interface shape
and location, and determining the velocities in the two liquids. Sub-
stituting Eqs. (2.3) and (2.4) into Eq. (2.7) gives,
(2'8)
Differentiating Eq. (2.6) and substituting into Eq. (2.8) and simpli-
fying the results using C + n = H, gives
^ ,i
1 nU _ Ap_ _n !n_ () Q\
8x ~ K p H 3x ^ '
Substituting Eq. (2.3) into Eq. (2.1) and eliminating £ using £ = H - n
gives
11
-------
The partial differential equation for n is obtained by substituting
Eq. (2.9) into Eq. (2.10).
.
.
3t 8x n p 2 2 ~ 3H
Introducing a conyected coordinate x - x - Ut , Eq. (2.11) is simplified
tQ
_
3t n p ~* 2 3H
Equation (2.12) may be reduced to an ordinary differential equation
using a similarity transformation.
= F(z) (2.13)
H/F
The resulting ordinary differential equation in z, the similarity
variable, is
2 23
2 dz ~ , 2 2 ~ 3
dz
Equation (2.15) is solved for an initially vertical interface by
assuming a power series in z, to give
n = ( H + ) (2.16)
This solution may not apply for small times, when the interface slope
is large and vertical velocities are not negligible.
The expression for the velocity in liquid 1 is obtained using Eqs .
(2.3) and (2.9).
u +
1 n p H
12
-------
Using Eqs. (2.4), (2.6) and (2.9) the expression for the velocity in
liquid 2 is obtained.
- IT - ^- (1 - Ih 8n (7
U2 ~ U n p (1 H} "3x~ (
From Eq. (2.16), the horizontal projection of the interface, L, is
L = 2H/f & I = 2H/F (2.19)
The position of the interface is given by x., the distance to the inter-
face from the x origin.
x£ = Ut +|y (2.20)
Results of Eq. (2.16) indicate that the interface is a straight line,
of decreasing slope, which is convected through the aquifer with a
velocity U. The solution is also valid for the case of no through
flow, U = 0.
If the density difference is due to the presence of some material in
solution (i.e., salt), and the density is a linear function of the
concentration of the material, c, then
p = P + Ap -- (2.21)
The concentration of material in liquid 2 is co. To obtain a one-
dimensional description of the variation of the liquid density, the
concentration, c, is averaged over the aquifer thickness.
H
I = M 2 c dy (2.22)
H J H
2
Using the indicated concentration distribution,
c = co * - I < y < - (f -
c = 0 , - (f - n) < y < f
(2.23)
13
-------
The average concentration is
n C T
c = c = ( x + 1 O ")L\
* IT T ^ o * \ *-*"/
which indicates a linear distribution of average liquid density,
p" = p + Ap |- (2.25)
o
The average flux of_dissolved material with respect to the convected
coordinate system (x,y) is
H
J = pcu dy (2.26)
H J g
" 2
"u = u - U
which is evaluated, using Eqs. (2.18) and (2.23), as
j^.coKAP./.TK.nlH
J co P2 n p U H; H 3x
Using Eq. (2.24), Eq. (2.27) is modified to the form
K Ap . TK 8c .
J " * P2 nF (1 "I)n ^ (2
For a dispersive mixing process, the flux is related to the concentra
tion gradient through a dispersion coefficient, E, as follows
(2.29)
The dispersion coefficient can thus be evaluated by equating Eqs.
(2.28) and (2.29) as
for P2 = p". The result of Eq. (2.30) suggests that the gross character-
istics of density induced mixing may be approximated by the one-dimen-
sional dispersion equation of the form
14
-------
H*'H-
(2.31)
where EI is a constant dispersion coefficient. Although Eq. (2.30)
implies a parabolic distribution of the dispersion coefficient, E,
over the length of the interfacial zone, satisfactory results can be
obtained by using the average value.
_ - _ H K Ap
h, h
1 6 n p
(2.32)
This dispersion coefficient, which characterizes the mixing due to
density differences* will be identified as the density dispersion co-
efficient. This expression for E is later used to estimate the im-
portance of the effect of the density difference on the mixing process.
The longitudinal distribution of average concentration can thus be
determined approximately by solving Eq. (2.31) subject to the initial
condition
t = 0
c = c
x>0
c = 0
x<0
which corresponds to an initially vertical interface at the origin
x = 0. The solution is
c
C
1,1 c,
- + - erf (
2 2
(2.33)
A comparison of this one-dimensional approximation with the exact re-
sult from Eq. (2.24) is shown in Figure 2.2.
In many field situations, it is the average water quality, over the
aquifer thickness, which is of importance. If the mixing process is
assumed to be due only to the density difference between the mixing
liquids, the one-dimensional model represented by Eq. (2.33) will ap-
proximately describe the mixing process, and predict the depth-averaged
water quality. The concept of a one-dimensional dispersion model to
represent a two-dimensional density induced mixing process is an im-
portant simplification of this phenomenon.
2.2 Grain Dispersion Effects
The immiscible displacement model is an idealization, since there is
always some degree of interfacial dispersion. The objective here is to
estimate the importance of grain dispersion in comparison to the density
15
-------
induced mixing.
The validity of the assumption of immiscible displacement depends on
the magnitude of the density difference Ap . In certain cases the
density difference will be large enough to completely govern the
dynamics of the mixing, (i.e., the velocity field is that predicted by
the immiscible model); however, the resulting model may not adequately
describe the variation of density near the interface. To obtain a
better description of the density variation, grain dispersion must be
considered.
The mass conservation equation which describes the distribution of
solute concentration c in a porous medium is [Bear, et al, 1968]
9c , 3c , 3c 3 ,_ 3cx , 3 ,_ 3cx / ~,s
^+U 3^+V^ = H °W + ^T (D2 37} (2'34)
where u and v are the horizontal and vertical seepage velocities respec-
tively, and D-. and D~ are the longitudinal and lateral dispersion coef-
ficients respectively. To simplify the solution of Eq. (2.34), the
following assumptions are made.
(i) The dispersion coefficients, D.. and D~ , are generally velocity
dependent, but for this analysis the velocxty changes are assumed small
and D and D are taken as constants and proportional to the net seepage
velocity U.
(ii) The vertical convective transport may be neglected when the Dupuit
approximation, which assumes negligible vertical velocities, is
applicable.
(iii) The convective velocity in the dispersed region can be approxi-
mated by the velocity of the interface as determined by Eq. (2.20).
3x.
u = u. = =- (2.35)
1 O L
The interface is idealized as the line of 50% concentration.
The three assumptions simplify Eq. (2.34) to give,
2 2
3c , 3c 9 c 3 c . ,,
3F + ui 3i " Di 7T + Dz 72 (2'36)
3x 3y
and the interface location is as given by Eq. (2.20). The interface
slope, S, and its reciprocal, 3, are defined as
(2-37)
16
-------
Equation (2.20) is rewritten as
x. = Ut + 3y
(2.38)
The following coordinate transformation will be used:
Xl =
(2.39)
Expressing the concentration c as a function of x , y, t, Eq. (2.36)
becomes
3c
3t
xi>:
^
- r> c + n 8 c
~ Dl , 2 + D2I 2
3x
Since the lines x = constant are parallel xvith the interface, it is
reasonable to assume that
= 0
which can be verified from the results (Eq. 2.44). The convective
dispersion equation then reduces to
. 0
o
c = 0 x<0
(2.43)
17
-------
The solution of Eq. (2.42) subject to these conditions is
c
c
(2.44)
If the dynamics of the mixing are given by the immiscible displacement
model, 3 is given by Eq. (2.19), and hence £ may be determined. From
Eqs. (2.19) and (2.41),
= 2/F
5 = D^ + 2D2T t
(2.45)
j. f
= 2 +2 erf
o
x - x.
(2.46)
The result of Eq. (2.46) indicates that as time increases, the impor-
tance of lateral dispersion increases. The approximate width of the
dispersed zone about the interface, 6, is defined as
6 =
3(c/co)
x=x.
1-
-1
)t
(2.47)
As the width of the dispersed zone, 6, becomes large, dynamic effects
of grain dispersion are expected, and the dynamics of the mixing are
not given by the immiscible model. The ratio 6/L, where L is given by
Eq. (2.19) is proposed as a criteria to estimate the importance of
grain dispersion.
1/2
7TD
(2.48)
The result given by Eq. (2.46) is expected to be valid if 6«L. As 6
approaches L it becomes necessary to consider the dynamic effects of
the grain dispersion. This analysis will be developed in the next
section.
2.3 Interaction Between Grain Dispersion and Density Induced Mixing
In Section 2.2 the additional mixing due to grain dispersion was deter-
mined under the assumption that grain dispersion has a negligible in-
fluence upon the dynamics of the process. As mixing due to grain
dispersion becomes large, the immiscible model of Section 2.1 is no
18
-------
longer adequate to describe the motion. To determine the influence of
grain dispersion on the fluid motion, the mass and momentum equations
must be solved simultaneously. Coupling the mass and momentum equations
for the entire flow field results in non-linear partial differential
equations, which are extremely difficult to solve. However, a useful
approximation can be developed by coupling these equations at the mean
interfacial position, x = x..
The governing equations are dynamic equations (Darcy equation)
u=-^f (2.49)
np 3x
(2'50)
where p is the pressure, k is the permeability, and y is the dynamic
viscosity. The equation for conservation of total mass is, under the
conditions of small density difference (Ap/p «1)
|^4^=0 (2.51)
3x 3y '
Using the same assumptions made in Section 2.2, the mass conservation
equation for a solute of concentration c is
2 2
3c 3c 3 c 8 c ,0 ,..
-3T + ui 3x~= Di7T + D27T (2-52)
3x 3y
The Dupuit approximation assumes a hydrostatic pressure distribution,
8p/3y = - pg, and thus from Eq. (2.49) it follows that
3u _ kg_ 3p _ K 1 3p
"^ ~ "^ "^ \^-
3y ny 3x n p 3x
where K is the hydraulic conductivity.
Equation (2.53) can be evaluated at the mean interfacial position,
x = x
3y
x=x.
i
(2.54)
If the width of the dispersed zone, 6 = AirC , is large in comparison
to the horizontal interface projection, L, the velocity gradient in
Eq. (2.54) can be approximated by
19
-------
_3_u
9y
3u.
x
(2.55)
x=x.
The range of validity of this approximation is evaluated in detail in
Appendix A. Using this approximation, Eq. (2.53) can be rewritten as
3u. 3 x.
K
x=x.
i
3y 3y3t n p 3x
Using Eq. (2.38), Eq. (2.56) may be rewritten in terms of 3
dg = K 1 3p
dt n p 3x
(2.56)
(2.57)
x=x.
The solution of Eq. (2.52) for the appropriate boundary conditions is
that obtained in Section 2.2.
c_
c
o
11 cf
=2+2 6rf
\
x-x,
J
(2.58)
dt
(2.59)
Differentiating Eqs. (2.21) and (2.58) and combining gives
_ Ap (X-XJ)2
(2.60)
To couple the mass and momentum equations at the mean interfacial posi-
tion, Eq. (2.60) is evaluated at x = x., and substituted into Eq. (2.57)
_d_B _ K A
dt = n P
(2.61)
Equation (2.61) must be solved simultaneously with Eq. (2.59) in dif-
ferential form.
(2.62)
The solutions of Eqs. (2.61) and (2.62) with the initial conditions
20
-------
3 = 5 = 0 at t=0
are, for the reciprocal slope 3,
»2 + i££-!l-* (2-63)
and for the dispersed zone width, 6,
The ratio of the dispersed zone width to the horizontal projection of
the mean interface is, from Eqs. (2.63) and (2.64),
Equations (2.63) and (2.65) may be approximated for both large and
small interface slopes, noting that D.. is usually an order of magnitude
larger than D . For large slopes, the lateral dispersion term is neg-
ligible and Eqs. (2.63) and (2.65) reduce to
:
If the interface slope is small (large 3), longitudinal dispersion is
negligible and Eqs. (2.63) and (2.65) become
(2.69)
It should be noted that, in the analyses of grain dispersion in Sections
2.2 and 2.3, the boundary conditions of no mass flux at the upper and
lower aquifer limits were not explicitly applied. The resulting solu-
tions are therefore strictly valid only for aquifers which are of infinite
vertical extent. It can be seen by direct evaluation of the vertical
21
-------
mass flux at the boundaries (y = ±H/2) that the boundary conditions are
not exactly satisfied. In Appendix B, some possible effects of this
boundary condition are explored. However, no explicit evaluation of the
effects was obtained, other than an indication that the boundary effects
may become important when 6/L is large.
2.4 Discussion of Theoretical Results
In the previous sections, theoretical models were presented to describe,
for a one-dimensional linear flow in a confined aquifer, the mixing of
two fluids due to the effects of density difference and grain dispersion.
A_lso introduced was the concept of the density dispersion coefficient,
E, the magnitude of which may be used to estimate the relative importance
of density difference in the mixing process.
In the first model (Section 2.1) only the effect of density difference
is included. An immiscible displacement process is assumed and results
describing the tilting of an interface as it is being convected through
the aquifer at the net seepage velocity are obtained. One character-
istic of this model is the horizontal projection of the interface, L,
given by Eq. (2.19). Previous investigators have obtained expressions
for L, for immiscible displacement in enclosed systems with no net flow^.
For the case of no net flow, Gardner, et al., (1962) obtained for small'
and large time respectively
L/H = /3 T
(2.70)
L/H = 2/F
(2.71)
The second result (Eq. 2.71) was originally developed by Dietz (1953)
and also agrees with that found in this study (Eq. 2.19). Gardner,
et al., (1962) also suggested an interpolation formula which encom-
passes Eqs. (2.70) and (2.71) in the form
12t
4+3T
(2.72)
In the development leading to Eq. (2,19) it was noted that the result
was independent of the net seepage velocity and thus the agreement
between Eqs. (2.19) and (2.71) is expected. The results of the immis-
cible displacement theory (Eq. 2.19) and Gardner's interpolation form-
ula (Eq. 2.72) are shown in Figure 2.3.
In the second model (Section 2.2), the effects of grain dispersion were
considered under the assumption that grain dispersion has no influence
upon the dynamics of the flow. Under this assumption, the horizontal
projection of the interface, now idealized as the line of 50 per cent
concentration, is the same as the immiscible displacement result,
22
-------
E/D,
A Upper^ Limit of Immiscible Displacement Model
for E/D = 87f/3 (5/1X1/2)
B Lowejr Limit of the Coupled Analysis
for E/D = 8ir/3 (6/L=2)
10.0
Immiscible Displacement
Eq. 2.19
Eq. 2.72
L/H
Coupled Theory (Eq. 2,63)
2
1.0
6/L <2
Increasing Effect of
Grain Dispersion
0.1
0.1
1.0 10.0 100.0
T = KApt/npH
Figure 2.3 Linear Theoretical Horizontal Interface Projection
1000.0
-------
Eq. (2.19). An estimate of the width of the dispersed zone, 6, about
the interface is given by Eq. (2.47), As indicated in Section 2.2, the
immiscible model (Eq. 2.19) is valid when 6/L is small compared _tp unity.
The ratio 6/L, as evaluated from Eq. (2.48) for three values of ₯/D,
is shown in Figure 2.4 as a function of the dimensionless time, T. From
Figure 2.4, it is seen that even for large values of E/D , mixing due to
grain dispersion will eventually become as important and influence the
rate of interface tilting. The exact limits of applicability of the
immiscible model can not be specified but reasonable estimates are
obtained by assuming that the limiting value of 5/L is a fixed number
less than one, say 6/L - 1/2. For a given value of E/D , the indicated
maximum value of dimensionless time, T, for which the immiscible solu-
tion applies, is found from Figure 2,4 or Eq. (2.48). For example, when
E/D.. = 8rr/3 the maximum T = 15, as shown in Figure 2.3.
In the third model (Section 2.3), effects of density difference and
grain dispersion are analyzed by an approximate coupling of the mass
and momentum equations based on large 6/L. The solution is given by
Eq. (2.63), which represents a family of curves relating g, the inverse
interface slope (g^L/H), and T for various values of E/D and D /D .
Equation (2.63) has two asymptotes, Eqs. (2.66) and (2.68), for large
and small interfacial slopes or small and large times ,_respectively.
In Figure 2.3, Eq. (2.63) is shown for four values of E/D.. , with D /D_
constant and equal to 10, As noted in Section 2.3, the coupled analysis
(Eq. 2.63) is valid for 6/L»l. The ratio 6/L from Eq. (2.65) is shown
in Figure 2.5 for four values of E/D . An approximate limit on the range
of application of Eq. (2,63) can be^ found by taking 5/L = 2 in Eq. (2.65)
or Figure 2.5. For example, when E/D. = 877/3, the limiting value is
T = 200 and the coupled analysis (Eq. 2.63) should be valid for times
greater thin this value as indicated in Figure 2.3. It can be seen that
the regions of validity of the immiscible theory (Eq. 2.19) and the
£pupled analysis (Eq. 2.63) are mutually exclusive for a given value of
E/D . In general, there will be an intermediate time range in which
neither theory is applicable. In the example E/D.. = 8ir/3, the range is
15 < T < 200. The actual behavior in this range is not known, but a
smooth transition between the models is expected. The portions of the
curves representing Eq. (2.63) for 5/L < 2 are shown as_dashed lines in
Figure 2.3. It is seen that, even for large values of E/D , lateral grain
dispersion will eventually influence the interface tilting. From Figure
2.3, it is seen that as E/D decreases, interface tilting becomes less
and as lateral dispersion becomes important, the rate of tilting decreases.
The ratio 6/L indicates the relative importance of the two mixing mech-
anisms, grain dispersion and density difference. When 6/L is small
(Section 2.2), the grain dispersion has no effect on the interface
tilting and mixing dynamics are governed by the immiscible solution.
When 6/L is large the combined effects of the two mixing mechanisms are
important and the coupled analysis governs the mixing. The coupled
solution does not reduce to the immiscible solution as D ,D« -> 0,
because the regions of validity for the two analyses do not overlap.
24
-------
Ul
10.0
6/L
1.0 -
0.1
10.0
= K _Ap_ t_
T n p H
100.0
1000.0
Figure 2.4 6/L as a Function of Time for the Case of No Dynamic Effects of Grain Dispersion
-------
10.0
6/L
0.1
0.1
1000.0
Figure 2.5 6/L as a Function of Time for Interaction Between Grain Dispersion and
Density Induced Mixing
-------
A numerical example for some typical field conditions will illustrate
the importance of the different mixing phenomena. Consider a sandy
aquifer 50 feet thick having a hydraulic conductivity, K, of 0.03 ft/sec
and a porosity, n, of 0.35. For a_one percent density difference, the
density dispersion coefficient is E = 7.15 x 10"^ ft /sec. For a hy-
draulic gradient of 0.01, the net seepage velocity from Darcy's law is
U = 8.6 x 10~^ ft/sec. The longitudinal dispersion coefficient is esti-
mated to be DI = 2.8 x 10~° ft/sec from results by Harleman, Melhorn
and Rumer (19o3) and D.. /D? = 10 is assumed. The ratio E/D is approxi-
mately 26, indicating that the mixing process is initially one of immis-
cible displacement. The immiscible model (Eq. 2.19) will be valid up to
a dimensionless time T = 57 which is determined from Eq. (2.48). Thus a
period of somewhat over one month would be required before the interface
tilting would be influenced by grain dispersion. After one month, the
horizontal projection of the interface, L, is 670 feet and the width of
the dispersion zone about the interface, as given by Eq. (2.48), is
300 feet.
The results of the analyses of the combined effects of grain dispersion
and density- induced motion can also be interpreted in terms of a total
dispersion coefficient, E , defined by
(2'73)
where j is the mass flux relative to the mean flow U, and c is the
depth-averaged concentration given by
fH/2
- 1
C = H
c dy (2.74)
-H/2
The total dispersion coefficient, E , can be evaluated directly from
the results of Sections 2.2 or 2,3 in the general form
.. 1 x-x.
^"2 + 2«f~ <2'75>
where x. = Ut + 3y (Eq. 2.38), 3 is given by Eqs. (2.19) or (2.63) and
C by Eqs. (2.45) or (2.64), The average concentration is determined
using Eq. (2.75) in Eq. (2.74) with the result
2 ?
-Q
[0erf9-6erf9+(e -e i)/v^T] (2.76)
- , n-=- -6 -Q
2211
where 0 = (x+3H/2)//4T , 9 = (x-3H/2)//4T and x = x - Ut. The gradient
of the depth-averaged concentration is then
27
-------
23H
(2.77)
The mass flux, j , is expressed as
H
H/2
PD
3c
rH/2
-H/2
pcu dy
(2.78)
-H/2
where u = u - U. Assuming small density changes, Eq. (2.78) is adequately
approximated by
rH/2
cu dy
(2.79)
-H/2
Using Eqs. (2.73) and (2.79), the total dispersion coefficient becomes
rH/2
Et -
cu dy)/ -r
(2.80)
-H/2
The convective flux integral in Eq. (2,80) was evaluated at x=Ut (x=0)
in Appendix B (Eq. B.6) as
rH/2
, u x jj
cu dy = OQU ir-
_!
H
-H/2
in which L = 3H,_5 = ATT? and $ is defined in Eq. (B.7).
Eq. (2.77) with x = 0 and Eq. (2.81) in Eq. (2.80),
VD1 3 F T
L. J. _ J ,£j . f"\
(2.81)
Introducing
(2.82)
where g(L/6) = $(L/6)/erf (/JrL/26). The right hand side of Eq. (2.82)
represents the relative increase of the total dispersion coefficient
over the longitudinal dispersion coefficient for the medium. The func-
tion g(L/<5) was evaluated numerically and is shown in Figure 2.6. For
a given value of the parameter E/D.. , L/6 is determined by the appropriate
model from Section 2.2 or 2.3 (Eqs. 2.48 or 2.65) as discussed earlier_
in this section. The results of this computation for three values of
are shown in Figure 2.7. The dashed portion of the curves represent an
estimated interpolation in the region where neither of the two theories
for <5/L is strictly valid. It is seen from Figure 2.7 that the total
28
-------
g(L/6)
o.i
100.0
Figure 2.6 g(L/6) for One Dimensional Linear Flow
100.or r-
I!
Di
i MI
E/D = 8^/3
Immiscible Thy
Transition zone
Coupled Thy
TT/6
Figure 2.7 Total Dispersion Coefficient, E , as a Function
of Time for Linear Flow
29
-------
dispersion coefficient, E , decreases with increasing time and that
eventually the additional dispersion due to density differences becomes
negligible.
Rose and Passioura (1971) reported experiments in which the additional
dispersion associated with density differences was determined from break-
through curves. These experiments were done with columns of circular
cross-section and therefore are not directly comparable with the present
results. However, the increase in dispersion coefficient indicated from
Eq. (2.82) for these experimental conditions is found to be of the same
order of magnitude as that observed.
30
-------
3. ANALYSIS OF DENSITY INDUCED MIXING IN RADIAL FLOW
In this section the methods of analysis developed for linear flow in
Section 2 are extended to the more practical case of radial flow from
a fully penetrating well in a horizontal confined aquifer. The approach
generally parallels that for linear flow including first an analysis of
immiscible displacement, next a superposition of dispersive effects on
that result and finally a consideration of the combined effects of den-
sity differences and grain dispersion. The analysis is generalized to
treat a recharge-storage-pumping sequence that would be typical of field
well operations.
3.1 Immiscible Displacement
As in Section 2.1, the problem of mixing two liquids of different density
in a homogeneous, confined aquifer may be simplified by neglecting grain
dispersion and molecular diffusion. In such a case, the mixing process
is described by an immiscible displacement model. Figure 3.1 illustrates
the displacement of a heavier fluid of density, p~, by a lighter fluid
of density, p7 , in a horizontal confined aquifer of constant thickness H.
If the two liquids have approximately the same kinematic viscosity, v,
the hydraulic conductivity, K, may be taken as a constant throughout the
flow.
The following analysis is based on the Dupuit approximation, which as-
sumes hydrostatic pressure distributions and thus essentially horizontal
flow. The analysis is strictly correct only when the interface slope is
small and hence the vertical velocities are small.
The continuity equations for the two fluids are
where t, is the vertical distance to the interface, in liquid 1, measured
from the top of the aquifer, and n is the vertical distance to the in-
terface, in liquid 2, measured from the bottom of the aquifer. The
seepage velocities of the two fluids are given by the Darcy equation,
31
-------
Fluid 2
/ \\y/L\\v/,\v//, \\v/.
Figure 3.1 Definition Sketch of a Radial Confined Aquifer
32
-------
Vertical velocities are neglected due to the Dupuit approximation.
The porosity of the aquifer is n, and the piezometric heads of the two
liquids are .. and -. The continuity expression for the entire flow
field is obtained by adding Eqs. (3,1) and integrating,
?un + nu_ = -r-* = H (3.4)
1 2 2irnr r
where A = Q/2TrnH. The convection of the interface through the aquifer
is due to the net seepage velocity, A/r.
Along the interface, the fluid pressure must be continuous (p1 = p?),
or
where
Ap = P2 - PI (3.6)
For small density differences Eq. (3.5) becomes
4>2 = 4>-, (C-H) (3.7)
in which the relative density difference is identified by Ap/p~= Ap/p.
The problem involves the determination of the interface shape and
location. Substituting Eqs. (3.2) and (3.3) into Eq. (3.4) gives
S^-i 9o n
-L^!^_ U f n n \
Differentiating Eq. (3.5) and substituting into Eq. (3.8) with H = n +
results in the expression
^ J%
1 _ An _ A_p n_ 9n. n QN
3r rK p H 8r ^ ;
Substituting Eq. (3.2) into Eq. (3.la) and eliminating t, using £ = H -
gives
33
-------
Replacing the derivative of with Eq. (3.9), the partial differential
equation for r) is found to be
where D = - ^- H.
n p
Recharge Cycle Because of the nature of the flow geometry, it is helpful
to express Eq. (3.11) in a volumetric coordinate system. During recharge,
with a positive flow rate Q into the aquifer, the transformations are,
₯ = (TrnH)r2 (3.12a)
^* = Qrt (3.12b)
which represent, respectively, transformed space and time coordinates.
Using these transformations, Eq. (3.11) becomes
3n , 9n _ 2D.-r) ,.. aw -§ai f-\ i -i\
o v*. J V A 11 n o~V~
*
Note that this substitution transforms the radially dependent coef-
ficient of the convective term into a constant. The equation can be
simplified further by transforming to a convective coordinate system:
₯ = ₯ - *a (3.14)
Eq. (3.13) now becomes
This recharge equation can be reduced to an ordinary differential equa-
tion using a similarity transformation
f = n/H = F(z)
(3.16)
z - Wer₯*
where
(3.17)
-------
Ar = Qr/27mH (3.18)
and the subscript "r" refers to the recharge cycle. The resulting
ordinary differential equation in z, the similarity variable, is
The parameter e indicates the relative importance of the density dif-
ference, represented by D, as compared to convection, represented by A.
For increasing e , the relative rate of tilting of the interface should
increase.
For small e , a perturbation solution for Eq. (3.19) of the form
f = fQ + E/-L + erf2+'" (3'20)
is assumed. For equations of 0(1) this results in a nonlinear equation
in f and z
o
This equation is solved for an initially vertical interface by assuming
a power series in z , to give
(3.22)
The equations for higher orders, e , are of the general linear form
A
dz
(3.23)
n = 1,2,3,'
For n = 1 the equation is
2 d dfl 1 2
Z) ± - 2-z~^ - 2f1 = - y(l-3zZ) (3.24)
dz
and the particular solution is
fl = I6(1"3z) (3<25)
35
-------
Similarly, the equation for n = 2 is
df
- 2z
' ' 2 dz 2 3_
dz
for which the particular solution is
1- ~ T^z2) (3.26)
Therefore, during recharge, the interface behaves as
2
£ 2
- (l - i Z2)+- (3.27)
The solutions of Eq. (3.23) will, in general, include two solutions of
the linear homogeneous equation and a particular solution which depends
on F. However, the homogeneous equation, which is closely related to
the Legendre equation, can be shown to have mathematical singularities
at z = ±1 and therefore the solutions of the homogeneous equation must
be excluded on physical grounds. The result represented by Eq. (3.27)
is shown in Figure 3.2 for two values of e . It can be seen that,
even for the relatively large value e = .65 , the departure from the
linear interface shape of the zero order solution is not large. Note
that the upper and lower limits of the interface (f = 1 and f = 0) ad-
vance as e increases, but that the arrival time of f = 1/2 has been
delayed relative to the average displacement.
The length of the interface in volumetric coordinates is expressed as
_. . 21 21
which can be determined from Eq. (3.27) by evaluating the end points of
the interface,
(3-30)
Thus from Eqs. (3.16), (3.28), (3.29) and (3.30)
VJ1 -
= 2E V1 - ^r + ' (3.31)
36
-------
f =
-1.0
+1.0
z = V/e ₯.
r
f = n/H
-1.0
0
z =
+1.0
Figure 3.2 Comparison of Interface Shapes for
Various Values of e During a Recharge Operation
37
-------
The interface location can be described by a coordinate -V-. , the volu-
metric distance from the origin. For the zero order solution (Eq. 3.27
with e =0) the interface position is given by
(3.32)
where
g = 2er
and y is the vertical distance from the bottom of the aquifer. The
radius to the interface is expressed as r. and the interfacial velocity,
u, is
3r. A
(3-33)
where
r± = y^/TmH = r*/l+g (3.34)
/2A^F (3.35)
Storage Cycle During storage of liquid 1 in the aquifer, there is no
net flow (Q=0) and the convective term of Eq. (3.11) becomes zero, thus
111 - n i m - -^1-2- r 9nl
3t ~ D r 8rU1 H;H r 3^J
This equation may be transformed to volumetric form using Eq. (3.12a)
|2. = 2D(2TrnH) |^ [§(1 - £) V |j] (3.36)
and may be further simplified by introducing the coordinate system
£ =₯-₯
o
where ₯ is the volume of liquid of density p.. which was injected
during recharge. Equation (3.36) then transforms to
= 2D(2,nH) - [d - )(^ W-] (3.37)
38
-------
When the length of the wedge is small compared to the distance from
the well, ₯ « ₯ , Eq. (3.37) is approximated by
|a - 2D<2»«* ^[iU - j}> ffr-l (3.38)
This equation may be reduced to an ordinary differential equation using
the similarity transformation
,
z = (3.39)
t+Cn
o 1
where C.. is an arbitrary constant added in order to fit the initial
interface shape at the start of the storage period. The resulting or-
'dinary differential equation is
dz dz dz
(3.40)
which has the solution (see Eqs. 3.21 and 3.22)
f = |(l+z') (3.41)
The result of Eq. (3.41) represents a linear interface in volumetric
coordinates and is implicitly limited to cases in which the volumetric
length of the interface is small compared to the initial volume ₯ .
The explicit features of this limitation will be discussed in sections
3.4 and 3.5, but it can be seen from Eqs. (3.39) and (3.41) that this
solution will become invalid for large times. It is also implied that
the initial condition for the storage cycle should correspond to a
linear interface at the end of the recharge cycle and thus that e be
small.
It would be desirable to obtain solutions during the storage cycle
which are not limited by the approximation in Eq. (3.38), i.e., small
e during recharge. Although the equation governing the storage cycle
(Eq. 3.37) is similar to that for the recharge cycle (Eq. 3.15), the
method of solution used for the recharge cycle does not apply because
a similarity form of Eq. (3.37) can not be obtained.
Withdrawal Cycle During withdrawal, fluid is pumped from the aquifer
at a rate 0 and the governing equation becomes, from Eq. (3.11),
39
-------
»«>
where A = Q /2imE. This equation may be transformed to a volumetric
convective coordinate system using Eq . (3.12a) and
o w
₯ =
where
t = time at the start of withdrawal
o
₯- = net injected volume at t
o J o
Equation (3.42) then reduces to
= - 2e2 ?- [f(l-f )($ + *) M] (3.43)
with e = /D/A . When the length of the interface is relatively small
compared to the distance from the origin to the front, V;,. » V, Eq . (3.43)
is approximated by
- 51 = e2 1- [f (1-f) % (3.44)
3(C2-₯* ) 3^ 3*
This equation has the solution
(3.45)
where C« is a constant to be determined from the initial conditions for
the withdrawal cycle. The length of the interface in volumetric coordi-
nates will be given by
(3.46)
and the location of the interface is expressed as
2e
j / /
(3.47)
40
-------
The solution for the withdrawal cycle, as given in Eq. (3.45), represents
a linear interface in volumetric coordinates. The range of applicability
of this solution will, as discussed in Sections 3.4 and 3.5, depend on
the initial condition from a storage or recharge cycle. It can be seen
from Eq. (3.46) that as ₯A approaches zero (i.e., the entire stored vol-
ume is withdrawn) , the ratio AV/V^ becomes large and the approximation
implied in Eq. (3.44) becomes invalid. A perturbation solution of the
exact equation (Eq. 3.43) was attempted using the method employed for
the recharge cycle. A similarity form of Eq. (3.43) can be found, but
in the withdrawal case, the initial conditions can not be satisfied by
the similarity form.
Matching Between Cycles The solutions for the storage and withdrawal
cycles (Eqs. 3.41 and 3.45) can be combined with the zero order solution
for recharge (Eq. 3.22) to represent general recharge-storage-pumping
sequences. The general approach in matching between the cycles will be
outlined here and explicit applications to aquifer operations will be
developed in Section 3.5. Because each of the three solutions repre-
sents a linear interface in volumetric coordinates, the matching between
the cycles can be expressed conveniently in terms of the inverse slope
of the interface in volumetric coordinates, 3, which is given by
(3-48)
The constants C.. and C_ appearing in the storage and withdrawal solutions
can be expressed in terms of the initial inverse slope, 3 , at the start
of the cycle. From Eq . (3.41) for the storage cycle
3 = 2/4TTnHDV t+CL/H
o 1
and if the time t is measured from the start of the cycle, the constant
C can be determined in terms of 3Q. The resulting expression for the
inverse slope during storage is therefore
16irnDV t .
3 = g (1 + - = )1/Z (3.49)
32H
o
Likewise for the withdrawal cycle the inverse slope is, using Eq. (3.46),
where 3 is the initial inverse slope when V. = V , the initial net
q ., * o
injected volume.
The zero order solution for recharge may also be generalized to include
41
-------
initial conditions corresponding to an interface of inverse slope g
and initial injected volume V . Noting that the zero order solution
during recharge is represented by Eq. (3.15) with V « V , and defining
VA=Q t + V with t measured from the start of this recharge cycle,
the solution can be written
f = f [1 +
(3.51)
The constant ,C_ ,is again evaluated in terms of 8 and ₯ to yield
3 o o
(3.52)
The inverse slope, B, determined from any of the three cases can provide
the initial value B for any subsequent cycle.
Using the above notation, the solutions in all three cycles (Eqs. (3.22),
(3.41) and (3.45)) can be represented in the general form
f = 2+ BH (3'53)
where 3 is given by Eq. (3.49), (3.50) or (3.52).
3.2 Effects of Grain Dispersion
The immiscible displacement models of section 3.1 are an idealization
of the actual situation, since there will always be some degree of in-
terfacial dispersion. The objective of this section is to estimate the
importance of grain dispersion, under the assumption that the dispersion
does not attenuate the rate of interfacial tilting. This corresponds
to the case in x^hich the density difference, Ap , is large enough to
completely govern the dynamics of mixing (i.e., the velocity field is
predicted by the immiscible model) , but in which the resulting model
does not adequately describe the variation of density near the inter-
face. To obtain a better description of the density variation, grain
dispersion must be considered. The analysis of the recharge cycle will
be developed in detail, and the results for the other cycles will be
obtained as a simple extension of this case.
Recharge Cycle The convective-dispersion equation, neglecting molecular
diffusion, is, during the recharge cycle,
42
-------
+'r-7?+7
-------
and 3 = 2e ^A/H, the concentration is expressed as a function of -V^,
V- and y.
9c . i 3c i 3c 1 r_ a c
Ar/F
+ »2~
,2
r 8 c
1 2
3y
V'
- 28
(3.59)
V'
Since the coordinates, ₯. = constant, are parallel with the interface,
it is reasonable to assume that
_§£
3y
= 0
The convective dispersion equation then reduces to
3c
,2, 9 c
C3.60)
When the dispersed zone around the interface is small compared with the
distance, r^, from the origin, the two convective terms on the left hand
side of Eq. (3.59) will cancel near the interface, r = r.. The second
term in brackets in the longitudinal dispersion term will also be neg-
ligible under these conditions because it involves a lower order of dif-
ferentiation than the first. This type of approximation can be formally
justified, as done by Gelhar and Collins (1970), using perturbation
methods and can be verified from the result of the present analysis. It
is also consistent to replace the coefficients (A/r) in the dispersive
terms by A/r.,.. The coefficient in this equation can be eliminated by
the' introduction of a new variable, £.
(3.61)
The function, £, is evaluated from the known flow field, based on the
immiscible displacement model, Eqs. (3.35) and (3.31).
2
al
[rri
3r
H
(3.62)
44
-------
With this new variable, Eq. (3.60) transforms to
(3-63>
The initial condition corresponding to a vertical interface at the
origin (V=0) is
t = 0 c = c ₯ > 0
o
c = 0 ₯ < 0
and the boundary conditions are
t>0 c=c ₯--><»
o
c = 0 ₯ < -°°
The solution of Eq. (3.63), subject to these conditions, is
c_ , 1 + 1 erf !L (3.64)
co 2 2 /45
The dispersed zone width, in volumetric coordinates A, may be defined as
-Tr e 2 r^
i- + -f-r-a2* (3.65)
K H
and the approximate radial width of the dispersed zone, 6, may be defined
as
fa., ir e 2 r^
= 2**\hr- + -^ -- ^ <3-66)
* V 3r* H2 5
As the width of the dispersed zone, 5, becomes large, dynamic effects
of grain dispersion are expected, and the dynamics of the mixing are
not given by the immiscible displacement model. The ratio, A/AV, where
AV is given by Eq. (3.31) and A is given by Eq. (3.65), is proposed as a
criterion to estimate the importance of grain dispersion in the radial
flow problem during recharge
2 2
A a IT 3e r .
(3.67)
45
-------
Equation (3.67) indicates that the lateral dispersion becomes important
for
_ 2 2
3c r* a
-2-->l (3.68)
The concentration distribution described by Eq. (3.64) is expected to
be accurate when A«A₯. However, as A approaches AV, it is necessary
to consider the interaction of grain dispersion and density differences
in the mixing dynamics.
Withdrawal Cycle Under the assumptions outlined at the beginning of this
section, the convective dispersion equation (Eq. 3.55) during withdrawal
for A = - A is
w
., A .2 A -2
9c . 9c w 8 c w 8 c ,
where the subscript refers to withdrawal. Using the same method of
analysis described for the recharge period, this equation can be trans-
formed to the volumetric coordinate system, V , measured from the inter
face.
-__ _ _ _ .
Vr 9V£ Vri 3V£ " Vr 3V " 3r
2 0 3V02
w X.
For withdrawal, VA is given by (see Section 3.1)
while V. and V. are defined as before, in terms of VA (Eqs. (3.57) and
(3.58)). Using the same approximations used to obtain Eq. (3.60) from
Eq. (3.59), Eq. (3.70) reduces to
2
If-'luS- (3'72)
at, d₯ /
where
(3-73)
46
-------
and C is an arbitrary constant to be determined by matching with the
concentration field at t=t . The appropriate solution of Eq. (3.72) is
of the form
where £ can be evaluated from the known value of 3 found for the with-
drawal period using immiscible displacement theory (Eq. (3.46) with
A₯=gH or Eq. (3.50)). Substituting this value of g into Eq . (3.73) yields
r "4 r TrnH 3/2 £W2 1/2 ,_ ** M . /o -7^
= - t ot ₯A + a V^ ' (C - )] + C (3.75)
If the form of C_ implied in Eq. (3.50) is used, and C^ is evaluated in
terms of the initial volumetric width of the dispersed zone, A = ATT? ,
when t=t and VA=₯ , the variable £ is given by
2
(irnH)
+ 1 (v 5/2_₯ )}] +T- (3.76)
j O ^ A-TT
As in recharge, the ratio A to AV represents the relative importance of
grain dispersion. The volumetric width of the dispersed zone A = /47r^
is given by Eqs. (3.75) or (3.76), and AV is obtained from Eqs. (3.46)
or (3.50). Thus the ratio A/AV can be expressed as
= r _ _ 1/2
l J
In this solution it has been assumed that the mixing dynamics are given
by the immiscible displacement theory. Therefore, A/AV must be small
for the solution to be valid.
Storage Cycle When the net convective motion in the aquifer is zero,
the additional mixing due to grain dispersion is negligible. Physically,
this is implied from the fact that the volumetric length of the inter-
face, AV, must be small compared with the mean volumetric distance to
the interface, V , for the solution during storage (Eq. 3.41) to be
valid. Since the width of the dispersed zone generally is proportional
to the square root of the distance traveled by the interface, and the
length of the interface is , under these implied limitations , small com-
pared with the distance from the origin, the additional dispersion during
47
-------
storage will generally be small. Thus it is consistent to assume that
the volumetric thickness of the dispersed zone around the interface does
not change during storage. For extremely large storage periods (e.g.
several years) it may, however, be necessary to include the broadening
of the dispersed zone as a result of molecular diffusion.
3.3 Interaction Between Grain Dispersion and Density Induced Mixing
In Section 3.2, the additional mixing due to grain dispersion was
analyzed, assuming that the dynamics of the mixing is described by the
immiscible displacement model (Section 3.1). As the mixing due to grain
dispersion becomes large, the immiscible displacement model no longer
describes the fluid motion. To determine the influence of the grain
dispersion on the motion, the mass and momentum equations must be solved
simultaneously. When the mass and momentum equations are coupled for
the entire flow field, a practically intractable system of nonlinear
partial differential equations is obtained. However, as in the linear
case, a useful approximation can be developed by coupling these equations
at the mean interfacial position, ₯=₯.,
The governing dynamic equations (Darcy equations) are
v = (3.78)
r ny 3r
in which v and v refer to radial and vertical seepage velocities re-
r y
spectively, p is the pressure, k is the permeability and p is the dynamic
viscosity. The governing mass conservation equation is
3c , 3c , 3c 1 3 /rv 3c,. , 9 ,_. 8cN /0 on^
K + Vr 3? + Vy 3? = 7 ^(Dlr 3?> + ^2 3^ (3'80)
which is simplified by the same assumptions appearing in Section 2.2.
3c . 3c IA| 32c , IA! 32C ,, _..
3F + Ui 3? = al r -2 + a2 r 7T (3'81)
9r 3y
The interface is assumed to be linear in the volumetric coordinate system
(Eq. 3.58) and 3 is an unknown time-dependent function which represents
the inverse volumetric interfacial slope.
The dynamic equations are manipulated in the same fashion as their linear
counterparts (Section 2.3). The Dupuit approximation implies a hydro-
static pressure distribution,
48
-------
9y
= -pg
and thus from Eq. (3.78) it follows that
r = kg_ _3p_ = K 1 _3p_
3y ny 3r n p 3r
(3.82)
where K is hydraulic conductivity.
Equation (3.82) can be applied at the interface
3v
3u.
i
r=r.
= 3y
when the thickness of the dispersed zone, A, is much larger than the
horizontal projected length of the wedge, AV. The analogous approximation
for the linear case is examined in Appendix A. The result of this ap-
proximation is
3u.
i
32r.
i
3y ~ 3y3t
3V.
3yvr. 3t
K _! ^p_
n p 3r
(3.83)
r=r.
Using the linearity between density and concentration changes (Eq. 2.21),this
becomes
3u. _, . 3c/c
i _ K A_p_, QS
3y n p 8r r=r.
(3.84)
Recharge Cycle The solution to Eq. (3.81) during recharge is obtained
by the approximate method used in Section 3.2, and is given by Eq. (3.64)
in which 5 is given by Eq. (3.61). However, in this case 3 is an unknown
function of time. Differentiating c/c with respect to r and evaluating
the result at the interface gives
3(c/cQ)
3r
2TtnH r.
i
r=r.
This result is used to evaluate the right hand side of Eq. (3.84), and
the left hand side of Eq. (3.84) is, after multiplying by r. , written
as
3r-
r. 3t3y
3r2
r. 3t"2r. 3t
i i
49
-------
where 3 = 9₯./3y has been introduced. It is consistent with the approxi-
mation of the convective dispersion equation (see Section 3.2, Eq. 3.60)
to approximate r. by rA in the last expression. Equation (3.84) then
reduces to
2
JIT ( V /o' ~ IZHT (3.85)
Equation (3.85) must be solved simultaneously with Eq. (3.61) in differ-
ential form,
(3-86)
In order to simplify the simultaneous solution of these two equations ,
a new variable, 9 = 3/₯,j.V2, is introduced. Equations (3.85) and (3.86)
then become
*-\
r - (3.87)
d*A
= [2a (7rnH)1/2 + - 2 r- 02]₯ (3.88)
* 2(7rnH)
which can be solved simultaneously to give
51/2 = [2ai(TrnH)1/2 6 + - "2 G3] (3.89)
1 6(TrnH)17/ 2e^
This result is -substituted into Eq . (3.87) in order to find the inverse
slope, 6 = VA 6, for the initial condition 3=0 at t=0,
A 4« 5/2
/ TTxl/2 .2 , a2 3 4 Er V* ., ...
a , OmH) 0 + - T7T w~ = T - 9 - (3.90)
The volumetric width of the dispersed zone, A = /4ir£, can be found from
Eq. (3.89)
2a a0r.
A. [-1 + ^-lB2]^ (3.91)
>- AII ^>
ov . e
* r
and the ratio of A to AV is
50
-------
AV - 3H
Equation (3.90) can be approximated for both large and small interfacial
slopes, noting that a., is usually an order of magnitude greater than a .
For larger slopes (3 small), lateral dispersion is negligible and Eq. fs.90)
reduces to
(3.93)
)-,/.
The volumetric width of the mixed zone is
/"~t ,1/2
^-3 (3.94)
and the ratio A to AV is,
A
A₯ 2 1/2
er V*
For small interfacial slopes (large 0), longitudinal dispersion is negli
gible and Eq. (3.90) simplifies to
(3'96)
2 7TH
In this case the ratio of A to AV is given by
AV ~ 3U2, U.1/2
H (irnH)
1/2 1/2
J V (3'97)
Storage Cycle As explained in Section 3.2, the additional mixing during
storage due to grain dispersion is negligible. Suppose that at the start
of the storage period the volumetric length of the interface, AV, is
much smaller than the width of the dispersed zone, A. Then the rate of
tilting of the wedge will be governed by a coupled analysis which accounts
for the interaction of the density difference and the given dispersed
zone.
In this case it is assumed that the volumetric width of the dispersed
zone remains constant and equal to its initial value, A , throughout the
storage cycle. From the definition of A,
51
-------
A - = -t--\ = .(^\
o 3₯vc ' 27rnHr 3rvc '
o o
and from this expression the right hand side of Eq. (3.84) is evaluated
to yield
8Ui _ K Ap_ 2rmH
3y ~ n p AQ ri (3'98)
The velocity gradient is represented by
3u. _ 3r. , ,, _
3y 3t3y 27rnH 3tV.
where B = 3V±/3y and ₯£ = VQ + 3(y-H/2). Equation (3.98) then becomes
(3.99)
rrnHr. 3t r. A
ix o
and by introducing the approximation analogous to that preceding Eq.
(3.85), i.e., r. £ r or ₯. = ₯ ,this reduces to
d3
dt ~ A
o
If 3 is the inverse volumetric slope of the interface at the start of
the storage cycle , t=t , then
4imW
6 = S0 + A °(t-t ) (3.100)
which remains valid as long as 0H = A₯- « ₯ .
Withdrawal Cycle The coupled solution for the withdrawal period utilizes
the convective dispersion equation found in Section 3.2 (Eq. 3.69) and
its solution (Eq. 3.74). The analysis parallels that described for the
recharge period, and yields two equations analogous to Eqs. (3.85) and
(3.86) for the recharge period. The governing equations differ only by
signs, and the details of the development are omitted. The result for
the equation describing the flow dynamics (the equivalent of Eq. (3.85))
is
2
/9
'
52
-------
and the equation expressing dispersion (equivalent to Eq . (3.86)) becomes
j>- a-iQ_ i A /r,
-irrfe-SV
In these equations VA is given by
V* = V - Q (t-t )
* o xw o
The simultaneous solution to these equations is found as in the recharge
case and
c. 0.103)
6
where C, is a constant used to match with the previous cycle. The in-
verse slope, 3 , is found by substituting Eq. (3.103) into Eq . (3.101)
and integrating
2 4
1/2 32 a2 84 2ew 3 4ew 3/2
a (irnH) H r-r - + 7 ~- C, = -- =- V + C
(3.104)
The constant, C7, is dependent on the initial slope, 3
When considering multiple cycles, it may be convenient to evaluate the
constants C, and C7 in terms of the initial dispersed zone thickness,
A , and the initial inverse slope, 3 . Using V = irnHr 2} the dispersed
zone thickness becomes
3
3
A - (*.« - 2V- - , + - *.A<> (3.105)
w * o ₯A V
and the general expression for the inverse interfacial slope, 3, is
o 24(irnH) V. V
K o
ir !/2'
O " V ₯., V
0*0
53
-------
3TTH
IV*3/2 - V] (3.106)
The concentration field is obtained by using Eq. (3.105) in the general
solution as given by Eq. (3.74). The interpretation and application of
these results are discussed in the following sections.
3.4 Discussion of the Theoretical Results
In the previous portions of this Section, theoretical models were pre-
sented to describe the mixing of two fluids in a radial flow system, due
to density difference and grain dispersion. The first model presented
in Section 3.1 was based on an immiscible displacement process in which
the only mixing mechanism was the density difference. A nonlinear
partial differential equation (Eq. 3,11) describing the interfacial shape
and position for either recharge, storage or withdrawal was derived. An
equivalent equation was used by Perrine and Gay (1964) in their investiga-
tion of secondary recovery of oil. However, they were primarily inter-
ested in the effects of large viscosity differences on the interface and
did little analysis of those situations in which the density difference
is the predominant effect.
A solution of Eq. (3,11) was first developed for the recharge process.
It was found that the interface is described by a straight line in volu-
metric coordinates when the volumetric length of the interface, A₯=3H,
is small in comparison to the volumetric distance of the convected front
from the origin, V^. This is the case when e (Eq. 3.17) is small com-
pared to unity. For an initially vertical interface at the well this
ratio, AV/₯^, remains constant throughout recharge (Eq. 3.31) and the
interface is adequately described by the straight line interface at all
times during the recharge period if e is sufficiently small. When there
is some amount of tilting at the start of recharge (Eq. 3.52), such as
would occur after a recharge, storage, withdrawal sequence, the ratio
AV/V^ decreases as VA increases, Thus the description of the interface
by a straight line will improve with time.
In order to improve upon the linear interfacial solution found for re-
charge, perturbation techniques were used to investigate the curvature
of the interface. The results of this perturbation analysis, as shown
in Figure 3.2, indicate that there is only a small departure from the
linear interface, even for relatively large values of E . It is seen
that for large e the lower end of the interface has advanced into the
aquifer farther than it would for a linear interface. In the subsequent
analysis of the recharge process, it will be assumed that the interface
is linear in the volumetric coordinate system.
One of the primary characteristics of the immiscible model is the volu-
metric length of the interface. For the linear interface (lowest order
form of Eq. (3.31)),
AV = SH = 2e V, (3.107)
r *
54
-------
2
is shown in Figure 3.3 as a dimensionless plot of gD/Q H=AVD/Q H versus
2e T. The parameter x is a dimensionless time,
= K Ap_ _t = Dt
n p H H2
It is also useful to consider the interface as it appears to an observer
in a single observation well at a radius r. This is often the type of
information which may be gathered in the field. The observer_would
record depth averaged relative concentration in the aquifer, c/c , as a
function of time
c_ = I
c H
H c
7- dy (3.108)
c
o
Because the fluids a_re considered immiscible for the purpose of this
model, the value of c/c is identical to the dimensionless height of
the interface, n/H_-_ A simple method of representing these results is
the time slope at c/c =n/H=0.5 (i.e., V=VA) which is, from Eq. (3.51),
9n/H
= (At)"1 = -Q /3H
n/H=0.5 r
where At is the inverse of the rate of change of relative concentration,
It is convenient to present this result in terms of the dimensionless
variable
as obtained from Eq. (3.107). This result is shown graphically in Fig-
ure 3.4.
Equation (3.11) was also solved for storage and withdrawal processes. Dur-
ing storage the interface is described by a straight line in volumetric
coordinates when AV is small compared to the initial volume, V (Eq. 3.41).
The requirement of small AV/V must be met for the initial condition at
the start of storage. However, unlike the recharge process examined
earlier, the ratio AV/V is a function of time, V is a constant and A^=0H
is increasing with time (Eq. 3.49). If storage is continued long enough,
AV will eventually become large compared to V and the approximation of a
linear interface may become inappropriate. A reasonable limit on the ap-
plication of the linear interface solution (Eq. 3.41) is the condition
AV/V <1. From the results for recharge, it would be expected that inter-
face curvature will become important for AV/V near unity. Analysis of
these curvature effects is suggested as a topic for future study.
The withdrawal phase also yielded a linear volumetric interfacial solution
when the volumetric interfacial length, AV (Eq. 3.50), is small compared to
the volumetric distance of the convected front from the origin, VA. A
necessary, but not sufficient, condition for this solution is that e be
small compared to unity (e «1). The initial condition inherited from the
previous recharge or storage process must also be considered. Thus P H/V
55
-------
O"
o"
CO.
10
10
10
H/a = 75, a/ct
Immiscible Thy, Eq. 3.107
Coupled Thy, Eq. 3.90, 6/L>2
Coupled Thy, Eq. 3.90,
6/L<2
A Upper Limit of
Coupled Thy, 6/L=2
B Lower Limit of Im
Thy, 6/L=.5
C Upper Limit of Im
Thy, 6/L=.5
D Lower Limit of
Coupled Thy, S/L=2
10
Figure 3.3 Inverse Interface Slope, 3, as a
Function of Time during Recharge
56
-------
Immiscible Theory,
Eq. 3.109
Grain Dispersion
Effects, Eq. 3.113, -
6/L > 1
Coupled Theory,
Eq. 3.114, 6/L > 2
Coupled Theory,
Eq. 3.114, 6/L < 2
10
10
10
10
10'
10
2e T
r
Figure 3.4 Inverse Rate of Change of Relative Concentration,
At, as a Function of Time during Recharge
57
-------
must be small compared to unity at the start of withdrawal. Finally, it
is important to examine the temporal conditions under which the linear
interfacial solution is valid. The ratio A₯/₯^ can be expressed as
(Eq. 3.50)
As VA decreases to zero during withdrawal , the ratio A₯/₯^ increases and
the linear interface solution eventually becomes invalid. As in the sto-
rage case, the limit for this solution may be taken as A₯/₯^<1. Using
this condition in Eq. (3.110) there is the implied lower limit on the ap-
plication of the linear interface result during withdrawal,
(6 H/V )2 + 4e
2
w
When VA approaches this limit, it would be expected that interface curva-
ture effects may become significant. If a linear interface is assumed,
the lower end of the interface reaches the well when A₯/2=₯A. For this
condition A₯/₯A>1 and thus significant errors may be introduced when the
linear interface theory is applied at the well. The possible effects of
interface curvature as the interface approaches the well should be ana-
lyzed in future studies.
The effects of grain dispersion were analyzed in Section 3.2 by assuming
the grain dispersion to have no influence upon the flow dynamics. The dy-
namics are given by the immiscible model (Section 3.1) where the inter-
face is now idealized as the line of 50% relative concentration. The
volumetric interfacial length, A₯, was assumed to be small compared to the
distance from the origin, ₯A or ₯ , and a linear volumetric interface was
obtained. An estimate for the volumetric width of the dispersed zone, A,
was found for recharge (Eq. 3.65) and withdrawal (Eq. 3.76). The ratio of
the dispersed zone width to the volumetric interfacial length, A/8H, was
proposed as the criterion to estimate the importance of grain dispersion
during radial flow. For recharge this is given by Eq. (3.67) and for with-
drawal by Eq. (3.77).
As described in Section 3.2, this model is valid when A/3H is small com-
pared to unity. Figure 3.5 shows the ratio, A/gH, versus dimensionless
time, T, for several values of e during a recharge operation with an in-
itially vertical interface. This graph was constructed by assuming a
typical ratio of dispersivities (a1/a«=10) and by taking the ratio of
aquifer height, H, to longitudinal dispersivity, ct, , as a constant equal
to 75, a typical field value. A/$H is large for small time indicating
that interaction between density difference and longitudinal grain disper-
sion may be important near the well. A/3H decreases at first, but it
58
-------
Ui
VD
10.0
1.0
A/BH
0.1 i-
0.01
0.01
0.1
1.0
10.0
100.0
1000.0
Tlgure 3.5 A/3H as a Hinction of Time During Recharge for the Case of No
Dynamic Effects of Grain Dispersion
-------
eventually increases as lateral dispersion becomes important (Eq. 3.67).
It can be expected that even for large values of e , mixing due to grain
dispersion will eventually become as important as density induced mixing.
The exact limits of applicability of the immiscible model cannot be spe-
cified, but as was demonstrated in Section 2.4 reasonable estimates may be
obtained by assuming the limiting values of A/BH at a fixed number less
than one. A/8H=l/2 will be chosen here. For a given value of e , the in-
dicated range of dimensionless time, T, for which the immiscible analysis
is valid is found from Figure 3.5 or Eq. (3.67). For example, when e =0.1
the minimum 2e T is shown in Figure 3.5 as 0.35 and the maximum is 25. Be-
tween these two values the immiscible theory can be expected to apply.
The dimensionless length of time,AT, as defined in Eq. (3.109), is found
by differentiating the depth averaged concentration (Eq. 3.108) with re-
spect to time. The depth averaged concentration is obtained by integrat-
ing Eq. (3.64),
_ 2 _ 2
p 1 v/ZT 1 ^9 ^1
^* -*-r-it"~3r c c i -^~ / ^~ J- \ i i /1 -i i i \
= -^[H ;rrr 1 y ~erf U9~u ..erf u H (e -e ) ) J (3.111)
O /IT
where
n
= [V+(-l)nBH/2]//4T
The derivative of (Eq. 3.111) with respect to time is evaluated at ₯j.=V
(V=0)-
3c/c
o
3t
_ = (At)-A = --i_lerf(^fi) (3.112)
v=o 3H ot 2 A
For the recharge case (9V.fc/8t=Q ) with an initially vertical interface, BH
is given by Eq. (3.107) and A/$S by Eq. (3.67). Thus AT becomes
J~ \1^
AT = -^--| = - 2e T[erf r^ r- ]~X (3.113)
n p H r ^ ct2 (2e^T)
where
_
a 3?r
This result is sho;ra graphically in Figure 3.4. Note that for some mid
range values of time there is effectively no increase in AT due to the
addition of grain dispersion over the value predicted by the immiscible
model.
During the storage cycle, the additional mixing due to grain dispersion
was neglected. This assumption is consistent with the conditions neces
sary to derive the immiscible model with a linear volumetric interface.
60
-------
The model requires that the dispersed zone width, A, be small compared to
3H at the beginning of any process to which it is applied. During storage
3H must also be small compared to V for the model to apply. The dispersed
zone width at the start of storage, A, is proportional to the distance
traveled by the interface, V . During storage the interface will tilt about
its midpoint, with the result that the interfacial endpoints would travel
the most and undergo the greatest amount of di^v-r=ion. But this travel
must still be small compared to V in order fo< A₯/V ratio to remain
small. Thus it is consistent to neglect this additional dispersion.
In the third model, presented in Section 3.3, the combined effects of den-
sity difference and grain dispersion were analyzed by an approximate coup-
ling of the mass and momentum equations based on large A/3H. It was as-
sumed that the interface remained linear for this development. The inter-
facial slope for a recharge process with an initially vertical interface
is given by Eq. (3.90) and can be represented by a family of curves relat-
ing the volumetric inverse slope, 3, and V^ for various values of e , a, /a-
and H/ex. . Equation (3.90) is shown in Figure 3.3 as a dimensionless plot
for several values of e with a./a=10, and the ratio H/o, = 75. As noted
in Section 3.3, the coupled analysis is valid when A/3H»I. The ratio
A/&H (Eq.(3.92) with Eq. (3.90) substituted for 3) is shown in Figure 3.6
for a number of values of e , The range of applicability of the coupled
model (Eq. 3.90) can be analyzed by setting an approximate limit on the
value of A/gH. For this discussion the condition A/gH>2 is assumed for
the coupled model to apply. For example, when e =0.1 the coupled model is
applicable up to 2e T=0.018 during which longitudinal dispersion is impor-
tant. Later when 2e T exceeds 1000, the coupled analysis again becomes
applicable with lateral dispersion as the predominant mixing mechanism.
Figure 3.3 indicates these regions of validity and if closely examined,
it can be seen that the regions of validity for the immiscible model (Eq.
3.107) and the coupled model (Eq. 3.90) are mutually exclusive for a given
value of e . In general, there will be intermediate time ranges during
which neither model is applicable. In the example, e .l, these inter-
mediate ranges are 0.018<2e r<0.35 and 25<2e i<1000. The actual behavior
of the wedge in these ranges is not known,but a smooth transition between
the models is expected.
The curves representing Eq. (3.90) in Figure 3.3 are shown as dashed lines
for A/8H<2. It can be seen that, even for large values of e , longitu-
dinal dispersion will influence the tilting near the well, and if recharge
continues long enough, lateral dispersion will influence the tilting. As
e decreases, interfacial tilting becomes less and as lateral dispersion
becomes important, the rate of tilting decreases. The point of intersec-
tion of the immiscible theory and the coupled theory can be found by solv-
ing equations (3.90) and (3.107) simultaneously. That expression can in
turn be solved for the value of e for which the coupled theory is tangen-
tial to the immiscible theory at the point of intersection, e =(2T,a1/H)
(3a2/2a )1'2=3.26xlO~2 for Figure 3.3. For values of e greater than
this, there will be an intersection of the two theories; for values of e
less than this, the immiscible theory is never valid during recharge.
61
-------
100.Or 1 r
10.0
A/BH
1.0 -
0.1
10
10.0
100.0
Figure 3.6 A/BH as a Function of Time During Recharge for the Case of
Interaction Between Grain Dispersion and Density Induced Mixing
-------
The ratio A/3H indicates the relative importance of the two mixing mech-
anisms, grain dispersion and density difference. When A/8H is small (Sec-
tion 3.2), grain dispersion has no effect on the interface and the dynamics
are given by the immiscible model. When A/gH is large, the combined ef-
fects of the two mechanisms are important (see Appendix A for analogous
development in plane flow) and the coupled analysis governs the mixing.
The coupled solution does not reduce to the immiscible solution when
a ,a~-*Q because the regions of validity of the two solutions do not overlap.
The length of time, At, defined by the inverse time derivative of c/c
(Eqs. (3.108) and (3.111))at V=0 is given by Eq. (3.112) for both the°coup-
led and immiscible dynamic theories. For recharge (3V./9t=Q ) in which
interaction between grain dispersion and density difference is significant
3H is given by Eq. (3.90) and A is given by Eq. (3.91). Then the dimen-
sionless form of At can be written as
3/rT k
AT = KAp_At = _|^. [erf : ri (3ell4)
n p H Q H 7 J
2k^ (3D/Q Hr
4(1+ £ ^ )
k j (2£rT)
where
2TT ^ 1 1 / /, ^ / Q
_ n .I... -L/ 4 o/ o
^7 ~~ ^rr" 7v 7»^"T
This expression is presented graphically in Figure 3.4 for several values
of e . For values of A/3H less than 2, the curves are shown by dashed
lines. Note that for large values of e there is effectively no increase
in AT, for some midrange values of T, over that predicted by the immis-
cible theory, Eq. (3.109).
The combined effects of density difference and grain dispersion were also
investigated for the storage and withdrawal cycles. During storage the
change in dispersed zone thickness, A, was assumed to be negligible and
the rate of interfacial tilting was described by Eq. (3.100). For with-
drawal, expressions were presented which described the inverse slope of
the interface (Eq. 3.106) and the volumetric width of the dispersed zone
(Eq. 3.105). Some general features of these results may be noted. If
lateral dispersion predominates during withdrawal (Eq. 3.106 with a =0),
then it can be shown that A/3H always increases with time and the coupled
analysis can definitely become applicable. The same trend is noted for
the immiscible theory (Eq. 3.77), thus indicating that a transition from
immiscible dynamics to coupled dynamics can occur. If longitudinal dis-
persion is of primary importance during withdrawal, A/3H may increase or
decrease depending on magnitude of V ,3 and A (Eq. 3.106 with a =0).
If V »3 H, the indication is that A?3H°increases.
o o
The analysis of the interaction between density difference and grain
63
-------
dispersion, as presented in Section 3.3, has certain implicit limitations.
One of these is that A/3H be large compared to unity. This restriction is
explained in Appendix A for the case of linear two-dimensional flow and
the same restriction is implied for radial flow. As pointed out in Appen-
dix B, there is implied, from these approximate solutions, a relaxation of
the no-flux boundary conditions at the upper and lower boundaries of the
aquifer when A/3H is large. The possible effects of these boundary condi-
tions should be investigated in future studies.
3.5 Applications of the Radial Theory
One of the objectives of this research is to make it possible to forecast
the amounts and quality of fresh water which can be removed after storage
in a saline aquifer under various operational schemes. Procedures by which
the theoretical results in the previous sections can be combined to simu-
late complex operational conditions will be outlined. These procedures
should also be applicable to deep-well waste disposal operations in which
there is a density difference between the injected wastes and the native
water.
The immiscible theory (Section 3.1) can be used to construct simple, ap-
proximate results for many types of applications. To illustrate the use
of this theory, a one sequence operation of recharge, storage, and with-
drawal cycles will be analyzed. Additional sequences may be analyzed by
matching the solutions between cycles described in Section 3.1.
Consider a recharge operation with a flowrate, Q , for time, t , storage
for time, t , and withdrawal at a rate Q . For these conditions, we want
to find the interface position as a function of time and the amount of
water recovered. Utilizing the set of solutions described by a linear
volumetric interface, the volumetric length of the interface during re-
charge, 6H, is described by Eq. (3.52) with g =₯ =0. At the end of re-
charge when a volume V =Q t has been injected, the volumetric length of
interface is 3H=2e V which provides the initial value of 3 for the stor-
age cycle. During storage the net flow is zero and the interface is des-
cribed by Eq. (3.49), with 3 =2e V /H and after storage for a period t.
o r o s
(3H/VQ)2 = 4£
Using this value of 3 as the initial value, 3 , in Eq. (3.50) for with-
drawal, we find
(3H/₯Q)2 = 4e2(l+2ts/tr) + 4E2[l-(₯^/₯o)2] (3.115)
which can be used to evaluate the volume of water recovered.
For this analysis it will be assumed that withdrawal ends when the depth
averaged relative concentration, c/c , reaches some specific level c/c =b,
0
-------
based on geometric considerations of the interface as it reaches the well,
From Figure 3.7 it can be seen that when c/c =b at the well,
VA = gH(l-2b)/2 (3.116)
The total volume of water recovered at that time is
(3.117)
Note that this recovered volume includes some contaminated water with
"c7c =b. The recovery ratio, R=₯_/₯ , is, from Eqs. (3.116) and (3.117),
o K o
R = 1 - (l-2b)(BH/2V ) (3.118)
o
The value of (3H when withdrawal ends is found by solving Eqs. (3.115) and
(3.116) simultaneously and with that result the recovery ratio becomes
e 2+e 2(l+2t /t ) .
R = 1 - (l-2b) ( S / )1/Z (3.119)
l+(l-2b) e
w
This simple result can form the basis for approximate comparisons of al-
ternative recharge, storage and withdrawal schemes. As indicated in Sec-
tion 3.4, some error may be introduced because of interface curvature
near the well. However, the comparisons with experiments in Section 5
indicate that Eq. (3.119) provides reasonable estimates of R which are
somewhat less than those observed.
Some characteristics of the result in Eq. (3.119) can be seen by consid-
ering t =0 and b=0, in which case the recovery ratio can be represented
graphically as shown in Figure 3.8. Note that lines of constant R are
roughly circular in the e , e plane, hence, a given value of recovery
ratio can be attained approximately by specifying a constant distance
from the origin in this plane. The graph also illustrates the importance
of the parameters e and e in determining the recovery ratio. Favor-
able recovery ratios are generally obtained with low values of e and e .
Some other features can be analyzed using the immiscible solutions. For
multiple sequences of recharge-storage-withdrawal, interface movements
and recovery ratios can be found by repeated application of Eqs. (3.49),
(3.50) and (3.52) with appropriate matching between cycles. As a spe-
cific example of multiple recharge-withdrawal sequences, consider a well
operation in which e and e are unchanged (i.e., the same flox^ rates)
for several sequences and t =0. In each recharge cycle the same volume
of water is injected and during withdrawal, water is pumped until the
wedge first reaches the well (b=0). Under these conditions we find,
using Eqs. (3.49), (3.50) and (3.52)with matching between cycles, the
following recovery ratios for each sequence:
65
-------
1.0
e 0.5
r
I I i 1 1 1 1
0.5
1.0
w
Figure 3.8 Recovery Ratio, R, as a
Function of e and e (t =0)
r w s
Figure 3.7 Sketch of the Interface at
the End of Withdrawal
66
-------
- V *1 "
R = 1 + £ T - £ , £2 = £2 .. + £2(l+2£ .)
n n-1 n' n n-1 1 n-1
where recovery ratio R for the n' sequence of cycles is defined as the
ratio of the volume of water withdrawn during that sequence to the volume
injected during recharge in that sequence. The following examples il-
lustrate how the numerical values of the recharge ratio increase during
repeated sequences:
Sequence no. 12345
!Ll = 0.141: 0.859 0.928 0.942 0.949 0.953
!L1 = 0.526: 0.475 0.608 0.643 0.661 0.672
It can be seen that there is a significant increase in recovery ratio
during the first few sequences but that an upper limit on the recovery
ratio is approached asymptotically after several sequences. However, it
should be recognized that dispersive effects are cumulative and thus are
likely to alter these immiscible results after several sequences.
When the effects of dispersion are included, prediction procedures are
more complex and general analytical representations are not possible.
As discussed in Section 3.4, there are two different theories - immiscible
dynamics with superimposed dispersion (Section 3.2) or coupled dynamics
(Section 3.3) - which may be applicable depending on whether the ratio of
dispersed zone thickness to interface length, A/8H, is respectively much
less than or much greater than unity. In addition, there will be some
transition region in which neither of the theories is strictly correct
and transition from one theory to the other could concievably occur sev-
eral times during a recharge-storage-withdrawal sequence.
For purposes of specific prediction during a sequence, it will be assumed
that the immiscible model applies when A/$H<1 and the coupled model when
A/0H>1. Thus, the values of A and $H are computed from the appropriate
equations in Sections 3.2 or 3.3 and final values from the previous cycle
provide the initial values of the dispersed zone width, A , and the in-
verse slope, 3 , for the current cycle. The ratio A/3H is monitored
throughout the computation and the appropriate model (immiscible or coupled)
is selected using the current value of that variable.
67
-------
Table 3.1 summarizes the various equations that are used in the computa-
tional Trooe-.luro for a recharge-storage-withdrawal sequence. The table
is esse .: 1 :.y complete, but those recharge equations marked with an as-
terisk -.-.j.o.iad for the particular initial condition of a vertical
interface ^ \..= 0 with A =0. More general forms would be required in
some cases. For example, when an initial dispersed zone width, A > is
specified during recharge with an initial radius, r , and volume, V ,
Eq. (3.65) is generalized to
16Tm
A
A =
2
16i;e a~
r 2
H2
2
₯ r -₯
62H2
U 4* 2
r
2
r
(3.120)
The computational procedure outlined above is used in Section 5.6 to
develop comparisons between the miscible theories and the laboratory
experiments.
68
-------
TABLE 3.1 SUMMARY OF EQUATIONS APPLICABLE TO A RECHARGE, STORAGE, WITHDRAWAL OPERATION
A/6H1
e
(AV/H)
A
(V^TO
A
6H
Immiscible Model (Sections 3.1 & 3.2) Coupled Model (Section 3.3)
Recharge
Eq. 3.52
Eq. 3.65*
or
Eq. 3.120
Eq. 3.67*
Storage
Eq. 3.49
A
o
A /Eq. 3.49
J
Withdrawal
Eq. 3.50
Eq. 3.76
Eq. 3.77
Recharge
Eq. 3.90
Eq. 3.91*
Eq. 3.92
Storage
Eq. 3.100
A
o
A /Eq. 3.100
Withdrawal
Eq. 3.106
Eq. 3.105
Eq.3.105/Eq.3.106
For the case V = A = 0; a more general form (e.g., Eq. 3.120) is required in some cases.
-------
4. LINEAR FLOW EXPERIMENTS
The objectives of this part of the experimental investigation were to
evaluate the theoretical results for linear flow and determine any
limits of their validity. Experiments designed to evaluate the immis-
cible displacement model were successfully conducted. Due to model
capabilities and the availability of materials for the porous media,
experiments to evaluate the theory for the interaction between grain
dispersion and density induced mixing were not conducted. Addition-
ally, experiments were run to obtain data on the behavior of the inter-
face for small T.
4.1 Experimental Equipment
The linear experiments were conducted in a lucite box which contained
the porous media in a section 120 inches long, 3 inches wide and
6 inches deep (Figure 4.1),The porous medium was placed in the test
section through the top of the box. A sponge rubber strip was placed
between the medium and the top plate. The sponge rubber could expand
to prevent short circuiting along the top in the event that the medium
settled with time. The medium was kept in place by brass screens at
the ends of the box. Beyond the screens, the inlet and outlet reser-
voirs, which were free of porous media, provided constant heads over
the cross section. Inflow and outflow of the inlet and outlet reser-
voirs, respectively, were connected to constant head devices in order
to maintain steady flow rates during experiments.
Piezometer taps were located in the side wall at distances of 4.75 cm,
149 cm and 296 cm from the upstream or inlet screen. The upstream
screen was taken as the origin of the coordinate axes, with the x-axis
horizontal and at the mid-height of the medium. Three electrical con-
ductivity probes were located at x-^ = 62.4 cm, X2 = 182.9 cm and
Xo = 274.2 cm. The probes consisted of two parallel stainless steel
rods, 1/8 inches in diameter, spaced 3/8 inches on centers. The probes
averaged the conductivity over the depth of the medium to give a one-
dimensional record of concentration,
4.2 Procedure
After a careful packing, the porous medium was saturated with water,
before the rubber strip and top plate were installed, to prevent trap-
ping of air in the medium. In place measurements of hydraulic conduc-
tivity were made using fresh and salt water before and after each series
of runs.
Preliminary experiments were conducted to determine the longitudinal
grain dispersion coefficient for each medium used. Salt solutions of
specific gravity 1.0005, were used as tracers. For the initial
condition
t = 0, c = 0 x>0
c = c x < 0
o
71
-------
N3
I
Constant Head Devices
>J tfl
Fresh
Water
.1 Inle
ife
"*"
Inlet Valves
x
Probe 1
Probe 2
01
-------
the dispersion coefficient could be determined using the solution
(4.1)
where U is the steady seepage velocity. The slope of the concentration
versus time curve, for a constant x, at t = x/U = t,-n, is used to deter-
mine D.. from
Dl =
47T
where
50
c 3t
o
t=t
50
To satisfy the initial condition, the experimental model was initially
filled with fresh water, c = 0, and the flow rate set by observing the
manometer and adjusting the constant head devices. Then, the fresh
water flow was shut off and the salt water inlet valve and the reservoir
outlet valve were opened simultaneously. The fresh water in the inlet
reservoir was replaced by the salt tracer solution. Next, the reservoir
outlet valve was closed and the flow rate valve was opened to start the
experiment. This procedure was successful in introducing an initially
vertical interface, corresponding to the initial condition.
The conductivity probes were connected through an AC Wheatstone bridge
to a recording oscillograph to produce a continuous time record of
conductivity. Since conductivity was approximately a linear function
of tracer concentration, the conductivity records were interpreted
directly as concentration records and the necessary information was
obtained from them.
For the actual experiments, the model initially contained salt water
of known density p« and salt concentration c . Again the flow rate was
set by observing the manometers and adjusting the constant head devices.
The initially vertical interface was introduced by the procedure used
for the grain dispersion experiments, except in this case dyed fresh
water replaced the salt water initially in the inlet reservoir. Records
of conductivity were obtained in the same manner as for the dispersion
experiments. For runs with a density difference, the interface posi-
was observed visually and recorded periodically on the front of the box
with a china marker. These visual markings were recorded and the hori-
zontal projection of the interface, L, was determined for each marking.
73
-------
A series of experiments was conducted with the box containing a crushed
quartz sand having a mean grain diameter, d,_n, of 0.069 cm, and a uni-
formity coefficient of 1.37. The inplace hydraulic conductivity was
0.10 cm/sec and the porosity was 0.35.
4.3 Experimental Results for Linear Flow
Preliminary experiments to determine the longitudinal grain dispersion
coefficient produced results indicating that local heterogeneities were
present in the medium. The dispersion coefficients determined from
conductivity records at probe 3 were larger than those determined from
records at probe 2 which were in turn larger than those determined from
records at probe 1, Visual observations indicated thin layers having
higher hydraulic conductivity than the surrounding media. Thus, the
measured dispersion coefficients reflected the effects of this layering.
The largest indicated dispersion coefficient was used to characterize
the longitudinal dispersion coefficient, D , for the medium.
A number of experimental runs were made for specified values of the
density difference, Ap, and the seepage velocity, U. The parameters
for these_experiments are given in Table 4.1 along with the computed
value of E/D for each run. Smaller values of E/D could not be ob-
tained because of limiting model capabilities. The horizontal inter-
face projections, as determined by visual observations, are shown as a
function of time on the dimensionless plot in Figure 4.2. Also shown
are the theoretical expressions for L/H versus T given by Eq. (2.19)
and the interpolation equation (Eq, 2.72) given by Gardner, et al.,
(1962). The result given by Keulegan (1954)
L/H = /2r (4.2)
to represent his experimental results for an unconfined system with no
net flow is also shown in Figure 4,2, It appears that the presence
of a free surface reduces the rate of interface tilting. Gardner,
et al., (1962) also reported experimental results for models with no
net flow which confirmed their theoretical result (Eq. 2.72),
The_experimental results of the current study show that, in this range
of E/D , the rate of interface tilting is not affected by the net seep-
age velocity in these experiments. The current experimental results
do not provide a direct confirmation of the coupled analysis in Section
2.3 (Eq. 2.63) but are ^consistent with that theory. That is, for each
value of the parameter E/D.. the range in dimensionless time T is such
that the immiscible solution (Eq. 2.19) should be valid as discussed in
Section 2,4 and indicated in Figure 2.3, It was found that, using this
porous medium, it is not possible to operate this experimental equip-
ment in the range where grain dispersion becomes important. A model of
greater length or larger dispersivity would be required to observe these
effects.
74
-------
1 II I II
0.1
0.1
1.0 10 100
Figure 4.2 Horizontal Projection of Interface, L, as a Function of Time
-------
u
(cm/sec)
.008
.0101
.0085
.0309
Ap_
P
0,030
0.030
0.030
0,030
E
,2, .
(cm /sec)
0,0218
0,0218
0.0218
0,0218
"E7D
£./ JJ,
13.7
10.8
12.7
3.5
a = 0.2 cm, D. = Ua , E = - &- f
1 '1 1' n p 6
Table 4.1: Range of Test Results for the Linear Model
One-dimensional concentration, c, was recorded as a function of time at
the three conductivity_j>robes, for each experimental run. A typical
experimental curve of c/c versus T is shown, compared with the approxi-
mate one-dimensional concentration model given by Eq. (2.33) in Figure
4.3. A comparison of the results indicates a reasonable degree of agree-
ment between theory and experiment especially when it is noted that
Eq. (2.33) was based on an average dispersion coefficient over the mixed
zone.
76
-------
1.0
0.9
.7
.6
C/CQ .5
.4
.3
.2
0.1
Probe 1 (x=63.4 cm)
E/I>1 =12.7
Ap/p = 0.03
Experimental Curve
Theoretical Curve
Eq. (2.33)
2.0
3,0
4.0
5.0
6.0
7.8
8.0
9.0
Figure 4.3 A Typical Linear Flow, Experimental Breakthrough Curve Compared
with Theory
-------
5. RADIAL FLOW EXPERIMENTS
There were two principal objectives of the radial flow experiments:
(1) to construct a model with radial flow characteristics capable of
modeling real aquifer field situations; and (2) to evaluate the theo-
retical results and to determine the limits of their validity. The
theoretical results initially were to be evaluated for the immiscible
displacement model including evaluation of the separate topics of in-
jection, storage and withdrawal. Other tests were needed to evaluate
the effects of grain dispersion on the immiscible model and the combined
effects of density and grain dispersion on the total mixing process.
5.1 Model Construction and Equipment
The construction and operation of the radial model was more complex
than the linear model because radial flow from a line source is more
difficult to achieve than uniform flow in a constant area model. Also,
experience during operation of the linear model indicated that the choice
of media and the method of packing were important for measurement of dis-
persive effects.
Considerable care was taken to determine the model dimensions and media
to be used. Several porous media including porous foam rubber, crushed
quartz sand, and a consolidated media formed from plastic beads mixed
with epoxy resin were tested in a small linear model. The plastic beads
mixed with epoxy resin forming a rigid porous media were found to be
most suitable considering homogeneity and ease of construction. The
dimensions of the model were determined from practical considerations.
The necessity of packing the model uniformly while still leaving a model
thickness large enough to allow accurate interface measurements indi-
cated that a model thickness of three inches was desirable. Using this
depth and considering estimated interfacial tilting rates, it was deter-
mined that the model should have an eight foot radius. To keep the
model manageable with this radius, it was necessary to construct a sec-
tor model with a 15 degree included angle. These dimensions enabled
the modeling of a variety of real aquifer situations by varying the
density difference of the fluids and the flowrate.
Once the dimensions and choice of media were determined, the model was
constructed as shown in Figures 5,1 and 5.2, The model consisted of a
1/8 inch thick plexiglass sheet fastened to a 1 inch thick piece of
plywood, 2-1/2 feet by 8 feet, which acted as the base. The two radial
sides of the model were 1/2 inch thick plexiglass sheets, three inches
high, that allowed visual observation of dyed fluids in the model. The
sheets were placed at a 15 degree angle with the apex having a 1/8 inch
opening, simulating the well. Piezometer taps, 1/8 inch in diameter,
were drilled at ten locations in the model, including one within the
well itself (see Figure 5,1). The piezometer taps were connected with
plastic tubing to a piezometer board and a manometer. The piezometer
board consisted of open-ended plastic tubes (1/4"ID) fastened to a board
with paper marked at 1 mm intervals as a background. The manometer was
79
-------
Constant Head Tanks
00
o
To
Fresh-Water
Reservoir
Rotating Tube
Valve to Initiate lij
Vertical Interface [Jj
Wastewater
Well radius = 0.5 inches
Media radius = 7 ft 10 in
Media angle = 15 degrees
Piezometer Board
Collector
Manifold
To Salt-Water
Reservoir
Well
Flushing
Valve
Wastewater
Note:
Outlet
Valve
Total model length = approx. 10 ft
Model thickness = 3 inches
,lje.g., Location of piezometer number 1
Piezometer
Number
Radius
I (cm)
123456 78 9 10
1.27 31.3 61.8 123.1 123.1 123.1 199.5 199.5 199.5 239.0
Figure 5.1 View of the Radial Model
-------
Figure 5.2 Photographs of the Radial Model
During Operation
SI
-------
graduated in hundredths of a foot and could be read to a thousandth of
a foot accuracy with a vernier scale. Since it was more accurate than
the piezometer board, the manometer was used in all tests where piezom-
eter difference was measured. The apparatus was designed so that the
manometer could be connected between any two piezometer taps. The
piezometer board was used to give an overall measurement of the piezome-
tric head and to simultaneously compare the head at different locations.
With the piezometer taps drilled, the model was ready for packing.
Glass beads about one millimeter in diameter were sieved through a U.S.
Standard Sieve Series. Only those beads held by the no. 20 sieve after
passing the no. 18 sieve (1000 micron mesh) were used. This insured
uniformity of the media and gave a mean diameter (d,-n) of about 1 mm
and a uniformity coefficient of approximately 1.0. The beads were mixed
with epoxy resin (1 part resin to 25 parts beads by volume) and hand-
packed into the model. While the epoxy was hardening, 21 conductivity
probes were inserted vertically into the media (see Figure 5.3). The
probes consisted of two hypodermic needles (24 gauge stainless steel
tubing) located 3/4 inch apart on a circumferential line at various
radii from the well. A resistance measuring circuit, consisting of a
carrier amplifier system (Sanborn Co. Model 350-1100C), was connected
to the probes and then to a recording oscillograph (Sanborn Co. Model
296). The probes were inserted through the entire thickness of media;
hence, the recorder gave the depth averaged change in conductance as a
function of time at the location of each probe. These probes are here-
after referred to as "averaging probes". One of these averaging probes
was inserted inside the well. This probe served to indicate if the
start of a test had a step input with a vertical interface. During with-
drawal tests this-probe indicated the time at x^hich saltwater first
reached the well. This information was used later in determining the
percentage of water recovered during the test.
What shall be referred to as "point" probes were also inserted into the
top of the model. These probes consisted of two very short hypodermic
needles (16 gauge) inserted just through the top coat of the model,
1/16 inch apart. They were located at the same radius and adjacent to
averaging probes #6 and #10. These point probes were also connected to
a resistance measuring circuit and then to a recorder.
The 3/4 inch spacing of the averaging probes x^as originally chosen from
consideration of two factors. First, it was felt that they should be
this far apart so that during insertion of the probes the media between
them wouldn't be highly disturbed. Secondly, this spacing yielded sat-
isfactory current between the probes for anticipated salt concentrations
(1% to 2.5%). Experimentation showed that it would have been better if
the averaging probes were moved closer together. There was apparently
little disturbance to the media during probe insertion. In addition,
salt concentrations much lower than anticipated (less than 0.1% by
weight) were used for the tests. This meant that the current was less
than expected. Finally, a third factor was the radial width of the
electric field generated between the probes. With 3/4 inch spacing
82
-------
00
OJ
G>-©
Well
1) e.g., Location of Probe Number 1
Reservoir
Probe Number 12345678 9 10 11
Radius (cm) 1.27 19.6 34.9 45.7 50.0 67.7 87.0 104.5 104.2 122.9 142.0
Probe Number 12 13 14 15 16 17 18 19 20 21
| Radius (cm) 142.0 165.0 165.0 165.0 167.0 193.5 193.5 218.0 218.0 218.0
Figure 5,3 Location of Conductivity Probes
-------
between the probes, the electric field width (area that is sensitive to
the approaching interface) became significant and the probes read not as
points but as zones. This caused what can be termed electronic disper-
sion of the recorder output. For example, a vertical immiscible inter-
face approaching and passing the probes would yield a characteristic
I!S" shaped dispersion curve output rather than a step change as desired.
If the probes x^ere moved closer together, the thickness of this zone of
influence would be lessened and electronic dispersion decreased. Prob-
ably 3/8 inch spacing of the probes would have been preferable. A po-
tential theory analysis of the electronic dispersion for the 3/4 inch
spacing, however, showed it to be approximately an order of magnitude
less than hydrodynamic dispersion and it was disregarded during most
data analysis.
With the probes in place, the top surface of the model could be con-
structed. It x^as made by spreading a 1/8 inch thick epoxy coat onto
the exposed media. Prior to application the epoxy was allowed to become
very viscous so that it wouldn't flow into the media. It hardened into
a translucent top surface, bound to the media, eliminating short cir-
cuiting of the fluid flow along the top. The model was to be 7.62 cm high,
but due to seepage of the epoxy top coat in certain small areas, the
height was somewhat less. Visual observation of the seepage through
the side plates made it possible to estimate an average height of 7.46 cm
for the model.
The well end of the model consisted of concentric plastic tubes with
slots cut into them that acted as a valve during tests (see Figure 5.4).
This system was designed to insure a vertical interface between the
fresh and saltwater at the start of a test.
The well tube was connected through a "T" fitting and valves to saltwater
and freshwater constant head tanks (see Figure 5.1). The constant head
tanks could be raised above or below the well level in order to simulate
injection or withdrawal from the well. The valves and "T" connection
allowed a rapid change from freshwater to saltwater feeding of the well,
or vice versa.
The end of the model opposite from the well consisted of the media end-
ing in a circular arc of 7 foot 10 inch radius followed by a large
reservoir (see Figures.1). The reservoir was connected through a mani-
fold to a constant head tank. During an injection cycle, water was
drained from this tank, but during withdrawal fluid was constantly sup-
plied in order to maintain the head. The system of constant head tanks
at either end of the model produced hydrostatic conditions within the
well and the reservoir.
5.2 Experimental Procedure
Preparation of the Model for Operation Once the model was constructed,
it was tilted at an angle to the floor and saturated slowly from the well.
By tilting the model, air within the media had a chance to escape through
84
-------
00
Ul
To Constant Head Tanks
Rubber 0-Ring
1/2" I.D.
3/4" O.D.
Rubber 0-Ring
Flow into the model occurs only when
the slots line up
3/4" I.D.
1" O.D. }
To Exit Valve
and
Wastewater Tank
x^ ~^ ^
Slot
3" high
1/8" wide
-*
-V
N
'
--"
~PH.
X.
^
Figure 5.4 Detail of Well Construction
-------
the reservoir end as the water entering slowly displaced it. After the
model was saturated, it was leveled and prepared for experiments.
Deaerated water was used to saturate the model and run tests. This was
done to minimize the likelihood of gas coming out of solution in the
model forming air pockets. Deaeration was accomplished by pumping the
water through a chamber under a vacuum until the dissolved oxygen level
was around 3 mg/1. The water was also sand filtered and passed through
a demineralizer that removed anions and cations from the water. Finally,
to avoid the growth of either bacteria or algae in the model, a disin-
fectant, zephiran chloride, was added to the water.
Model Operation during Tests To start a test, a solution of sugar, salt
and water that had the desired density and conductivity was pumped into
the model. The freshwater constant head tank was set to the proper
height for the desired flowrate. About 30 seconds before the start of
the test, valves were opened allowing the freshwater to go from the con-
stant head tank into the well, flushing the saltwater out through an
exit valve (see FigureS.l). At time t = 0 the valve within the well was
turned so that the slots lined up and flow into the model could occur.
At the same time, the well exit valve was closed and a valve at the
reservoir end of the model was opened to allow flow through the model.
Flow was continued as long as desired, stopped for a period if neces-
sary, and could be reversed by changing the relative position of the
constant head tanks. This procedure was used to model recharge, storage
and withdrawal cycles.
Methods of Data Collection There were three types of measurements made
during tests: (1) measurement of flowrate; (2) visual measurement of
the interface shape and location; and (3) measurements made with the
conductivity probes. Flow measurements were made by using a 1000 ml
graduated cylinder and a stopwatch. Several readings were taken during
a test and the average of the readings was used as the flowrate.
Visual measurement of the interface consisted of periodic marking of
the interface location with a china marker on the clear plexiglass
sidewalls. The time of each marking was recorded. It was found that
as long as the density difference was in the range of principal interest
(from 1% to 2.5%), the dyed fluid made a distinct interface that could
be observed and whose location could easily be marked. When the den-
sity difference between the fluids was less than 1% by weight, the inter-
face was no longer sharply defined. Apparently, as the density difference
between the two fluids decreased, the relative importance of dispersion
and the effects of slight inhomogeneities in the model construction
increased. Since the experiments normally had higher density differences
than 1%, visual observation of the interface was possible.
In addition to the visual observation of the interface position, mea-
surement of the interface was made with the averaging conductivity
probes. In order to obtain a linear relationship between salt concen-
tration and output from the recorder, salt concentrations less than 0.1%
86
-------
by weight had to be used (see Figure 5.5). In order to achieve the
desired density difference of 1% to 2.5%, sugar was also added to the
salted water. The sugar was a nonelectrolyte and hence increased the
water density without affecting its conductivity. The sugar, however,
did have some affect on the viscosity of the denser fluid. Viscosity
increases of 2.5 to 7 percent are found for sugar concentrations of
from 1% to 2.5% by weight.
5.4 Determination of Model Parameters
The porosity of the model was determined in two ways. First, a value
of 0.337 was determined using a small test cylinder that had been packed
at the time of model construction. Its porosity was determined from
comparison of the volume of water required to saturate the cylinder
versus the gross volume of the cylinder. Since the test specimen was
packed in a similar manner to the model, with the same bead and epoxy
mixture, its porosity should be the same as that of the model. The
second method for determining the porosity was done using the in situ
model media. By measuring the travel time of fluid flow between con-
ductivity probes at different radii, an average porosity of 0.34 x^as
determined. This value is used in subsequent calculations.
The hydraulic conductivity of the media was also determined in two ways.
First, a cylindrical permeameter (1-1/2" diameter, 10" long) which had
been packed at the time of model construction was tested. Hydraulic
conductivity tests were made using this cylinder with the flow perpen-
dicular to the packed surface. The hydraulic conductivity from these
tests approximated the results that would be obtained by vertical flow
through the actual model. Consequently, the hydraulic conductivity of
the test cylinder was an approximation of the vertical conductivity of
the radial model. The average vertical conductivity from these tests
was 0.48 cm/sec at 60°F.
The horizontal conductivity of the model was determined from tests run
on the in situ media. The tests consisted of measuring piezometric
head and flow rate in order to calculate the hydraulic conductivity
from the well equation for a confined aquifer between two constant head
reservoirs. The equation is:
(5.D
where: K = horizontal hydraulic conductivity (cm/sec)
Q = flowrate (cm /gee)
s ^ piezometric difference between the radii
r and r2 (cm), (r,>r?)
H ^ thickness of the aquifer (cm)
The average horizontal conductivity derived from the use of this equation
was 0,43 cm/sec at 60°F (see Table 5.1). The fact that the conductivity
87
-------
CO
CD
V-i
QJ
4-1
n)
13
cn
H
p
C
H
ra
to
C
cu
CJ
t-l
0)
.10
.08 ~
.06
.04
.02
100 200 300 400 500
Recorder Deflection
600
700
Figure 5.5 Calibration of Averaging Conductivity Probes
-------
TABLE 5.1
HORIZONTAL HYDRAULIC CONDUCTIVITY
Note: The governing equation for the model:
r\ '
K =
Temperature of the water for all runs = 75°F
From the table below, the average K = 0.52
Average K = 0.43
sec
@ 75°F
sec
@ 60°F
Location
(Number refers
to piezometer)
1 to 10
1 to 8
1 to 5
1 to 3
1 to 2
2 to 8
2 to 5
3 to 8
Discharge*
Q' (cm-Ysec)
0.569
2.4
3.63
6.14
2.23
3.93
1.56
4.80
1.18
4,18
2.03
4.92
2.34
3.36
4.82
2.34
3.36
4.82
2.34
3.36
4.82
K
(cm/sec)
0.502
0.514
0.500
0.49
0.509
0.501
0.515
0.495
0.504
0.488
0.496
0.474
0.563
0.606
0.531
0.572
0.628
0.532
0.563
0.569
0.548
*This is the sector model discharge; i.e., Q = 24Q'
89
-------
measured in the test cylinder was very close to the in situ horizontal
conductivity of the model indicated that the media was reasonably
isotropic.
Since the model was to be used to study the effects of dispersion on
the mixing process between the two fluids at the interface, a measure-
ment of the media's dispersivity was required. The dispersivity was
determined through the use of the conductivity probes already mentioned.
As the interface between the two fluids passed a probe, the reading from
the recorder plotted the change in salt concentration versus time at that
point. There are several theories that mathematically describe the con-
centration c/c at a particular radius for flow from a well. One recent
study (see Gelnar and Collins, 1971) can be used to estimate the longi-
tudinal dispersivity of the media by considering the slope of the curve
c/c versus time at the point c/c =0.5 for the case of no density dif-
ference between the two fluids. Gelhar and Collins derived the following
equation assuming that the dispersivity, a , is a constant:
lrl
2u
1 r-.~.f /_
T err\
,16_
^ 3-,
2
r - i
3v ]
2
) r 4 r 4
m , * co . , 1/2
}]
(5.2)
where; c = concentration of tracer at radius r
c ^ original concentration of the tracer in the model
ou ^ longitudinal dispersivity
t = time
r = radius to point in question
rA = radius to the 50% point (=/2At )
r = radius of the well
D = coefficient of molecular diffusion in the porous medium
A = ^ where n = porosity, H = depth of aquifer, Q = flowrate
" /TTnH
The effects of molecular diffusion and the radius of the well are small
in magnitude compared to other terms and Eq. (5.2) reduces to
r2 - r 2
C 1 *^
co 2 ,16 3.1/2
(-^ v*}
(5.3)
Taking the derivative of Eq. (5.3) with_respect to time and evaluating
at c/c = 0.5 remembering that r - r. = 0 at c/c = 0.5,
o * o
3(c/c
c/cQ=.5
1677 a,
= S
50
(5.4)
Since t = rA /2A, we can solve for a in terms of the slope at the
90
-------
point c/c = 0.5 and r.
o
Experimentally, rA, is the radius to the probe in question and the slope
at c/c = 0.5 is determined from the recorder output.
Using this method of analysis, the dispersivity was determined from
several of the averaging probes as summarized in Table 5.2 and Figure
5.6. The longitudinal dispersivity generally falls in the range between
0.1 and 0.2 cm although at the larger radii, a few observations indicate
higher values. These large values of the dispersivity may be associated
with inhomogenieties of the media. The very small density difference
(see Table 5.2) may also have affected the results; however, for these
values of e this effect is indicated to be negligible using the results
of the coupled analysis in Section 3.3. The so-called electronic dis-
persion associated with the averaging probes may have introduced some
error in the dispersivity determination for some probes near the origin,
but, as discussed in Appendix C, this effect is unimportant over most
of the model.
A few determinations of dispersivity made using the more closely spaced
point probes were found to be of the same magnitude (0.075 cm < 04 < 0.18 cm)
as those found with the averaging probes. In all subsequent calculations,
the longitudinal dispersivity is taken to be constant throughout the
medium and equal to 0.15 cm.
5.5 Model Verification
The model was tested to see if radial flow was in fact occurring. Obser-
vations from the piezometers were used to obtain the distribution of draw-
down versus radius from the well. The results are plotted in Figure 5.7
and compared with the solution of the well equation, Eq. (5.1). Com-
parisons of piezometric head for the taps at the same radius showed no
measurable differences thus indicating a purely radial flow. The obser-
vations of piezometric head distribution in the model generally indicate
that the desired radial flow was attained and that the model was homo-
geneous with respect to hydraulic conductivity.
To further check the flow situation, observations of travel times were
made with practically zero density difference between the two fluids. A
small amount of salt was added to one fluid (0.02% salt by weight) in
order to obtain a measurable output from the conductivity probes. By
comparing the readings from probes located at the same radius, it could
be seen that the saltwater was reaching radial points at approximately
the same time. Also visual observation of the semiclear top coat of
the model during density tests showed that the injected fluid moved
91
-------
TABLE 5.2
DETERMINATION OF DISPERSIVITY FROM AVERAGING PROBES
Run Number
Discharge Q'cm-Vsec**
Density difference Ap/p
e
Probe
Number
2
3
5
6
7
8
9
10
11
12
14
17
18
21
Radius
cm
19.6
34.9
50.5
67.7
87.0
104.5
104.2
122.9
142.0
142.0
165.0
193.5
193.5
218.0
C-l
2.27
0*
0
1-2
4.66
0
0
1
4.94
0.001
0.039
3
2.96
0
0
5
3.40
0.001
0
Longitudinal Dispersivity a - cm
0.153
0.108
0.110
0.151
0.150
0.293
0.318
0.105
0.144
0.148
0.293
0.095
0.163
0.172
0.14
0.175
0.355
0.115
0.142
0.125
0.106
0.162
0.152
0.173
0.2
0.217
0.0927
0.146
0.157
0.405
The condition Ap/p = 0 indicates cases in which the density of
the liquids was matched using sugar and salt solutions to within
±0.0002 relative density difference.
**This is the sector model discharge; i.e., Q = 24Q'
92
-------
0,4
0,3
Run Number
A i
V 3
O 5
D 1-2
O c-i
LO
H
w
Q)
ex
tfl
H
O
no
0.2
0.1
I
1
0
0
V
1
20
60 80
Radius (cm)
100
120
140
Figure 5.6 Media Dispersivity Calculated From the Averaging Probes
-------
6 -
2.0
3.0
4.0
5.0
6.0
7.0
Figure 5.7 Piezometric Head Distribution in the Model
-------
radially into the model with a velocity independent of angular position.
To correlate the effectiveness of visual observation of the interface
location with the depth-averaging conductivity probe measurements, the
following comparison was made. The interface position was marked on
both of the side walls as the interface passed a conductivity probe. A
comparison was made of n/H as determined from the visual markings of the
interface with the relative concentration, c/c , from the probe output.
These values were recorded at several different times as the interface
passed the probe. The results for a run with a relatively large density
difference (Ap/p=0.020) are shown in Figure 5.8. This figure indicates
a good correlation between the probe reading and the visual recordings.
5.6 Experimental Results for Radial Flow
A series of experiments was undertaken to investigate the behavior of
density induced mixing in aquifers using the previously described visual
and electronic methods of data collection. Table 5.3 summarizes the ex-
perimental parameters along with the computed values of e and e for
each run. In some experiments, only recharge tests were conducted, while
in other experiments complete recharge, storage and withdrawal cycles
were investigated. When an item is not given in the table, it was not
measured during that particular run. A summary of model characteristics,
as used in data reduction, is given for convenience:
Model Height H = 7.46 cm
Porosity n = 0.34
Longitudinal Dispersivity a. = 0.15 cm
The hydraulic conductivity, K, was determined for each run using the
average viscosity of the two fluids (see Table 5.3).
The parameters e and e were found in Section 3 to characterize the re-
charge and withdrawal processes. In Figure 5.8 the results of a recharge
run with a relatively large value of e (=0.422) using both visual and
electronic interface measurements are compared with the predicted inter-
face position of the immiscible theory with a linear interface in volu-
metric coordinates (Eq. 3.107). The theory, which is based on small e ,
compares well with the experimental data. The indication is that the
interface tilting is governed predominantly by immiscible displacement
under these conditions and that the curvature of the interface is unim-
portant even for this relatively large value of e .
Electronic data were obtained in the form of breakthrough curves (see
Figure 5.12) for various conductivity probes during recharge and with-
drawal. A method was developed for interpreting these data_based on the
measured rate of change of depth averaged concentration at c/c =0.5.
This information was used to characterize the various breakthrough
curves and was expressed in terms of the parameter At:
95
-------
TABLE 5.3 EXPERIMENTAL CONDITIONS FOR THE RADIAL MODEL
Run
1
2
1-1
A-l
A- 2
A- 3
A- 4
A- 5
A- 6
A- 7
A- 8
A- 9
A-10
B-l
B-2
B-3
B-4
C-l
Ap/p1*
0.001
0.001
0.005
0.005
0.005
0.005
0.0215
0.0215
0.021
0.021
0.021
0.020
0.020
0.018
0.018
0.018
0.018
0.00014
Q '*
r
3
cm
sec
4.94
3.87
2.49
3.70
3.70
0.76
4.79
1.32
0.48
2.43
2.30
0.78
0.26
2.55
1.64
5.04
1.86
2.28
Q '*
w
3
cm
sec
5.74
5.55
6.07
4.70
2.70
1.86
5.39
1.96
Fluid
Temp.
°F
-71
=71
75
70
71
69.5
^70
^70
= 70
71
69
70
67
69.5
-69
^69
68
67
Average
V
centi-
stokes
0.995
0.995
0.940
0.985
0.977
0.985
1.012
1.012
1.012
0.995
1.020
1.012
1.050
1.012
1.010
1.010
1.020
1.050
Average
K
era
sec
0.49
0.49
0.52
0.49
0.50
0.49
0.48
0.48
0.48
0.49
0.47
0.48
0.46
0.48
0.48
0.48
0.47
0.46
e
r
0.020
0.043
0.123
0.098
0.099
0,216
0.177
0.337
0.556
0.248
0.250
0.422
0.726
0.222
0.277
0.158
0.257
0.020
e
w
0.079
0.080
0.155
0.179
0.216
0.260
0.153
0.251
Recharge
Time
sec
1500
1500
6000
2423
1440
1440
1440
2100
Storage
Time
sec
75
49
240
157
480
720
1200
10900
=1.000 gm/cm for all runs; - denotes estimated
temperature; Q ',Q =model
r w
flowrates, i.e.
-------
1.0
0.9
n/H 0.8
or
_. 0.7
c/c0
0.6
0.5
0.3
0.2
0.1
. | . | I | . | !
x .1 t_ i nni / T* _ O oo\
Immiscible lay (hq . J . Zz.)
~ t \ *\ \ ^
i * \ \ \
_ i x \
i \ \
\ v \ 1 \
; i ^ \ \ ^x
1 \ A V
k \ " \
- ^ \ \ X.
" \ \ X \
\ ^ v
)> Y \
, .A i v\ i , ix,
2.0 4.0 6.0 8.0
i i
O Probe 5
A ii ii
V Probe 6
T
A Probe 7
A II II
D Probe 9
O Probe 10
^k ii ii
Run A-9:
X
\
\x
^8
\
i i ^ ^ i
10.0 12.0
1 ' 1 '
Visual
Electronic
Visual
Electronic _
Visual
Electronic
Visual
Electronic
Visual -
Electronic
e = 0.422
r
-
-
X
e ^
! , I r^
14.0 1 h . 0 J 8 . 0
(KApt/npH)
L i.
Figure 5.8 Comparison of Visual and Electronic
Measurement of the Interface with
Immiscible Displacement Theory
-------
(At) 1 =
3(c/co)
at
c/c =0.
o
5
For a linear interface in volumetric coordinates c/c =0.5 coincides with
₯ =3s?* (r =rA) , where the subscript "p" refers to the°probe at which the
measurement is being made. The dimensionless form of At (Eq. 3.109)
found from the recharge experiments is shown in Figure 5.9 as a function
of dimensionless time, T (Eq. 2.14) and is compared with the theoretical
results presented in Section 3 (Eq. 3.109). The range of e and T is
such that the ratio of the dispersed zone width to interface length in
volumetric coordinates, A/$H, is small. Therefore only the immiscible
theory should apply for the most part. At low 2er however, there may
be some effects of dispersion as noted by the dashed line for e =0.1.
Figure 5.9 shows excellent agreement between theory and experiment,
although some systematic differences are noted in the early part of those
runs with larger values of e . Figure 5.10 shows a comparable visual
measurement of the length of the interface in volumetric coordinates,
A₯=j3H (Eq. 3.28), versus time for three runs with large e . Only the im-
miscible theory (Eq. 3.31) is shown, as e and T were in the range in
itfhich only density difference governs the mixing dynamics. Excellent
agreement between theory and experiment is again found, although some
systematic decrease in gH for smaller values of 2e T is indicated. A
glance at Figure 3.3 or Figure 3.4 indicates that this difference in At
or H for lower values of 2c T during runs of relatively high e cannot
be attributed to dispersive effects. Similar differences were observed
for the linear experiments presented in Section 4.3 and are attributed
to the influence of vertical velocities for small times, and the limited
applicability of the Dupuit assumption under those conditions. Figures
5.9 and 5.10 also show that, for this range of e and T, the rate of
interface tilting is not affected by dispersion and thus provide a di-
rect confirmation of immiscible displacement theory (Sections 3.1 and
3.2) and are consistent with the coupled theory (Sections 3.3).
Figures 3.3 and 3.4 of Section 3.4 indicate that density difference and
grain dispersion will interact to govern the mixing dynamics for small
time when longitudinal dispersion will be important or for very large
time when lateral dispersion will be important. It was found that it is
not possible to operate this experimental equipment in the range of 2e T
where lateral grain dispersion becomes important for recharge. A model
of greater length or larger dispersivity would be required. However,
several experimental runs were conducted for small e over the range of
times (2e T) in which the coupled theory should apply. The results are
presented in Figure 5.11; good agreement between theory and experiment
is indicated. The value a =0,15 cm was used for the plots of the coupled
theory. It is interesting to note that these runs were made with a very
small density difference (Table 5.3), and that several of them were also
analyzed to determine longitudinal dispersivity by ignoring this small
density difference (Table 5.2). In light of the importance of density
difference as indicated in Figure 5.11, it is possible that several of
the larger dispersivities found in Table 5.2 may be due in part to the
density effect.
98
-------
40 -
10 -
1
0.1 _
0.01
O A-l.A-2 0.098,0.099
0.216
0.177
0.337
0.556
0.248
0.250
0.422
Eq. 3.109
Eq. 3.90, e =0.1
.01
.1
10
2e T = 2e - ^- -
r r n p H
Figure 5.9 Comparison of Theory and Experimental Time
Dependent Concentration Slopes During Recharge
99
-------
100
Ci
<
10 _
Cf
Q
CQ.
1 -
O A-
A A-9
D A-10
.1
1 10 100
0 0 K Ap t
2e T = 2e
r r n p H
Figure 5.10 Volumetric Length of the Interface, BH=AV,
As a Function of Time
100
-------
1.0
0.1
0.01
10"
10
-4
10
Immiscible Theory, Eq. 3.109
-Coupled Theory. Eq. 3.114
/
-4
10" J 0.01
2e T= 2e K Ap t
r r n~H
0.1
Figure 5.11 Comparison of the Coupled Theory and
Experimental Time Dependent Concentration
Slopes During Recharge
101
-------
It has been shown that simple direct comparisons between experiment and
theory are possible for the recharge cycle with an initially vertical in-
terface and that excellent agreement is found. However, during storage
and withdrawal direct comparisons are more difficult because of the com-
plex dependence on conditions in previous cycles of the recharge-storage-
withdrawal sequence. The depth averaged concentration was monitored at
the various probes and at the well during withdrawal. Visual observa-
tions proved to be inadequate for these runs, which generally involved a
low value of e . In Figure 5.12 some typical breakthrough curves at the
well are shown in dimensionless form where ₯ is the net injected volume.
The curves are skewed, and the 50 percent point occurs earlier than
would be expected with a linear interface in volumetric coordinates.
These features may reflect some effects of interface curvature and local
mixing in the well.
A series of breakthrough curves for several probe locations during a
complete recharge-storage-withdrawal sequence^ (Run B-2) is shown in Fig-
ure 5.13. The depth average concentration, c/c , is plotted in terms of
a dimensionless time coordinate, (^n~^*)/^0> where V is the pore volume
between the origin and the probe. The increasing length of the mixed
zone is obvious for the recharge portion of the sequence, but during the
withdrawal process this trend is not clear.
Figures 5.14, 5.15 and 5.16 present a comparison of several observed re-
charge and withdrawal breakthrough curves with theoretical results from
Section 3 at various radii. The theoretical curves in Figures 5.14 and
5.15 represent recharge with an initially vertical interface. The in-
itial period during which the coupled analysis applies (A/8H>1) is short
and was neglected in the calculation. Therefore, the theoretical curve
represents the immiscible displacement theory with superposed dispersion
(Eq. 3.111 with 0H given by Eq. 3.107 and v^ given by Eq. 3.65). There
is excellent comparison between the theory and experiment at Probe 7
(Figure 5.15), especially for the leading edge of the curve (V^
-------
o
U)
c/c
-1.0 -0.8 -0.6
0.4 0.6 0.8 1.0
Figure 5.12 Experimental Breakthrough Curves at the Well
-------
c/c
1.0
0.9 h
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
-0.4
-0.3
Run B-2
// / Probe 3 Recharge
Probe 6 Recharge
Probe 3 Withdrawal
Probe 2 Withdrawal
-0.2
-0.1
0.1
0.2
0.3
(V -V.)/V (V = rrr nH, the pore volume to the probe)
p * o p p
0.4
Figure 5.13 Experimental Breakthrough Curves During Recharge and Withdrawal
-------
c/c
-.0] -.« -.01
Km *-7 rtou 1
o.jw -* . o.wai
I . I
Figure 5.14 Comparison of Experimental and Theoretical Breakthrough
Curves During Recharge
c/c
Figure 5.15 Comparison of Experimental and Theoretical Breakthrough
Curves During Recharge
c/c
I - 1 - 1 - 1
/
/
I -*^ I / I L.
HI A-7 F10BI 7
H« t O.U9 s1 0.*IS
1 . I 1 11.
Figure 5.16 Comparison of Experimental and Theoretical Breakthrough
Curves During Withdrawal
105
-------
1.0
₯
THEORETICAL CUKVE NEGLECTING DISPERSION
EXPERIMENTAL POINTS
O RECHARGE
A WITHDRAWAL
0 r.l 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0,4 .
0.2
Figure 5.18 At as a Function of Time During a Recharge-Storage-
Withdrawal Sequence, Run A-7
0.4
niiiiiiiir
I ' 1 1 1 1 i r
^0.15 cm Theoretical
a, - 0.30 cm CurV6S
0.3
-QAt
0.2
0.1
EXPERIMENTAL POINTS
O RECHARGE
A WITHDRAWAL
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Figure 5.17 At as a Function of Time During a Recharge-Storage-
Withdrawal Sequence, Run A-2
106
-------
of the dispersed zone. However, it is observed that At decreases system-
atically. A possible explanation for this decrease is the interaction of
density difference and grain dispersion. As shown in Section 3.3, this
interaction could decrease the rate of interface tilting, thus lower At.
This possibility was examined for Run A-2 in Figure 5.17. The longitu-
dinal dispersivity, a,, was doubled to a value of 0.3 cm and the theo-
retical curve recalculated. For this value of a, and for the given
parameters of Run A-2, fj gH remained at a value greater than one for the
entire recharge, storage and withdrawal process reflecting the importance
of longitudinal dispersion. The theoretical analysis is now based on
the coupled theory (Eq. 3.114) and is shown as a dashed line in Figure
5.17. It appears that the additional dispersion makes little difference
in the general trend of the theory, and that the decrease in At during
withdrawal must be due to some other factor.
It is of interest to evaluate the theoretical results of the immiscible
displacement dynamics with and without superimposed grain dynamics. Fig-
ure 5.18 presents two theoretical curves, one accounting for the superpo-
sition of dispersion (solid line) and one neglecting dispersion (dashed
line) for Run A-7. In fact, since both predict the same interface length,
it should be expected that the superposed curve would give higher values
of At. In both cases, predicted values of At during withdrawal are sub-
stantially larger than those observed. Other experimental values of At
are tabulated in Table 5.4.
These discrepancies between theoretical and experimental data during the
withdrawal cycle may be related to several factors. However, it is un-
likely that dispersive effects can account for the difference, for those
effects were evaluated in Figure 5.17 and were shown to be relatively
small. Most likely, the primary factor is the increasing importance of
the curvature of the interface as withdrawal proceeds. This is reflected
in Figures 5.17 and 5.18 in which the difference between theory and ex-
periment increases with time. The potential weakness of the theory with
a linear interface in volumetric coordinates was identified and discussed
in Section 3.4.
One of the ultimate purposes of this research is to predict the amount of
freshwater that could be recovered in a recharge-storage-withdrawal
sequence. The observed recovery ratios obtained by analyzing the break-
through curves at the well will be compared with the theoretical results
from Section 3.5. For an analysis of various breakthrough curves, it
was found that 10% relative concentration at the well was a reasonable
choice for the value determining the end_of the withdrawal cycle. The
observed recovery ratio, R=V /₯ , using c/c =0.1 at the x^ell is shown
in Table 5.5 along with the recovery ratio predicted by the immiscible
displacement theory (Eq. 3.119, b=0.1). Comparisons with the miscible
theories (Sections 3.3 and 3.4) are also included. These computations
were based on the matching procedure described in Section 3.5. Also in-
cluded in Table 5.5 are experimental results reported by Kumer (1968).
The predicted recovery ratios show trends similar to those observed in
this study and by Kumer, but are somewhat lower. These results indicate
107
-------
TABLE 5.4
OBSERVED VALUES OF -QAt/V DURING WITHDRAWAL
Probe
Number
1
2
3
5
6
7
8
9
10
Run Number
A-l
0.162
0.166
0.199
0.211
0.211
A-2
0.162
0.164
0.180
0.182
0.196
0.216
A- 6
0.693
0.785
0.760
A- 7
0.376
0.300
0.478
B-2
0.557
0.494
0.554
B-3
0.336
0.283
0.310
108
-------
TABLE 5.5 A COMPARISON OF THEORETICAL RECOVERY RATIOS WITH EXPERIMENTS
Run Number
Present ,
Experiments
A- 7
B-l
B-2
B-3
B-4
Kumar (1968)
Experiments^
2
3
6
e
r
0.248
0.222
0.277
0.158
0.257
0.334
0.459
0.647
e
w
0.179
0.216
0.260
0.153
0.251
0.296
0.309
0.647
Ap/p
0.021
0.018
0.018
0.018
0.018
0.373
0.223
0.092
Qr
(cc/sec)
58.32
61.20
39.24
120.96
44.64
9.98
3.05
0.8
QW
(cc/sec)
112.80
64.80
44.64
129.36
46.80
12.68
6.73
0.8
Average
K
(cm/sec)
0.49
0.48
0.48
0.48
0.47
0.0081
0.00785
0.0099
Recovery Ratio, R1
Observed
0.790
0.795
0.728
0.820
0.482
0.647
0.578
0.085
Immiscible
Theory
(Eq.3.119)
0.748
0.717
0.631
0.762
0.291
0.545
0.468
0.141
Miscible
Theory
0.700
0.680
0.595
0.720
0.265
0.515
0.440
0.132
1 2
Recovery based on c/c = 0.1 at the well, Recovery based on c/c
0.03 at the well.
-------
that dispersion actually decreases the predicted recovery ratio and does
not account for £he relatively large observed recovery. It appears that
the interface curvature may be the important effect in these cases. If
the behavior shown in Figure 3.2 for the recharge also occurs during
withdrawal, the arrival of the wedge at the well would be delayed. This
phenomenon during withdrawal could substantially increase the recovery
ratio and account for the decreasing values of At observed in Figure 5.17.
110
-------
ACKNOWLEDGMENTS
This investigation was conducted at the Ralph M. Parsons Laboratory
for Water Resources and Hydrodynamics of the Department of Civil
Engineering at Massachusetts Institute of Technology under a grant
from the Office of Research and Monitoring, Environmental Protection
Agency, Research Grant No. 16060ELJ.
The research was supervised by Dr. Lynn W. Gelhar, Associate Professor
of Civil Engineering, and while Dr. Gelhar was on sabbatical leave by
Drs. Donald R. F. Harleman and Keith D. Stolzenbach.
Mr. Roy Milley assisted in the construction of the laboratory models
and Mr. Edward McCaffrey assisted with electronic instrumentation.
Mrs. Pat Swan assisted in editorial work and typed the manuscript.
The use of the M.I.T. Information Processing Center is acknowledged.
The support and assistance of the Office of Research and Monitoring,
Environmental Protection Agency through Mr. Jack W. Keeley, Project
Officer, is gratefully acknowledged.
Ill
-------
REFERENCES
1. Bachmat, Y. and J. Bear, (1964), "The General Equations of Hydrodynamic
Dispersion in Homogeneous Isotropic Porous Mediums", Journal of Geophys-
ical Research, Vol. 69, No. 12, pp. 2561-7.
2. Bear, J. and 0. Dagan, (1964), "Moving Interface in Coastal Aquifers",
Journal of the Hydraulics Division, Proc. ASCE, Vol. 90, HY 4,
pp. 193-215.
3. De Josselin de Jong, G., (1960), "Singularity Distributions for the
Analysis of Multiple-Fluid Flow Through Porous Media", Journal of
Geophysical Research, Vol. 65, No. 11, pp. 3739-58.
4. Esmail, 0. J. and 0. K. Kimbler, (1967), "Investigation of the Tech-
nical Feasibility of Storing Fresh Water in Saline Aquifers", Water
Resources Research, Vol. 3, No. 3, pp. 683-95.
5. Gardner, G. H., J. Downie, and H. A. Kendall, (1962), "Gravity Segre-
gation of Miscible Fluids in Linear Models", Society of Petroleum
Engineers Journal, Trans. AIME, Vol. 225, pp. 95-104.
6. Gelhar, L. W. and M. A. Collins, (1971), "A General Analysis of Longi-
tudinal Dispersion in Nonuniform Flow", Water Resources Research,
Vol. 7, No. 6, pp. 1511-21.
7. Green, D. W. and R. Cox, (1968), "Storage of Fresh Water in Underground
Reservoirs Containing Saline Water, Phase II", Kansas Water Resources
Institute, Contribution No. 36.
8. Guymon, G. L., V. H. Scott and A. Hermann, (1970), "A General Solution
of the Two Dimensional Diffusion-Convection Equation by the Finite
Element Method", Water Resources Research, Vol. 6, No. 6, pp. 1611-17.
9. Hamrick, J. M., (1970), "Density Induced Mixing in Linear Flow Through
Porous Media", M.I.T. Department of Civil Eng., MS Thesis.
10. Harleman, D. R. F., P. F. Melhorn, and R. R. Rumer, (1963), "Dispersion
Permeability Correlation in Porous Media", Journal of the Hydraulics
Division. Proc. ASCE, Vol. 89, HY 2, pp. 67-85.
11. Harleman, D. R. F. and R. R. Rumer, (1963), "Longitudinal and Lateral
Dispersion in an Isotropic Porous Medium", Journal of Fluid Mechanics,
Vol. 16, Part III.
12. Henry, H. R., (1960), "Salt Intrusion Into Coastal Aquifers", Inter-
national Association of Sci. Hyd., Pub. No. 52, pp. 478-87.
113
-------
13. Keulegan, G. H., (1954), "Ninth Progress Report on Model Laws for
Density Currents; An Example of Density Current Flow in Permeable
Media", National Bureau of Standards Report No. 3411, Washington, B.C.
14. Kohout, F. A., (1970), "Reorientation of Our Saline Water Resources
Thinking", Water Resources Research, Vol. 6, No. 5, pp. 1442-48.
15. Kumer, Anil, (1968), "Dispersion and Gravity Segregation of Miscible
Fluids in Porous Media for Stratified Radial Flow Systems", MS Thesis,
Louisiana State University, Baton Rouge.
16. Kumer, A. and 0. K. Kimbler, (1970), "Effect of Dispersion, Gravity
Segregation, and Formation Stratification on the Recovery of Freshwater
Stored in Saline Aquifers", Water Resources Research, Vol. 6, No. 6,
pp. 1689-1700.
17. Li, W. H. and G. T. Yeh, (1968), "Dispersion at the Interface of Miscible
Liquids in a Soil", Water Resources Research, Vol. 4, No. 2, pp. 369-77.
18. List, E. J. and N. H. Brooks, (1967), "Lateral Dispersion in Saturated
Porous Media", Journal of Geophysical Research, Vol. 72, No. 10,
pp. 2531-41.
19. Mercado, A., (1967), "The Spreading Pattern of Injected Water in a
Permeability Stratified Aquifer", Symposium of Haifa, International
Assoc. of Sci. and Hydrology, Pub. No. 72.
20. Moulder, E. A., (1970), "Freshwater Bubbles: A Possibility for Using
Saline Aquifers to Store Water", Water Resources Research, Vol. 6,
No. 5, pp. 1528-31.
21. Perrine, R. L. and G. M. Gay, (1964), "Radial Motion of a Surface
Between Fluids Within a Porous Media", Journal of Canadian Petroleum
Technology, Vol. 3, No. 1, pp. 33-38.
22. Reddel, D. L. and D. K. Sunada, (1970), "Numerical Simulation of
Dispersion in Ground Water Aquifers", Colorado State Univ., Hydrology
Paper No. 41.
23. Rose, D. A. and J. B. Passioura, (1971), "Gravity Segregation During
Miscible Displacement Experiments", Soil and Science, Vol. 3, No. 4,
pp. 258-65.
24. Rumer, R. R. and D. R. F. Harleman, (1963), "intruded Salt-Wedge in
Coastal Aquifers", Journal of the Hydraulics Division, Proc. ASCE,
Vol. 89, HY 6, pp. 193-220.
25. Schiedegger, A. E., (1961), "General Theory of Dispersion in Porous
Media", Journal of Geophysical Research, Vol. 66, No. 10, pp. 3273-8.
26. Shamir, V. Y. and G. Dagan, (1971), "Motion of the Seawater Interface
in Coastal Aquifers; A Numerical Solution", Water Resources Research,
Vol. 7, No. 3, pp. 644-57.
114
-------
27. Shamir, V. Y. and D. R. F. Harleman, (1967), "Numerical Solutions for
Dispersion in Porous Mediums", Water Resources Research, Vol. 3, No. 2,
pp. 557-81.
28. Theis, C. V., (1967), "Aquifers and Models", Proc. National Symposium
on Ground Water Hydrology, American Water Resources Assoc., pp. 138-48.
115
-------
DEFINITION OF SYMBOLS
The primary symbols are defined below. A few additional
special symbols are defined where used.
Symbol Definition Dimension
A Q/27rnH (L2/T)
A Q /2-rrnH (L2/T)
r,w r,w v
b n/H at the well
c Concentration of material
c Concentration of material in liquid 2
c Depth average of concentration
C Matching constant between cycles (n=l,2,.--)
d,-n Mean grain diameter (L)
D KApH/np (L2/T)
2
D Longitudinal grain dispersion coefficient (L /T)
2
D9 Lateral grain dispersion coefficient (L /T)
2
Dm Molecular diffusion coefficient (L /T)
2
E Density induced dispersion coefficient (L /T)
2
EI A constant dispersion coefficient (L /T)
2
E Average density dispersion coefficient (L /T)
2
E Total (gross) dispersion coefficient (L /T)
f Relative interface height (n/H) (L/L)
F(z),F(2T) Dependent similarity functions (n/H)
2
g Gravitational constant (L/T )
H Aquifer height (L)
j(x,t) Constant of integration
2
j Mass flux of material relative to mean (M/TL )
Linear flow
_ 2
J Mass flux of material relative to linear (M/TL )
mean flow, for immiscible displacement
2
k Intrinsic permeability (L )
K Hydraulic conductivity (L/T)
L Interface length for linear flow (L)
n Porosity
2
p Pressure (M/T L)
117
-------
Symbol Definition Dimension
o
Q Recharge (+) or withdrawal (-) rate in (L /T)
radial flow
3
Q ,0 Recharge rate, withdrawal rate magnitude (L /T)
r Horizontal radial coordinate (L)
r. Radius to interface (L)
r^ Radius to convective front (L)
r Radius to probe (L)
r Well radius (L)
w
R Recovery ratio (=₯./₯ )
K O
s Piezometric difference between two points
S Interface slope for linear flow (=L/H)
S,-n Slo£e_ of concentration breakthrough curve (T~l)
5 at c/c =0.5
o
t Time (T)
t Time at the start of a process (T)
trQ Time from the start of a process till (T)
c/c =0.5 at radius r
o
u Horizontal seepage velocity in linear flow (L/T)
u1 0 Horizontal seepage velocity in liquid 1,2 (L/T)
1 sz
u. Velocity of the interface (L/T)
u Seepage velocity relative to mean flow (L/T)
U Mean linear convective seepage velocity (L/T)
v Vertical seepage velocity in linear flow (L/T)
v Horizontal radial seepage velocity (L/T)
v Vertical seepage velocity in radial flow (L/T)
y 2 i
V Volumetric coordinate (=7rnHr ) (LJ)
V. Volumetric distance to convective front (L3)
V = V - VA (L3)
V. Volumetric distance to the interface from (L )
1 V origin
₯£ = V - Vi (L3)
V Net injected Volume at the start of a (L )
0 process
o
V Volumetric distance to a probe in the (L )
P model
V Volume of liquid recovered (L )
R
118
-------
Symbol Definition Dimension
x Horizontal coordinate in linear flow (L)
x. Linear distance to interface (L)
x-L = x - x± (L)
3c = x - Ut (L)
y Vertical coordinate (L)
z,z' Independent similarity variables
a Longitudinal dispersivity (L)
<* Lateral dispersivity (L)
3 Inverse slope of the interface = L/H in (L2)
linear flow, = AV/H in radial flow
2
3 Inverse slope at the start of a process in (L )
radial flow
6 Width of the dispersed zone about the (L)
interface
A Volumetric width of the dispersed zone (L )
A Volumetric width of the dispersed zone at (L3)
the start of a process
At = 1/S5Q (T)
AV Volumetric length of the interface (L3)
Ap Density difference between liquids 1 (M/L )
and 2
AT = KApAt/npH
t, Distance to interface from top of aquifer (L)
n Distance to interface from bottom of (L)
aquifer
1/7 1/2
= B/V. ' in radial flow (L
'*
y Dynamic viscosity (M/LT)
v Kinematic viscosity (=p/p) (L2/T)
£ Transformed time coordinate
p Density (M/L3)
t>l}2 Density of liquid 1,2 (M/L3)
p Depth average density (M/L3)
T Dimensionless time (=KApt/npH)
e Parameter characterizing radial flow (=/D/A)
e =/D/A
r,w r,w
<£ Piezometric Head in liquid 1,2 (L)
'> ^
119
-------
APPENDICES
Page No.
A. Range of Validity of Approximation in Eq. 2.55 120
B. A Note on Boundary Conditions 123
C. An Estimate of Electronic Dispersion for the Averaging
Probes 126
121
-------
Appendix A
Range of Validity of Approximation in Eg. (2.55)
In Section 2.3 it was assumed that the interface would remain linear
(Eq. 2.38), and that the net effect of grain dispersion would be to
decrease the rate of interfacial tilting. Thereafter, the primary
approximation used to develop the coupled solution involved the appli-
cation of 3u/3y (Eq. 2.53) at the mean interfacial position in order
to determine 3u./9y (Eq. 2.56). The approximation
_3u
3y
x=x.
3u.
i
3y
(A.I)
was introduced in Eq. (2.55). The following analysis indicates when
this approximation is valid.
If the interface x.(y,t) is of a linear form as given in Eq. (2.38)
with 3 unspecified, the analysis of Section 2.2 applies and the density
gradient can be obtained from Eqs. (2.21) and (2.48) as given in
Eq. (2.58). Thus Eq. (2.53) can be written
_
3y
K
n p
(A.2)
In this equation x. = x.(y,t) is the only function of y, and the equation
has the integrated form
u =
K Ap 1
n p 23
erf
x-Ut
- erf
x-x.
1
/4T
(A.3)
where j(x,t) is a constant of integration.
From the continuity equation (Eq. 2.51)
f y a.
v = -
j
r dy
3x J
y=-H/2
which satisfies the condition v = 0 at y = -H/2. The requirement that
v = 0 at y = H/2 implies that
3
3x
.H/2
-H/2
u dy
= 0
122
-------
Expressing the velocity in terms of the convective coordinate system,
u = u - U
this flux condition becomes
l-H/2
'-H/2
u dy = 0
This condition is used to evaluate j(x,t), with the result
= K Ap 1
U n p 28
H
+H/2 x-x.
-H/2
i . X-X.
erf dy - erf
(A.4)
Along the interface (x=x.), the_interfacial velocity relative to the
convective front is defined at u.,
u. = u. - U H u
i i
x=x.
(A.5)
which is combined with Eq. (A.4) to yield
_
u .
1
K Ap
n p
1
2BH
r+H/2
J-H/2
s±
erf
A
:r dy
£
ix=x.
(A.6)
Equation (A.6) is evaluated using
-z
erf z dz = z erf z +
const
and can be expressed as
p (-282H)
(A. 7)
where
- B(y-H/2)
- B(y+H/2)
(A. 8)
123
-------
Differentiating this equation with respect to y yields
^i K Ap -1
(A.9)
The approximation in Eq. (2.55) involves replacing quantity
jhi
3y
x=x.
i
= K A£
n p
(A.10)
from Eqs. (2.54) and (2.60) by the expression in Eq. (A.9). This will
be valid when
(A.11)
in which case, using erfy . = y.,
3y
n p
/47T£
8U
3y
(A.12)
x=x.
i
The maximum absolute values of y and p_ occur at y = ±H/2 and thus the
restrictions on y.. ,y» are seen to be equivalent to the requirement that
/4T
« 1
Since the approximate width of the dispersed zone about the interface,
5, is given by ATT £ (Eq. 2.47) the restriction implied by the approxi-
mation in Eq. (2.55) is L « 6. Therefore, the result obtained in Sec-
tion 2.3 is expected to be valid when the dispersed zone is wide in
comparison to the wedge length.
124
-------
Appendix B
A Note on Boundary Conditions
For a confined aquifer of constant thickness, the boundary condition
to be satisfied by the solution of the convective dispersion equation
at the upper and lower boundaries is
= 0, y = ± H/2 (B.I)
The approximate solution to the convective dispersion equation for the
appropriate boundary conditions considered in Sections 2.2 and 2.3 was
as given by Eqs. (2.44) or (2.58). For this solution
R (x-x )
-> (B'2)
where x-x. = x-gy. At the lower boundary, y = -H/2, and Eq. (B.2)
becomes
The consequence of not satisfying the boundary conditions is the crea-
tion of an artificial mass flux through the aquifer in the vertical
direction. At y = -H/2 this flux is proportional to Eq. (B.3).
The possible local effects of not satisfying this boundary condition
can be evaluated by comparing the magnitude of the ficticious vertical
flux to the horizontal dispersive mass flux at a given point. The
ratio of these two fluxes is, from Eqs. (2.58) and (B.3),
D2 If /D1 If = 6VD1
Thus, it is seen that, even though the ratio D/D.. is generally small
"-*- or less), the implied vertical flux of mass at the boundary will
eventually become equal to or larger than the horizontal flux because
3 increases with time as indicated in Eqs. (2.63) and (2.68). The im-
plication is that the effects of the no-flux condition at the boundaries
may become important when 3 > D-,/!^-
Another method of estimating the possible effects of the boundary con-
ditions is to compare the horizontal and vertical components of the
total mass transport through the aquifer. Hamrick (1970) computed the
artificial vertical flux and compared it with the total horizontal flux
125
-------
evaluated at its maximum, x = 0. The total vertical dispersive mass
flux is given by
dx =
(B.5)
and the maximum horizontal masj^ flux relative to the mean motion is,
assuming p to be constant, at x = 0,
r+H/2
'-H/2
8x
fH/2
_ dy + p (cu)
x=0 '-H/2
(B.6)
o 1
/7rL, 3 E A/L
) +-_*(-
where the integral
r
erf
TJ =
H
Lw/rT
(B.7)
may be evaluated numerically. However, for small . , the error function
can be approximated by the first term of its series expansion and $ is
given approximately as
(B.8)
At its maximum, within the aquifer
26
Under these conditions the error function in Eq. (B.6) is approximated
by erf(/TrL/2<5 ^ L/6. Thus, Eq. (B.8) is reasonably accurate for the
range of 6/L for which the interaction solution is applicable (6/L»l).
As a consequence J can be approximated by
X
T - PCoDlrL 1 E ,1.2.
Jx ~ 3 16 + 2 D(T) J
(B.9)
The ratio of the absolute values of the total mass flux is given by
126
-------
J
JL
J
x
2 6
(B.10)
Again, it is seen that the implied ratio of vertical to horizontal mass
transport will eventually become large, in this case because both 6/L
and 3 increase with time. Thus, from consideration of both local fluxes
and total mass transport, there is the implication that the boundary
conditions may, for sufficiently large times, influence interfacial
tilting and overall mixing for the coupled analysis (Section 2.3).
127
-------
Appendix C
An Estimate of Electronic Dispersion for the Averaging Probes
The following analysis is based on the solution for the concentration in
a uniform linear flow with no density difference between the two fluids,
, 1,1-,X-Ut
c/c = -=- + -= erf (
o 2 z /rtr
(C-l)
in which x is the longitudinal coordinate, u is the seepage velocity and
D is the longitudinal dispersion coefficient. It will be assumed that
a probe at x = x1 measures the average concentration in a region x.>a
around the probe because of the generated electronic field. The indi-
cated average concentration, c, from the probe readout becomes
=x +a
- = (l/2a)(l+erf
'4D t
dx
(C-2)
and after integration
4a
x+a-ut
,+a-ut
erf
x.. a-ut x..-a-ut
erf
/4D t
4D t
2 2
-(x +a-ut) -(x -a-ut)
4D-t
VTT
4D,t
e 1
(C-3)
From Eq. (C-3) we find that
c
(At)
-1
r-
= - - erf
X;L=ut
(C-4)
and as a-K), Eq. (C-4) reduces to
(At)
9
-1 Co
3t
u
(C-5)
128
-------
xjhich corresponds to a point measurement with no electronic dispersion.
The effects of electronic dispersion can be assessed by considering the
ratio of the indicated and actual rates of change of concentration (Eqs.
(C-3) and (C-4))
At/At = (6/2a)erf
where 6 is the dispersed zone thickness defined by
(C-6)
i£
3x
,-1/2
x=ut
Figure C-l shows this ratio as a function of a/5. It can be seen that
when a/6 < 0.2, the indicated and average rates of change differ by
less than 5%. The possible error in At is of importance because the
values of At from experiments are used to determine dispersion
coefficients.
This analysis can be adapted to the observations in radial flow if we
determine 6 during recharge from Eq. (5.3)
-1/2
and estimate that 2a is equal to the spacing of the probe electrodes
(a=3/8 in.). Using these relationships and the estimated dispersivity
from Section 5.4 (a =0.15 cm), the radius rA corresponding to a given
value of a/6 can be determined as shown in Figure C-l. Probe locations
are also shown.
From Figure C-l it can be seen that during recharge the distortion due
to electronic dispersion will be significant for only the first few
probes near the wall. For tests with a density difference or measure-
ments during withdrawal, this distortion will be unimportant because the
thickness of the dispersed zone will be larger under these conditions.
129
-------
1.0
0.8 -
0.6 -
At/At
Equation C-6
0.4
0.2
c
o
H
u o\
ra =>fc
o
o o
0)
o
=«=
=«=
J.
_L
0.1
0.2
0.5
0.6
0.3 0.4
a/6
Figure C-l Estimate of Electronic Distortion During Recharge
0.7
-------
SELECTED WATER
RESOURCES ABSTRACTS
INPUT TRANSACTION FORM
/. Report No.
2.
DENSITY INDUCED MIXING IN CONFINED AQUIFERS
7. Author(s) Gelharj L} wilson, J.L., Miller, J.S., and
Hamrick, J.M.
9. Organization
Massachusetts Institute of Technology, Cambridge,
Department of Civil Engineering
3. Accession No.
w
5. Report Date
6.
8. Performing Organization
Report No.
10. Project No.
EPA. WQO. 16060ELJ
12. Sponsoring Organization
15. Supplementary Notes
11. Contract/Grant No.
13. Type of Report and
Period Covered
Final Report, EPA, WQO, Project #16060ELJ, 1972.
137 p, 36 fig, 7 tab, 28 ref, 3 append
16. Abstract .
Analytical techniques are developed to describe the mixing between two fluids of dif-
ferent density in a confined aquifer, in which one fluid is introduced to the aquifer
by well recharge. The immiscible displacement process in both linear and radial flows
is analyzed and the effects of longitudinal and lateral dispersion are included using
a boundary layer approximation. The theoretical results demonstrate the effect of hy-
drodynamic dispersion in retarding gravity segregation due to density differences. The
theoretical results are compared with observations of aquifer mixing in linear and ra-
dial flow laboratory models. During recharge excellent agreement between the theore-
tical predictions and experimental results is found and the predicted retarding effect
of longitudinal dispersion are verified. During withdrawal some systematic difference
between the theory and observation are noted. Theoretical predictions of recovery ef-
ficiency during a recharge-storage-withdrawal sequence show trends similar to those ob
served but are typically 5 to 10% lower than those observed. Direct theoretical pre-
dictions of recovery efficiency during single or multiple sequences of recharge-
storage-withdrawal are developed for an immiscible system, and similar developments
outlined for miscible displacement. The results are suggested for use in developing
preliminary estimates of operating conditions in the field. This report was submitted
in fulfillment of Project Number 16060ELJ under the sponsorship of the Water Quality
Office, Environmental Protection Agency. (Gelhar, MIT)
17a. Descriptors
Groundwater Recharge, Water Storage, Dispersion, Density Currents,
Aquifer, Water Quality, Recharge Wells, Underground Storage
*Confined
17b. Identifiers
17C.COWRR Field & Group Q4B 05E 02F
18. Availability
19. Security Class.
(Report)
20. Security Class.
21. No. of
Pages
22. Price
Send To:
WATER RESOURCES SCIENTIFIC INFORMATION CENTER
U.S. DEPARTMENT OF THE INTERIOR
WASHINGTON. D. C. 20240
Abstractor Lynn W. Gelhar
^Institution Massachusetts Institute of Technology
WRSIC 102 (REV. JUNE 1971)
6PO 9 J 3.2«f
------- |