UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
ANALYSIS OF FRESHWATER LENS FORMATION IN SALINE AQUIFERS BY JOHN PATRICK GLASS A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1979 ACKNOWLEDGEMENTS This research project was initiated by Dr. B. A. Christensen and Dr. Hillel Rubin. It was supported financially by the Office of Water Research and Technology of the U.S. Department of the Interior. The author extends his gratitude to Dr. B. A. Christensen for his continued support and advice and for serving as Chairman of his Supervisory Committee. Special thanks are due to Dr. Joel Melville whose daytoday suggestions and assistance contributed substantially to the work. Above all, I am thankful for the assistance of my wife, Claudette, who participated in the experimental work, helped with the preparation of the figures, typed the entire manuscript and supported me financially throughout my graduate study. TABLE OF CONTENTS ACKNOWLEDGEMENT . LIST OF TABLES LIST OF FIGURES . . . LIST OF SYMBOLS . . . ABSTRACT . . . CHAPTER 1 INTRODUCTION . . 2 LITERATURE SURVEY . . General . . Studies Dealing with Interface Location Studies Dealing with Dispersion . Studies Concerned Directly with Freshwater 3 BASIC EQUATIONS AND FLUID PROPERTIES . Basic Equations . . Darcy's Law . . Continuity . . Well Flow in Leaky Aquifers . The Dispersion Equation . Boundary Conditions at an Interface . Fluid Properties . . Density . . Viscosity . . 4 THE HELESHAW MODEL . . General . . The Viscous Flow Analogy . . Model Scaling . . Model Construction . . Properties of the Viscous Fluids . Operation of the Model . . . . 52 Page ii . . . . v Lenses . . * . * . Analysis of Test Results TABLE OF CONTENTScontinued Page 5 SIMPLE CASES INVOLVING STEADY FLOW . . General . . Freshwater Injection in a Coastal Aquifer . Saltwater Intrusion . . Injection Into a Coastal Aquifer . Confined coastal aquifer . . Phreatic coastal aquifer . . Axisymmetric Flow in a Leaky Aquifer . 6 A FINITE ELEMENT MODEL FOR THE STEADY FLOW CASE . General . . . Description of the Galerkin Finite Element Scheme . The Fundamental Concept . . The Element and its Shape Functions . . Finite Element Formulation of the Governing Equation Development of the Element Matrices . . Computer Implementation . . Program Verification . . . Freshwater Lens in a Leaky Saline Aquifer . 68 69 69 70 74 80 82 94 94 94 94 95 99 105 117 118 122 7 NUMERICAL SIMULATION OF LENS GROWTH . .. 127 General . . Description of the INTERCOMP Model . . Model Verification . . . Lens Formation in Aquifers With Weak Vertical Flow . Lens Formation in Aquifers With Strong Vertical Flow Dimensionless Representation of Lens Growth . The Influence of Dispersion on Lens Shape . Changes in Lens Shape During Recovery ... 8 CONCLUSIONS . . ... .... 149 APPENDIX . . . . . 152 FINITE ELEMENT PROGRAM DOCUMENTATION . ... 152 REFERENCES . . ... . 168 BIOGRAPHICAL SKETCH . . . . 172 LIST OF TABLES Table Page 4.1 Summary of Similitude Requirements for Vertical HeleShaw Model With Two Immiscible Fluids ........... 37 6.1 Shape Functions for the 20 Nodes Element ............... 101 6.2 Partial Derivatives of the Shape Functions ............. 108 7.1 Fluid and Aquifer Properties Used With the INTERCOMP Model ...................................... 133 LIST OF FIGURES Figure Page 3.1 Definition Diagram for Well Flow in a Leaky Aquifer ...... 15 3.2 Temperature vs. Density for Water with Various Salt Concentrations ......................................... 22 3.3 Water Density vs. Salt Concentration at Various Temperatures ......................................... 23 3.4 Viscosity vs. Temperature for Freshwater and Saltwater ... 25 4.1 Schematic of Gap Between Plates ........................ 28 4.2 Glass Plates with Infusion Reservoirs and Seal .......... 40 4.3 Assembled Model on Its Stand ............................ 41 4.4 Back Side of Model with Constant Head Tanks .............. 42 4.5 Infusion Pump Connected to Model ........................ 43 4.6 Temperature vs. Viscosity for the Two Viscous Fluids ..... 45 4.7 Test Run With Horizontal Model. Injection Rate = 20.44 cc/min ......................................... 48 4.8 Test Run With Horizontal Model. Injection Rate = 41.04 cc/min ......................................... 49 4.9 Instability Developing After 50 Seconds. Injection Rate = 81.72 cc/min ....................................... 50 4.10 Unstable Interface After 3 Minutes. Injection Rate = 81.72 cc/min ....................................... 51 4.11 Interface After 1 Minute. Injection Rate = 20.44 cc/min ......................................... 53 4.12 Interface After 3 Minutes. Injection Rate = 20.44 cc/min ......................................... 54 4.13 Interface After 5 Minutes. Injection Rate = 20.44 cc/min ......................................... 55 4.14 Interface After 7 Minutes. Injection Rate = 20.44 cc/min ......................................... 56 LIST OF FIGUREScontinued Figure Page 4.15 Interface After 9 Minutes. Injection Rate 20.44 cc/min ......................................... 57 4.16 Interface After 11 Minutes. Injection Rate 20.44 cc/min ......................................... 58 4.17 Interface After 10 Seconds. Injection Rate = 8.16 cc/min Model Depth Reduced by half ............... 59 4.18 Interface After 1 Minute. Injection Rate = 8.16 cc/min Model Depth Reduced by Half ............... 60 4.19 Interface After 2 Minutes. Injection Rate = 8.16 cc/min Model Depth Reduced by Half ................ 61 4.20 Interface After 3 Minutes. Injection Rate 8.16 cc/min Model Depth Reduced by Half ............... 62 4.21 Interface After 4 Minutes. Injection Rate 8.16 cc/min Model Depth Reduced by Half ................ 63 4.22 Interface After 6 Minutes. Injection Rate = 8.16 cc/min Model Depth Reduced by Half ................ 64 4.23 Interface After 8 Minutes. Injection Rate = 8.16 cc/min Model Depth Reduced by Half ................ 65 4.24 Comparison of Experimental Results With Equations for Interface Rotation ..................... ......... 67 5.1 Shape of Saltwater Wedge (Dupuit Parabola) ............... 71 5.2 Injection Into a Confined Coastal Aquifer ................ 72 5.3 Freshwater Lens Produced by Injection Into a Confined Coastal Aquifer ......................................... 77 5.4 Dimensionless Representation of Interface Toe, Q' = Q/ = 10.0 ......................................... 78 5.5 Dimensionless Representation of Interface Toe, Q' = Q/R = 5.0 ........................................ 79 5.6 Injection Into a Phreatic Aquifer ...................... 81 5.7 Steady Freshwater Lens in a Leaky Saline Aquifer ......... 84 5.8 Differential Wedge of Leaky Aquifer ..................... 87 LIST OF FIGUREScontinued Figure Page 5.9 Integral Curves for Equation (5.47) ...................... 90 2 5.10 Q6/ (27Kb ) vs. rt/B for Three Values of a/b ............. 92 5.11 Section Through Steady Freshwater Lens for Two Values of Leakage, B,: a = 0 ................................ 93 6.1 ThreeDimensional Element Showing Node Numbers ........... 96 6.2 Distribution of Two Quadratic Shape Functions Over a Two Dimensional Quadrilateral Element .................... 98 6.3 Local Coordinate System for the ThreeDimensional Element ........................................... 100 6.4 Sample Points for Gaussian Quadrature of Equation (6.23).. 112 6.5 Sample Points and Weight Factors for Gaussian Quadrature of Equation (6.25) ..................................... 114 6.6 Sample Points for Gaussian Quadrature of Equation (6.25) a = 0.57735; Weight Factor = 1 ........................ 116 6.7 Division of Flow Domain Into Finite Elements ............. 120 6.8 Comparison of Solutions for Flow from an Injection Well in a Bounded Leaky Aquifer (Constant Head Boundary, s = 0 at r = 200m) ......................... 121 6.9 Flow Domain for Numerical Location of the Interface in a Leaky Aquifer ........................................ 123 6.10 Comparison of Solutions for a Steady Lens in a Leaky Aquifer ................................................ 125 6.11 Comparison of Solutions for a Steady Lens in a Leaky Aquifer ....................................... .. .... 126 7.1 Finite Difference Grid Used With the INTERCOMP Model ..... 132 7.2 Comparison of the Concentration Distributions Predicted by the INTERCOMP Model and the Approximate Solution of Hoopes and Harleman ........................ 136 7.3 Lines of Constant Concentration After 40 Days of Injection ........................................... 137 7.4 Lines of Constant Concentration After 40 Days of Injection ............................................. 139 viii LIST OF FIGUREScontinued Figure Page 7.5 Lines of Constant Concentration After 40 Days of Injection ......................................... 140 7.6 Movement of the 50% Concentration Line With Continuous Injection ................................. 142 7.7 A Dimensionless Representation of Gravitational Segregation During Lens Growth ........................ 144 7.8 Comparison of Lens Boundaries After 30 Days of Injection With Different Dispersivities ........................ 145 7.9 Movement of 50% Concentration Lines During Production of Freshwater From a Lens ........................... 147 7.10 Movement of 50% Concentration Line During Production of Freshwater From a Lens ........................... 148 LIST OF SYMBOLS The following symbols are employed in this dissertation. Dimensions, where applicable, are also given in terms of length L, mass M, time T in square brackets. A = Area [L2] a = depth of the top of the confined aquifer below sea level, [L] B = leakage coefficient, [L] b = thickness of a confined aquifer [L] b' = thickness of a leaky layer, [L] bx = body force in the xdirection, [ML/T2 b = body force in the ydirection, [ML/T2] bz = body force in the zdirection, [ML/T2] c = fractional concentration of salinity D = the dispersion tensor, [L2/T] D = coefficient of longitudinal dispersion, [L2/T] Dm = coefficient of molecular diffusion, [L2/T] Dt = coefficient of transverse dispersion, [L2/T] d = spacing between plates of the HeleShaw model, [L] F = an arbitrary function f. = the jth term of the forcing matrix, [L3/T] g = acceleration due to gravity, [L/T ] h = vertical thickness of the flow region, [L] I0 = modified zeroth order Bessel function of the first kind Il = modified first order Bessel function of the first kind LIST OF SYMBOLScontinued i = unit vector in the xdirection [J] = the Jacobian matrix Jx = xdirection component of head gradient J = ydirection component of head gradient Jz = zdirection component of head gradient j = unit vector in the ydirection K = hydraulic conductivity for isotropic porous media, [L/T] K' = hydraulic conductivity of a leaky layer, [L] KO = modified zeroth order Bessel function of the second kind K0 = modified first order Bessel function of the second kind K.. = components of the hydraulic conductivity tensor for J (i, j = x, y, z), [L/T] Ki. = terms of the global stiffness matrix (in Chapter 6), [L2/T] k = intrinsic permeability, [L2] A k = unit vector in the zdirection Kij = terms of the element stiffness matrix (in Chapter 6), [L2/T] L = a linear operator Nth Ni = the it shape function n = porosity of the porous medium P = pressure [M/(LT2)] Pf = pressure on the freshwater side of an interface, [M/(LT2] Ps = pressure on the saltwater side of an interface, [M/(LT2)] Q = volumetric rate of injection into a well, [L3/T] Q' = dimensionless injection rate Qx = uniform horizontal flow per unit width of aquifer, [L2/T] q = discharge velocity vector, [L/T] LIST OF SYMBOLScontinued qt = longitudinal component of discharge velocity, [L/T] qx = xdirection component of discharge velocity, [L/T] qy = ydirection component of discharge velocity, [L/T] qz = zdirection component of discharge velocity, [L/T] r = the radial coordinate of a cylindrical coordinate system, [L] re = radius of a constant head boundary, [L] rt = radius of the toe of the interface, [L] S = the boundary surface, [L2] Ss = specific storage, [1/L] s = drawdown or pushup, [L] t = time, [L] t' = dimensionless time u = xdirection component of flow velocity, [L/T] V = volume, [L3] v = ydirection component of flow velocity, [L/T] Wn = weighting factor for the nth sample point w = zdirection component of flow velocity, [L/T] ws = flow rate through the leaky layer, [L/T] X = projection of the interface on the horizontal, [L] Xi = x coordinate of node i, [L] x = a horizontal coordinate of a cartesian coordinate system, [L] x' = a dimensionless distance xt = distance to the toe of the interface, [L] Xw = distance from the injection well to the sea coast, [L] x' = dimensionless position of the injection well S= coordinate of node i, L] Yi = y coordinate of node i, EL] LIST OF SYMBOLScontinued y = a horizontal coordinate of a cartesian coordinate system, [L] y' = a dimensionless distance Z = z coordinate of node i, [L] z = elevation above a datum, or the vertical coordinate, [L] z. = elevation of a point on the interface, [L] 1 a = coefficient of aquifer compressibility, [LT2/M] a = a coefficient (in Chapter 5 and 6) at = longitudinal dispersivity, [L] at = transverse dispersivity, [L] 6 = coefficient of fluid compressibility, [LT2/M] S= a constant (in Chapter 5), [L] 2 = a modified leakage coefficient (in Chapter 6) [T] y = unit weight of a fluid, [M/(LT)2] 6= the GhybenHerzberg coefficient C = the third dimensionless coordinate in an oblique coordinate system n = the second dimensionless coordinate in an oblique coordinate system 0 = angular coordinate in a cylindrical coordinate system S= dynamic viscosity, [M/(LT)] v = kinematic viscosity, [L2/T] S= the first dimensionless coordinate in an oblique coordinate system 3 p = fluid density, [M/L3] p = average density, [M/L ] D = the Girinskii potential (in Chapter 5), [L3/T] S= nodal values of for n = 1 .20, [L] xiii LIST OF SYMBOLScontinued = piezometric head, [L] = an approximation to the piezometric head, [L] (0 = undisturbed piezometric head in a static aquifer, [L] f = piezometric head of freshwater measured in terms of fresh water, [L] ps = piezometric head of saltwater measured in terms of saltwater, [L] Subscripts f  i  m  p  r  s  t  x  y  z  and Superscripts freshwater or lighter fluid on the interface model prototype ratio saltwater or heavier fluid at the toe of the interface xdirection ydirection zdirection Abstract of Dissertation Presented to the Graduate Council of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy ANALYSIS OF FRESHWATER LENS FORMATION IN SALINE AQUIFERS By John Patrick Glass June, 1979 Chairman: B. A. Christensen Major Department: Civil Engineering The injection of freshwater into a saline aquifer can, under favorable conditions, result in displacement of the saltwater and crea tion of a distinct body of freshwater called a freshwater lens. It may be desirable to do this either as a method of storing the freshwater for future use or in order to form a barrier against saltwater intrusion. In either case, it is important to know beforehand whether the lens will be stable or the injected water will be lost. The purpose of this study is to develop methods for analyzing the hydrodynamic factors involved in the formation and maintenance of such lenses. A series of laboratory tests using vertical HeleShaw models was performed to determine the shape and rate of movement of the initially vertical interface between two fluids of different density. The results of these tests were compared with the predictions of two formulas proposed by previous investigators. It was found that for relatively thick aquifers it is not correct to neglect the vertical components of flow as has formerly been done. An analytical method is presented for the location of the boundary of a steady freshwater lens in a coastal aquifer. A vertically integrated potential function is used to superimpose the flow from the injected well on the preexisting uniform seaward flow in the aquifer. This results in a quasithree dimensional solution for the position of the coastal interface. A pair of coupled ordinary differential equations is derived which describe the steady flow of freshwater injected into a leaky saline aquifer under the Dupuit assumption. The nonlinear system of equations is solved numerically by isolating the unique integral curve which results in the only physically possible boundary conditions. A design chart is presented by which the maximum possible size of freshwater lens can be determined for various aquifer conditions. A three dimensional finite element model is developed which can be used to locate the boundaries of a steady freshwater lens under a wide variety of conditions. The model uses isoparametric elements of the quadrilateral family having twenty nodes each. This permits very close approximation to the curved flow boundaries typical of freshwater lenses. Inhomogeneous and anisotropic aquifers can be treated. A finite difference model developed for the U.S. Geological Survey is used to analyze the axisymmetric flow in transient freshwater lenses. Simulation of lens formation in a variety of saline aquifer configura tions indicates that previously available methods of analysis which neglect vertical components of flow may be reliable in very thin aquifers. For thicker aquifers, however, vertical flow can predominate resulting in a much more rapid rate of gravitational segregation of the fluids. xvii CHAPTER 1 INTRODUCTION Groundwater flow is an important part of the hydrologic cycle which nature uses to even out the temporal and spatial variations inherent in rainfall. It has long been man's practice to take water from groundwater reservoirs with the assumption that the water removed will be replaced by natural processes. In many cases, this arrangement works well, with nature putting water into the ground and man taking it out. In some areas, though, the increased use of groundwater supplies coupled with decreased natural recharge due to manmade changes in the land surface is either draining the freshwater aquifers or opening them to saltwater intrusion. In recognition of this, efforts have been made to augment natural recharge by such means as spray irrigation, water spreading, seepage ponds and deep well injection. This signals a new conception of aquifers as a form of underground storage which can be used in much the same way as surface reservoirs. As an extension of this, it has been proposed to make use of the storage capacity available in saline aquifers by injecting freshwater into them in such a way that the resident salt water is displaced and a lens of freshwater is created. The variety of possible flow patterns that may result from the injection of freshwater into a saline aquifer can be viewed as a contin uous spectrum bounded by three idealized extremes. One possibility is that the freshwater will disperse rapidly and lose its identity as a separate flow, merely serving to dilute the resident saltwater. In 2 another case, a flow pattern may develop which transports the freshwater away from the point of injection as a plume. A third alternative is that the freshwater will tend to remain near the point of injection creating an identifiable lens of freshwater which grows or diminishes accordingly as fluid is added to or removed from it. If the freshwater is injected in anticipation of its subsequent retrieval and reuse, the formation and maintenance of such a lens adjacent to the point of injection is essential. The major factors influencing the formation of freshwater lenses are: (1) buoyancy, (2) mode of injection, (3) the properties of the porous medium, (4) dispersion, (5) flow of the resident saltwater. The relative intensities of these'individual influences determine which of the previously described flow patterns will predominate. Favorable conditions for the development of freshwater lenses require that buoyancy, dispersion, and preexisting flow be relatively weak influences. In previous studies of freshwater lenses such ideal conditions have been assumed. In reality, these effects are never entirely absent. It is important, then, to determine at what point their effects become substantial and to have techniques for including their influences in the analysis of freshwater lenses. CHAPTER 2 LITERATURE SURVEY General The analysis of freshwater lenses is one of a general class of problems involving the stratified flow of two miscible liquids through porous media. It is closely related to such other problems of that class as saltwater intrusion, interface upcoming beneath wells, and deep well disposal of liquid wastes. A wealth of published information is available on each of these topics but much of it has no direct application to the study of freshwater lenses. Many of the methods developed for locating coastal interfaces, for instance, depend on conformal mapping techniques or otherwise on the theory of complex variables, which is limited to two dimensional problems. Likewise, much of the work done with deep well disposal deals with the movement and dispersion of plumes, emphasizing the rapid departure of the waste fluid from the scene of injection. Nonetheless, some of the work in each of these fields can be useful in lens analysis. It is only these studies and those which are directly concerned with fresh or hot water lenses that will be mentioned here. Studies Dealing with Interface Location The rotational movement of an initially vertical interface between two fluids of different density was studied byGardner, Downie and Kendall (1962) because of its applicability to recovery processes in petroleum reservoirs. The initial motion of the interface and the long term motion were found analytically using the assumption that the porous medium was thin and vertical flow could be neglected. These two solutions were then joined to give one equation which was assumed accurate for intermediate times as well. It was assumed in these derivations that the interface remained straight during its rotation toward the horizontal. Laboratory tests using models of square and circular crosssection packed with glass beads were performed which seemed to agree with the equation. Dagan and Bear (1967) treated the upward movement of a freshwater saltwater interface beneath a partially penetrating well using a pertur bation expansion to linearize the interface boundary condition. They checked their results with a sandpacked model and found them to be reasonable when the disturbance of the interface was small. The range of applicability of this method, however, appears to be quite limited. Shamir and Dagan (1971) developed a transient numerical model to predict the landward movement of a saltwater wedge in a coastal aquifer. The Dupuit assumption was used to formulate the governing equation which was then solved by an implicit finite difference method. Dispersion was neglected. The numerical results were checked against an analytical solution similar to that of Gardner, Downie and Kendall and with a HeleShaw model. The comparison with the analytical solution showed good agreement that the interface moved as a straight line. Comparison with the HeleShaw model, however, was not quite as good. The authors attributed this to the fact that vertical velocities were neglected in both the analytical and numerical solutions. Strack (1976) studied the influence of inland production wells on the steady state shape of the saltwater wedge in coastal aquifers. His analysis makes use of a vertically integrated potential function which can be defined so as to be singlevalued, throughout the aquifer. Lines of constant potential can then be found in the horizontal plane which lead to a threedimensional description of the interface based on the GhybenHerzberg principle. A similar problem was studied by Kishi and Fukuo (1977). They considered the effect of inland pumping on the saltwater wedge in a fan shaped coastal aquifer. The Laplace equation in radial coordinates is integrated vertically to make it two dimensional. Using the Ghyben Herzberg principle the dependent variable is transformed so that it will be singlevalued in both the freshwater region and the interface region. The resulting equation is solved by the method of Green's functions to give a threedimensional description of the interface. Studies Dealing with Dispersion Gardner, Downie and Wyllie (1962) considered the mixing across the interface between methane gas and nitrogen used as a cushion gas in a natural gas reservoir. A solution to the convective dispersion equation in radial coordinates was combined with the gravity segregation formula found by Gardner, Downie and Kendall. Harleman and Rumer (1962) derived an equation for interface rota tion which is the same as that found by Gardner, Downie and Kendall for the latter stages of movement. This equation also is based on the assumption of a straight interface. The same report also deals with longitudinal and lateral dispersion in both stationary and moving coastal interfaces. Laboratory tests using sandpacked models to simulate a coastal aquifer showed that in the absence of tidal fluctua tions the mixing zone was quite thin, especially near the toe of the interface. An equation was presented which predicted the thickness of the dispersion zone accurately in this region. When changing conditions in either the saltwater or the freshwater caused a shift in the equi librium position of the interface the effects of longitudinal dispersion became pronounced. It was concluded that such shifts are responsible for most of the dispersion found in natural coastal aquifers. Wooding (1963) obtained approximate solutions for the lateral dispersion across a stationary interface by transforming the governing equations into a natural coordinate system based on the shape of the interface. Using order of magnitude reasoning these equations were reduced to the form of Prandtl's boundary layer equations. For these, the usual similarity solution was obtained which relates the thickness of the mixing layer to the distance from the stagnation point and the slope of the interface. Hoopes and Harleman (1965 and 1967) analyzed the dispersion of tracers injected into groundwater. Their work consisted of analytical and numerical solutions to the convective dispersion equation and laboratory experiments with a sandpacked model to verify these solu tions. The relative effects of convection, lateral and longitudinal dispersion, and molecular diffusion were compared for the case of axisymmetric radial flow from an injection well and for the flow between an injection well and production well. It was shown that the effects of lateral dispersion on the steady flow between two wells was almost negligible except near the line joining the two wells. Studies Concerned Directly with Freshwater Lenses Glazunov (1967) proposed to use a cluster of three to five wells operated in various patterns of injection and production in order to form freshwater lenses of the desired shape. His analysis neglects both dispersion and gravity segregation, treating the convection of an interface in two horizontal directions by a method introduced by Muskat (1946). In this method the interface is represented as a material surface which is transported along the streamlines by the flow while exerting no effect on it as a boundary. Esmail and Kimbler (1967) studied the effects of radial flow geometry on gravity segregation and longitudinal dispersion in labora tory models composed of glued sand grains. The solution to the convec tive dispersion equation given by Gardner, Downie and Wyllie was verified for cyclic storage in freshwater lenses for several cycles of injection and production in the model. A dimensionless parameter analysis of the influence of radial geometry and interfacial mixing on the rate of gravity segregation was conducted resulting in an empirical equation for the rate of interface rotation. It was found that there is inter action between the processes of dispersion and interface rotation which retards gravity segregation. The authors concluded that the horizontal projection of the interface is proportional to the square root of the density gradient across it and that gravity segregation in their model had more effect on the recovery efficiency than did dispersion. Under these conditions, interface dispersion may actually improve the recovery efficiency by retarding gravity segregation. Kumar and Kimbler (1970) continued the research with glued sand models, finding that the recovery efficiency improves after several cycles of injection and recovery. Whitehead (1974) developed a computational method for predicting the recovery efficiency of a single well cyclic storage system in thin horizontal saline aquifers. His method was based on the empirical equation of Esmail and Kimbler. Whitehead verified the results of his computational scheme using the glued sand grain model and found them to be accurate to within 10% for three storage cycles. The recovery efficiencies achieved in the laboratory tests were in the range of 70100%. It was concluded that, given favorable subsurface conditions, the use of saline aquifers for freshwater storage could be economically feasible. In a series of reports (Agrawal, 1975; D'Amico, 1975; Tate, 1976; Kimbler and Whitehead, 1977) the results of Whitehead were extended to include the effects of aquifer dip and viscosity difference between the resident and injected fluids. These studies also used glued sand models. The conclusions reached from these experiments were that the best conditions for cyclic storage occur when the saline aquifer is horizontal and the injected fluid has the same viscosity as the resident fluid. It is widely recognized that a preexisting uniform flow in the saline aquifer is detrimental to the recovery efficiency of a lens storage system. Langhtee (1974) proposed to surround the storage area with bounding wells which would be pumped in such a way that a stag nation zone would result inside the boundary. The required pumping rates for these wells were found by superimposing the streamline and equipotential pattern for the individual wells and solving for the pumping rates which would give the best leastsquares fit to the desired piezometric head on the boundary. Molz and Bell (1977) also considered the use of boundary wells for head gradient control. They used a linear programming technique to find the optimum set of bounding wells which would give the desired stagna tion in the storage area with the least expenditure of energy for pumping. The shape and movement of naturally occurring freshwater lenses on Grand Cayman Island were studied by Childley and Lloyd (1977). They used a modified GhybenHerzberg ratio based on field measurements of salinity distribution in the mixing zone of the lens. The long term transient effects of seasonal variations in recharge rate were modeled numerically using the idea of a succession of steady states. Comparison of the numerical results with records of actual groundwater conditions on the island indicated that their model was useful in predicting long term changes due to water supply wells and changes in recharge condi tions. They concluded that, while it may take many years for a pumped lens to reach equilibrium, an overpumped lens can collapse locally in a matter of a few months. Smith and Hanor (1975) report on a field test of the freshwater lens storage concept. A 10 centimeter fully penetrating well was drilled into a confined saline aquifer 36 meters thick to be used for injection and retrieval of freshwater. After the construction of this well it was discovered that a uniform flow of approximately 15 centi meters per day existed in the aquifer. Recognizing the unfavorable conditions, the investigators proceeded with two tests, each using about 757,000 liters of freshwater. In the first test, the water was stored for two hours and 41% was recovered. In the second test the storage period was six days and the recovery efficiency was 24.5%. Werner and Kley (1977) conducted a field test of the storage of heated water in a water table aquifer 3.35 meters thick. A total of 430 cubic meters of water at a temperature of 45C was injected and the movement of the heat was monitored with buried temperature sensors for a period of 64 days. The heated water remained close to the injection well throughout the test period but no attempt was made to recover it. Molz, Warman and Jones (1978) reported a field test with a partially penetrating well in a confined aquifer 21 meters thick. A group of 14 observation wells were constructed to monitor the movement of the hot water injected. In this test, storage of 7570 cubic meters of hot water for 36 days resulted in an energy recovery factor of 0.69 which was considered promising. Papadopulos and Larson (1978) applied a finite difference numerical model developed for the U.S. Geological Survey (INTERCOMP, 1976) to simulate the test results of Molz, Warman and Jones. The predicted recovery efficiency was 75% versus a value of 69% measured in the field. Considering the accuracy of the measured aquifer properties, this prediction was considered quite good. CHAPTER 3 BASIC EQUATIONS AND FLUID PROPERTIES Basic Equations Darcy's Law The groundwater flow to be dealt with in this report is limited to laminar flow in porous media with only primary porosity. The equation of motion for this flow is Darcy's law which can be expressed, in the cartesian coordinate system, as qx = Kxx x + Kxy Jy + Kxz J q = x J + Ky + K J qz = Kzx J + Ky Jy + Kz J q X9 qy, q (3.la) (3.1b) (3.1c) = components of the discharge velocity in the x, y, and z directions, respectively = components of the hydraulic conductivity tensor for (i, j = x, y, z) ax ay where z z zz S = piezometric head If it is assumed that the coordinate axes are parallel to the principal directions of hydraulic conductivity then Darcy's law can be simplified to = K K j K l k (3.2) xx ax yy ay zz a where i, j, and k are unit vectors in the x, y, and z directions, respectively. This assumption is not a major restriction to the use of equation (3.2) because in practice it is usually difficult to distin guish any directional hydraulic conductivities except the vertical and the horizontal. If the soil is isotropic or only one dimensional flow is being considered equation (3.2) reduces to q =K V) (3.3) The value of the hydraulic conductivity, K, or of its directional components depends on the properties of both the soil and the fluid. Since the properties of saltwater and freshwater are somewhat different it may be necessary in the general case to consider the hydraulic conductivity as a function of the salt concentration. The flow resis tance of the soil is then described in terms of its intrinsic permea bility, k, which is related to hydraulic conductivity by K = kg/v (3.4) where v is the kinematic viscosity of the fluid and g is the accelera tion due to gravity. The piezometric head, q, appearing in equations (3.2) and (3.3) is defined as S= z + P/y (3.5) which depends on the unit weight, y, of the fluid concerned. Since y is also a function of salt concentration it is often more convenient to deal directly with pressure, P, in the equation of motion. Continuity The continuity equation for groundwater flow is V q = (3.6) in which Ss, the specific storage, is defined as Ss = pg(a + nB) (3.7) where a = coefficient of soil compressibility B = coefficient of compressibility for water n = porosity p = fluid density The combination of equations (3.2) and (3.6) gives the governing equation for groundwater flow Kxx + K I + K S (3.8) ax xx x ay 3yy y az zza z s t .x j i /Z[ S 38 For steady flow in homogeneous and isotropic soil this reduces to the Laplace equation += 0 (3.9) ax 2 y az Well Flow in Leaky Aquifers In dealing with flow to and from wells, it is customary to denote the head gradient in terms of drawdown, s, which is defined as s = ( Qo (3.10) where (0 is the head above the leaky layer (see Figure 3.1). The flow from an injection well causes an increase of head, in which case s will be called "push up" instead of drawdown. A leaky aquifer is one which is confined by a layer of material which has low hydraulic conductivity, but not so low that it can be considered impervious. The governing equation for flow in a leaky aquifer is V2s s/B2 s as (3.11) Kb at The leakage coefficient, B, is given by B2 = Kbb'/K' (3.12) where K' = hydraulic conductivity of the leaky layer b' = thickness of the leaky layer b = thickness of the aquifer The Dispersion Equation The mixing of saltwater with freshwater at the periphery of a freshwater lens is described by the convective dispersion equation ac 1 +t nq Vc = V D Vc (3.13) in which c is the concentration of the fluid being followed and D is the dispersion tensor. For isotropic porous media this tensor has only two independent terms, D. and Dt, corresponding to longitudinal and transverse dispersion respectively, relative to the direction of the mean flow. The values of these coefficients depend not only on the properties of the porous medium but also on the flow velocity and the properties of the fluid involved. It is often assumed that these coefficients can be represented by = Dm + a q /n (3.14a) Dt = Dm + t qq/n (3.14b) where Dm = coefficient of molecular diffusion in the porous medium aC = longitudinal dispersivity at = transverse dispersivity qp = discharge velocity in the direction of mean flow The dispersivity coefficients a, and at are properties of the porous medium only. Their magnitudes are related to the intrinsic permeability and the pore size distribution, but the actual values must be found experimentally. The longitudinal dispersion coefficient, D,, is usually 10 to 20 times greater than Dt but the ratio of the two depends on the flow velocity, being close to unity when the flow is very slow and the influence of molecular diffusion predominates. For axisymmetric radial flow of a tracerlike substance from an injection well, equation (3.13) can be simplified to ac ac nr 3c c = Q ) a c (3.15) it 2nbnr ar X 2abnr @2 V because there are no components of either velocity or concentration gradient except in the radial direction. This equation has been solved numerically for the case of steady flow from the well and con tinuous tracer injection beginning instantanously (Hoopes and Harleman, 1965; and Shamir and Harleman, 1966). Freshwater injected into a saline aquifer does not behave as a tracer because of the difference in density, so equation (3.15) does not apply even though the flow may be asixymmetric. Because vertical components of both velocity and concentration gradient are present, equation (3.13) must be used. The numerical solution of this equation will be the subject of Chapter 7 of this report. Boundary Conditions at an Interface When the mixing zone separating saltwater and freshwater is relatively thin, it is often convenient to neglect dispersion and assume that the two fluids are separated by a sharp interface. If such a boundary between the two fluids exists, a separate governing equation can be written for each zone and the two flows will be coupled by their common boundary condition at the interface. The interface is a material surface which is always composed of the same fluid particles. The essential features of such a surface are that both pressure and the velocity component normal to the interface must be continuous across it. For a stationary interface the normal velocity component is zero so the boundary condition is concerned with pressure only. The pressure on the freshwater side of the interface is Pf = Yf (f zi) (3.16a) and on the saltwater side PS = YS (s zi) (3.16b) where zi is the elevation of a point on the interface and the sub scripts f and s stand for freshwater and saltwater, respectively. Requiring that these pressures be equal results in an expression for the boundary condition at a stationary interface Ys Yf z. Y Y f (3.17) i AY s AY f where A = Ys Yf For the special case where the saltwater is at rest and a hydro static pressure distribution prevails in the freshwater zone this reduces to the GhybenHerzberg relationship Yf Az = AAf (3.18a) Ay 19 or Az = 6APf (3.18b) where 6 = Yf/AY is the GhybenHerzberg ratio. This states that a change in the piezometric head of the freshwater causes a change in the eleva tion of the interface that is opposite in sign and greater in magnitude by a factor of 6. If the interface is not stationary the condition of continuity of the normal velocity component must also be used. The interfacial surface can then be described by some function, F, such that F (x, y, z, t) = 0 (3.19) This surface is convected by the flow so that F+ q vF = 0 (3.20) at n where q = discharge velocity of interfacial particles n = porosity of the aquifer In terms of the velocities adjacent to the interface on both sides of it this becomes n  + q VF = 0 (3.21a) for the freshwater, and aF 0 n F + q VF = 0 (3.21b) at s for the saltwater. In terms of the elevation of a point on the inter face, zi, the function, F, can be written as F(x, y, z, t) = z zi (x, y, t) = 0 (3.22) Substituting this, together with Darcy's law, into equations (3.21) gives az. n KVqf v(z z.) = 0 (3.23a) at 3Z. n Kv v(z z ) = 0 (3.23b) at 1 Performing the vector operations indicated in these equations results in the boundary conditions for the freshwater zone and the saltwater zone. The interface boundary condition for the freshwater zone, then, is Yf a f Ys 4s Yf 2 Y nn K (V~ f) + K (f s) AY at y at Ay f Ay f K 0 (3.24a) az and for the saltwater zone, Yf 4, Ys s 2 Yf nf f s + K (s ) + K K (vqf ) AY at Ay at a S Ay f K s = 0 (3.24b) az Fluid Properties Density Up to this point the fluids with which this report is concerned have been referred to simply as freshwater and saltwater. These conven ient terms are used only in a relative sense and the actual properties of the fluids deserve some further elaboration. In the ensuing analy sis the primary difference between the resident and injected fluids will be their density. In general, this is a property that depends on pressure, temperature and salt concentration. In this report, however, the change in density due to compressibility will be neglected because it is so small, and the thermal effects, while not necessarily small, will be eliminated by assuming that the temperature of the injected fluid is the same as that of the resident saltwater. Figure 3.2 shows the relationship of density to temperature for water with several different salt concentrations. With the exception of the 3.5% concentration line, this figure is based on data for solutions of pure sodium chloride in water. The 3.5% concentration line shows the temperature versus density relation for seawater having a dissolved solids concentration of 3.5%. Saline groundwater is commonly associated with the intrusion of seawater into coastal aquifers. However, it is also found in many locations which have no apparent geological association with the oceans and where its concentration may be much greater or less than that of seawater. Figure 3.3 shows the relationship of density to sodium chloride concentration. It is worthwhile to note that in this range the relationship is apparently linear. 1.04 1.03 3. 5 1.02 E S 1.01 I'" =7% 1.00 0% 0.99 0.98 0 10 20 30 40 TEMPERATURE (DEGREES C.) FIGURE 3.2. TEMPERATURE VS. DENSITY FOR WATER WITH VARIOUS SALT CONCENTRATIONS DP C? 23 U1 r , c:a ru cr ___ ___ __  HL LOU ~  r)I j    i' UU _,W/6)____SN30 a) LU ` ~  ~ ~ LU (~w/6.l) __ AIISNO D LU  I~IL. NT~j " ` Viscosity It was shown in equation (3.4) that the hydraulic conductivity of an aquifer is partially determined by the kinematic viscosity of its water. Viscosity is a function of both temperature and salt concentra tion but within the range of present interest temperature is the more important factor of the two. The temperature dependence of kinematic viscosity for both pure water and seawater is shown in Figure 3.4. It is interesting to note that the kinematic viscosity of pure sodium chloride solutions, as given by Kaufmann (1960) for concentrations up to 4%, differ so little from the kinematic viscosity of pure water that they could not be distinguished from it if plotted at the scale of Figure 3.4. Since differences in temperature are not being considered in this study and the influence of salt concentration on viscosity is relatively weak, viscosity differences will be neglected also. 1.4 1.3 1.2 in o I > FRESHWATER 0.8 0.7 S10 20 30 4 TEMPERATURE (DEGREES C.) FIGURE 3.4. VISCOSITY VS. TEMPERATURE FOR 0 gg FRESHWATER 0.8 0.7 0 10 20 30 40 TEMPERATURE (DEGREES C.) FIGURE 3.4. VISCOSITY VS. TEMPERATURE FOR FRESHWATER AND SALTWATER CHAPTER 4 THE HELESHAW MODEL General Valuable insights into the flow phenomena associated with the injection of fresh water into saline groundwater can be obtained by studying certain simplified aspects of the problem. This chapter will describe the use of a physical model to analyze the movement of the interface between the two fluids in the absence of dispersion. The HeleShaw analog was chosen for this work because it is simple and inexpensive to construct and the test results are easily observable without elaborate instrumentation. The HeleShaw analog consists of one or more viscous liquids flowing in the gap between two closely spaced flat plates. Its usefulness in groundwater studies is based on the fact that the NavierStokes equations governing the viscous flow in this gap can be reduced to a form similar to Darcy's law in two dimensions. The Viscous Flow Analogy If u, v, and w are the velocity components in the x, y, and z directions, respectively then the NavierStokes equations for incom pressible flow in a cartesian coordinate system are u + u + v u w b + vV 2u (4.1a) at x + y+ az x p ax + u @v + v + w b + vv (4.1b) at ax ay az y P y aw aw cw aw ) aP 2 w+ + v + w W b P + v2w (4.1c) t ax 3y z z p Tz where bx, b and bz = components of acceleration due to body forces in the x, y, and z directions, respectively p = density of the fluid v = kinematic viscosity of the fluid P = pressure Applying these equations to the flow between the plates shown in Figure 4.1, it can be assumed that because of the low Reynolds number (Re < 1 for this model) the flow is laminar and there is no velocity component, v, in the direction normal to the plates. Note also that, since the velocity at the surface of the plates is zero, the velocity gradients normal to the plates are very large compared to those in the x and z directions so that the latter can be neglected. Finally, the only body force acting on the fluid is gravity, so bz = g and bx = b = 0. Incorporating these observations into the NavierStokes x y equations and multiplying by p results in aP a2u =x 2 (4.2a) ay DP 0 (4.2b) .P 82u Y + z 2U (4.2c) Y 3y 2 "y where y = pg is the unit weight of the fluid and p = pv is its dynamic viscosity. Recalling that the piezometric head, 4, is defined as 28 x U1 ui i ui 1= a SI 0 /01 Iz: LU U '=4 N __ / 3 >1w { = z + P/y (4.3) it follows that equations (4.2) when expressed in terms of q and divided by y become ac p a2u (4.4a) ax y 2 S= 0 (4.4b) ay = a2u (4.4c) az y 2 From the second of these equations it is seen that at any point (x, z) in the plane of the model, the piezometric head is not a function of y. Therefore, the first and third of these equations can be integrated with respect to y giving y au + C (4.5a) ax Y ay y = 2w + C (4.5b) az y ay The constant of integration, C, is found to be zero by observing that at the center of the flow region, where y = 0 both and are also zero, so that y = u (4.6a) ax y ay y = w (4.6b) az y Dy Integrating these equations again with respect to y gives Y i U + C1 2 Dx Y 1 2 2 az 2 az (4.7a) (4.7b) p w + C1 Y 1 Here, the constant of integration, C1, is found by requiring that, at the surface of the plates, where y = d/2, u = w = 0. With these boundary conditions the expressions for the velocity distribution across the gap are found to be u = (y2 d2) 2 S(y 2 w = (y2 d2 ) 2p 4 az where d is the thickness of the to find a velocity analogous to through porous media, equations thickness of the gap as follows q= d/2 d/2 d/2 d/2 d/2 gap between the plates. Now, in order Darcy's discharge velocity for flow (4.8) must be averaged across the dy d Y 4 y 12 P ax dy = 12 12 P az (4.9a) (4.9b) If the constant coefficients in equations (4.9) are used to define an effective hydraulic conductivity, Km, for the model Km 12 (4.8a) (4.8b) then the governing equations for flow between the plates become q Km a S m ax qz = Km  qz m Dz in which the resemblance to the obvious. In further analogy to ability for the HeleShaw model, so that two dimensional form of Darcy's law is groundwater flow, an intrinsic permea km, can be defined as km = d2/12 Km = km (4.12) (4.13) The continuity equation for the flow between the plates is au + w _ ax az (4.14) since there is no velocity in the y direction. Averaging these veloc ities across the gap, as before, results in (4.15) ax (Kmax + (Km )= 0 which has the same form as the governing equation for two dimensional flow in an incompressible porous medium. The model hydraulic conductiv ity, Km, is a constant depending only on the spacing, d, of the plates and the unit weight and viscosity of the fluid. Consequently, for parallel plates and homogeneous fluids the flow is homogeneous and isotropic, so that equation (4.15) reduces to the Laplace equation in (4.11a) (4.11b) two dimensions S + 2 = 0 (4.16) ax az2 Since Km is a function of physical parameters which are relatively easy to measure and control, the HeleShaw model is a very convenient tool for the scale modelling of two dimensional groundwater flow. Model Scaling In the preceding section it was shown that two dimensional ground water flow can be modeled by the HeleShaw analog because the governing equations for the two cases are of the same form. If quantitative results are to be obtained from the model, however, the scaling ratios between model and prototype must be chosen so that the governing equations are actually identical. In order to achieve this, the model laws which fit the HeleShaw model to the prototype case of groundwater flow must be developed. Since the present case involves the modeling of two immiscible liquids separated by an interface the equations on which the model laws depend are Darcy's law, the continuity equation and the boundary condition for interface flow. Letting the subscript m denote model parameters and p denote those in the prototype, Darcy's law for the model is i q m= Km (4.17a) 1 m a7b) qzm Km (4.17b) where the superscript : I : for the lighter fluid for the heavier fluid The corresponding equations for the prototype are q = xp ai1 q = K  zp p azp Let the scaling ratios between be denoted by the subscript r, model parameters and prototype parameters so that Kr = Km/Kp (4.19a) (4.19b) xr = x/xp (4.19c) etc. Then, substituting these relationships into equations (4.18) gives K X41 S i Km r m q xm/qxr K ri X qm m m zm zr Kr 3azr r r r If these are to be identical to equations (4.17) then it must be ture that Ki qi r xr Xr (4.20a) (4.20b) (4.21a) K  p p (4.18a) (4.18b) K 1 qr rr (4.21b) zr zr Equations (4.21) provide one of the model laws which allow results from the model to be transferred to the prototype. The continuity equations for isotropic flow in the model and prototype are 2i 2i + =0 (4.22) Xm m and 2 i 2i S+ = 0 (4.23) ax 2z p p Substituting the appropriate parameter ratios into equation (4.23) and requiring that it be identical to equation (4.22) results in xr = zr (4.24) The third similitude requirement is that the boundary conditions at the moving interface be given by identical equations in both model and prototype. These interface equations were given in Chapter 3. For the model they are Sf s s f 2 s n m mm n m m K (V ) + K (Vm .Vs)  mm Y t m y t m Ay m m Am m mn am Ym m m m m m = 0 (4.25a) K zm and f f s s s 2 f m m Ym mm Ym m f s n m n ( K Vm) + K (Vm Vs)  Aym atm m ay t m Ay m s Kmiz 1 =0 (4.25b) The corresponding equations for the prototype are the same as equations (4.25) except that the m subscripts are replaced by p. When the param eter ratios are substituted into the prototype equations and equality is required with the model equations, 16 groupingsof parameter ratios are generated, which must all be equal to unity. Many of these group ings are repetitive and, since they provide no new information, they will not be listed here. The groupings which are useful are 1 AY, t 1 AY t 1 Ay xr2 1 z 1 z r r r r r r r r f2 s (4.26) nr Yr r nr Yr r K Yr r Kr r Kr r The equality of the first and third terms of this equation gives the time scale for the model nr xr tr r (4.27) K 4 Equality of the first, second, and fourth terms of equation (4.26) together with equation (4.27) results in ff s A yr = r r = AYr zr (4.28) Finally, the equality between the last two groupings in equation (4.26) requires that f s =s r (4.29) which, in turn, means that f s Y = Y, (4.30) The complete list of similitude requirements for matching the prototype to the model is given in Table (4.1). Model Construction The design of the HeleShaw model involved the following con sideration: 1. The plates must be transparent, relatively inflexible and demonstrably flat. 2. Provision must be made for the controlled introduction of the viscous fluids and the release of entrapped air. 3. The model must be easily disassembled for cleaning between tests. 4. It must be possible to rotate the model from horizontal to vertical so that the initial shape of the interface can be controlled. The availability of materials and the requirement for transparency dictated that the material to be used for the flat plates be either 1/4 inch plate glass or 1/2 plexiglass. Since 1/2 plexiglass is more than 2.5 times more flexible than 1/4 inch plate glass, it was decided to use the latter. TABLE 4.1 SUMMARY OF SIMILITUDE REQUIREMENTS FOR VERTICAL HELESHAW MODEL WITH TWO IMMISCIBLE FLUIDS Velocity Ratios: i r ri q  xr Xr fh where i = S, and 1 Kz r q zr r for the lighter fluid for the heavier fluid Geometric Ratios: x =z r r Time Ratio: 2 nr Xr t r f Kr r Head and Density Ratios: ff s s Y, and s Y= r It was originally thought that all plate glass was very flat and that its flatness must be controlled to very close tolerances during its manufacture. Conversations with production engineers at several glass factories revealed, however, that most emphasis is placed on the uniformity of thickness of the glass and not on its flatness. Essen tially, all standard plate glass is now produced by the flotation process. This gives it a very smooth surface and uniform thickness and it is simply assumed to be flat enough for most purposes. Numerous pieces of commercially available plate glass were optically tested for flatness, and it was found that some pieces are much flatter than others. The tests were run by placing the carefully cleaned surfaces of two glass plates together and illuminating them with monochromatic light from a sodium vapor lamp. The uniformity of the interference fringes produced gives a qualitative indication of the flatness of the two glass panes. Three pieces of special plate glass were obtained which were manufactured by the traditional process of mechanical grind ing and polishing. The manufacturer claimed that this glass was flat to within two interference fringes per inch, where one interference fringe represents a deviation from the mean surface of approximately 2.75 x 104 millimeters. When these pieces of plate glass were observed under the sodium vapor lamp, their superior flatness was evident. Selected pieces of the commercially available glass, however, were nearly as good. These were used in the construction of the HeleShaw model. In order to provide for the controlled infusion of the viscous liquids into the gap between the plates, it was necessary to glue a section of machined plexiglass to each end of both glass plates. In this way, an infusion reservoir was created on either end of the model so that the liquids could be injected uniformly into the test section. This is shown in Figure 4.2. The gap between the two glass plates was maintained at a uniform thickness by a series of brass spacers 1/16 of an inch (1.59 mm) thick placed at intervals around the edges of the model. The plates were held in position by Cclamps which pressed them tightly against the spacers. The model was sealed by a length of surgical rubber tubing which was placed around the edges, just inside of the spacers, and crushed flat by the action of the Cclamps. The wall thickness of this tubing was approximately 1/32 of an inch so that it could be crushed between the plates without exerting much force against the glass. The assembled model, mounted on its stand is shown in Figure 4.3. The length of the test section is 76 centimeters and its height is 38 centimeters. With the rubber seal in place, however, the actual height available for flow is 35.6 centimeters. The viscous liquids are supplied to the infusion reservoirs on either end of the model through plastic tubes. The plexiglass end pieces on the back plate are each supplied with four fittings for 1/2 inch (1.27 cm) plastic tubing. The heavier fluid is conducted to the model through this tubing from constant head tanks mounted on a ring stand. This is shown in Figure 4.4. On the front plate, the end pieces are provided with air vents on one end, and four nipples for the connection of 1/4 inch (0.64 cm) plastic tubing on the other. The lighter viscous fluid is injected through this tubing into one of the infusion reservoirs by a positive displacement syringe pump. This setup is illustrated in Figure 4.5. 40 V)u V) vJ uL j Lu V)A w  C C, Q ca co LU VA uLu CC w LLU crc '4 (A a:D .f (A (A F VD Cj V)u LL c LJ < ::D __j V) CL < I C)r c ULL)Iv ) 0 < > LLII I~ I. V 0 C I CC oc LL t) I II UJ C I fi LL 42 LU 1' 4 in ov U Lii V) Io mr LCC U j A w Properties of the Viscous Fluids Dow Corning Corporation Series 200 silicone oil was used as the heavier fluid in the model. This is a clear dimethyl siloxane which is characterized by oxidation resistance and low vapor pressure. Its specific gravity is 0.977 and its kinematic viscosity is nominally 500 centistokes at 250C. For the lighter fluid in the model, a mixture of 90 SAE and 30 SAE nondetergent motor oils was used. This mixture was adjusted until its kinematic viscosity was the same as that of the silicone oil at a temperature of 24.50C. The temperature versus viscosity relationships for the two fluids in the neighborhood of 250C is shown in Figure 4.6. This data was obtained with a CannonZietfuchs viscometer with calibration traceable to the National Bureau of Stan dards. The specific gravity of the light oil mixture was 0.890 at 250C. These two fluids are not readily miscible, yet the surface tension effects between the two are not so great that capillary action at the interface is noticeable. Operation of the Model The HeleShaw model is a versatile analytical tool and its mode of operation depends somewhat on the type of test being run. In the present case the object of the tests was to observe the shape and rate of rotation of the initially vertical interface between the two fluids of different density. The idea was that the behavior of such an interface would be similar to that of the interface formed by the injection of freshwater from a well into a confined saline aquifer, neglecting dispersion. Of course, the two processes are not the same because the flow away from a well is radial while the HeleShaw model MOTOR OIL MIXTURE DOW 200 FLUID 2 23 24 25 2 TEMPERATURE (DEGREES C.) FIGURE 4.6. TEMPERATURE VS. VISCOSITY FOR THE TWO VISCOUS FLUIDS 540 530 520 510 500 is limited to twodimensional flow in cartesian coordinates. Direct application of the HeleShaw results to a prototype well situation is, therefore, not possible. With increasing distance from the well, however, the two flows become more and more similar so that the model tests are useful for the qualitative description of interface movement. The results of the model tests are also useful for comparison with numerical models which can be made to work in either radial or cartesian coordinate systems. In setting up the model for a test, it is necessary first to be sure that all parts of it are clean and dry. This can only be done by completely disassembling it and washing each part. A commercial pre paration of trichloroethane cleaning solution was found useful in washing the oils from the glass and plastic parts. After the model had been cleaned and reassembled, with the plates rotated to the horizontal position, the heavier oil was introduced at the downstream end; the upstream end being the one into which the lighter oil would be injected. During this process the upstream end was raised a few centimeters so that all of the air in the test section would be expelled. Air entrapped in the downstream infusion reservoir was released by loosen ing the brass plugs provided for this purpose. The introduction of silicone oil was continued until the oil reached the end of the glass and was just ready to run into the upstream infusion reservoir. At this point, the flow was stopped and the oil front was held at the edge of the glass by proper adjustment of the constant head supply tank. The 1/2 inch (1.27 cm) plastic tubes on the bottom side of the upstream end were then clamped off tightly with hose clamps. Three of the 1/4 inch (0.635 cm) tubes from the syringe pump were then hooked up and the upstream infusion reservoir was filled with motor oil. The fourth nipple was left open for the escape of air. After all of the air had escaped and the model was entirely filled with the two viscous fluids, the fourth tube was attached and the interface between the two oils was driven into the model to the starting point for the test. It was found helpful to let the interface sit in this position for several hours so that any irregularities in it would be straightened out by gravity. The elevated end of the model was then lowered to a level position and the test begun. The test was started by simultaneously rotating the model to the vertical position, turning on the syringe pump and starting the timer. The movement of the interface could then be photographed until the downstream end of the model was reached. Several variations of this procedure were used. The first test was run with the model left in its horizontal position so that the straightness of the advancing interface could be checked. It was found that, with higher rates of injection of the lighter fluid, the interface moved faster in the center of the model. This indicated that the higher pressure required to drive the viscous fluids across the model could cause the glass plates to bend, changing the thickness of the gap. Figures 4.7 and 4.8 show the shape of the interface when advancing at rates of 3.65 cm/min and 7.28 cm/min, respectively. Since the straightness of the interface at the lower speed was considered acceptable most of the remaining tests were run at this rate of injection. Figures 4.9 and 4.10 show a test which was run at a higher injec tion rate with the model in the vertical position. The bowing out of the glass plates caused the center section of the interface to advance 48 It hi. 0 i .i (J N i t L.. o. 49 ~~CD uj ~' tLw CD I C) F LLJ iy ::Do LL S C 0 f ~C'  I('4 Co II H C, LU r) C, L~,LU tc LA rr~j !gta ~ 1 "~1~ 1""a:: lrLii H r, IILL F a C J ' ~~cc aI I: U) S, ' I i i ~i~b LsParC ~ I ~r~02~ ' k' ; LL 51 II ' I co CD A 17) 7 mmm.L' ~i 7p,) L) L4 Lit!i C Li I t r . 'wE f more rapidly than the ends. Rotation of the upper part of the interface past the vertical results in instability, which is clearly evident here. The sequence of Figures 4.11 through 4.16 shows a test that was run with an injection rate of 20.44 cubic centimeters per minute. This is the same injection rate that resulted in the relatively straight interface shown in Figure 4.7 when the model was horizontal. The sequence of Figures 4.17 through 4.23 shows a test in which only the lower half of the model was used. This was accomplished by inserting a barrier composed of rubber splicing compound and silicon rubber sealant between the glass plates to reduce the width of the test section by half. The injection rate used here was 8.16 cubic centi meters per minute. Analysis of Test Results The rotation of an initially vertical interface between two fluids of different density has been studied in several previous investigations. Gardner, Downie and Kendall (1962) proposed an equation to predict the rate of interface rotation which, for an isotropic aquifer and fluids of equal viscosity, has the form ()2 16 (t/t')2 (4.31) 3 1 + t/t' where X = projection of the interface on the horizontal b = thickness of the aquifer t = time t' = 4/3 nbp KAp p = average density of the two fluids. *1~t !  I ~ ~ 4  I ' 1 ~~1 jI, ji 4. 7 .i13 r4 II C 0 I' LL_ I: I U I I LU NE (I i ~' i '." I~ sl * 1 2 It * 1 i . hi 1 1 it I, I * 1~ * 4 , i I' " *.~ 7. ' 1'L~lr I 4 4 z C) d cx 0 F LU < a: I 0 z: I cC LL Ii ' IItr jr t 4 tz St, .. I'' , t , C k ...: . "t~fi e I 4 j / ;l', ir ! i [ i ' I ~ ;~ rP4 t . t .. ' .LU r 4 j ,r , F 1 1 V , 4. I I .C 4..., II  *1 *. J I ; l I :, I ' iiin i V  ...... L   4 If.. i ' ii . : .: 1. 1; I 4 1 ". 1 .. ,. I .. 1 . ,,,"4 m* I I a:J '~ LI Cr I I j74.i ~~dls 1ik 4iL a: IYC'LFlr LU H r t Iv  I~~* >40c )fi ~t< C' 1.2 3 1* * IL J . I I '1, t t 1.1. "I a'.i :.1 I #sJS i b0y jL~ i i,, I' ^ , S Cj *1 C)' ItI 1~ Ci  I t I :Fi 44U41 t rS , CP LL =D . . ,t  p 58 4' c1 L U1 iN Vc S I 41  Cf) LUJ V) LUJ jo LUJ IL C0 F;~j i iLL. 59 00 C, C) I (L, < m 14 IV z LU 7o LL !C) CL LLUj i' I 4. LUJ D il i Cl, UL .L LJ 60 co ,T A rII O~LL ~ Cf. P IcI F < LLu Li u Lu J CK U ui u L 2 E i~ Lu Lu Li'.. 61 00 1 II ILU 2t7T9.IILL. V; u Lii ` `i.  ~lLLJ QI (!J Li) i' r Fi' cc I I' I c W i I !I i ~ Uo 1: i~i L 62 II W. co II rc C)2 LU  '4 C) LU Z LU L: IL rll cOtf F 4r Ar,, 63 C6) C) LU co Li.. fi *D LUJ rl LL LUI Li.. LUJ LUI LL 64 z 06 *LUJ CA) LU liJJ 'ti~l~~1'4 t LLJ '4 oJ C'3 LU t =D L J a: LLJ L J .j r L Li o LL r, s i r'll~ = 1.4 65 i;L.) c tI '" "! liLl L~LU tii 4,i F t J LL uj CD g ii~: z > r .;7:~ cLl Z;S 4' ~ r~tY~dig1138 ~ s i LL tl~~ If~g~~ i I vCY ZD 't) ICD LL. Equation (4.31) was constructed intuitively by joining the solutions to the extreme cases of nearly vertical and nearly horizontal interfaces. The solutions to these extreme cases were based on the assumptions of no vertical flow and no curvature of the interface. Another equation for the movement of an initially vertical inter face was given by Harleman and Rumer (1962). Once again, assuming a straight interface and no vertical flow, they obtained X 2t S = 8/3 (4.32) This equation is of the same form as the one for nearly horizontal interfaces from which equation (4.31) was derived. Figure (4.24) shows the comparison of equations (4.31) and (4.32) with data taken from three runs of the vertical HeleShaw model. Runs number two and three are shown in Figures 4.11 through 4.16 and 4.17 through 4.23, respectively. Run number one was performed in the same manner as run number two except that the injection rate was twice as high. The third run gave the best fit to equation (4.31). Runs number one and two, while fairly consistent between themselves, do not fit either equation very well. Since the main difference between these two runs and the third run is in the depth of the flow region, it is most logical to suspect that it is caused by the existence of vertical velocity components which were ignored in the derivation of the equations. The greater curvature of the interface which is apparent in the first two runs tends to support this view. 67 3 I I 2 S.8 I "Y  T 1  i t . .... .. / " ' .. .. . HARLEMAN & RUMER .4 GARDNER, DOWNIE & KENDALL .2 SRUN 1 o RUN 2 .I RUN 3 .04 .06 .08 .1 .2 .4 .6 .8 1 t/t' FIGURE 4.24. COMPARISON OF EXPERIMENTAL RESULTS WITH EQUATIONS FOR INTERFACE ROTATION CHAPTER 5 SIMPLE CASES INVOLVING STEADY FLOW General It was shown in Chapter 3 that the general problem of locating a moving interface between two fluids of different density is very complex even when the standard simplifying assumptions involving homogeneity, isotropy and immiscibility are used. It requires the solution of two coupled flow problems, one in the freshwater and one in the saltwater, joined by a nonlinear boundary condition on the interface, of which the position is unknown. The only methods presently available for getting a solution to such a problem are by physical and numerical modeling. However, if the problem is limited to the special case of a stationary freshwater lens in a static saline aquifer, then it becomes much more tractable. In this special case, it is only necessary to solve the continuity and motion equations in the freshwater region and the location of the interface, while still not initially known, can be treated by the GhybenHerzberg principle. This makes it possible to formulate some problems which can be solved analytically. Steady flow in a freshwater lens is not possible unless the aquifer possesses a constant head boundary which permits freshwater to flow out of the lens at the same rate that it is being injected. One such constant head boundary is found at the sea coast where injection wells may be used to create freshwater lenses in aquifers affected by saltwater intrusion. Another case is that of injection into a leaky saline aquifer where the leaky layer serves as a boundary which permits efflux of freshwater from the lens. Even in these cases, steady flow may only occur during a small part of the history of the lens. It represents a condition of equilibrium in which buoyant forces are counteracted by the loss of energy due to the flow resistance of the porous medium. It is useful to study this condition because it indicates the limiting size of lens that can be achieved under the given flow conditions. It should also be noted that, in the operation of a storage lens, certain advantages can be gained by maintaining the lens in a steady state through continuous freshwater injection. The first is that continuous injection is necessary to preserve the potential energy of the lens so that the maximum amount of freshwater can be withdrawn when it is needed. Note that as soon as injection is stopped, the lens begins to decay into a thin layer of freshwater at the top of the aqui fer which could not be selectively retrieved even if fresh and saltwater did not mix. The second advantage is that the continuous injection of freshwater to the lens tends to wash out the saltwater which had re mained in the pores and adhered to the particles of the porous medium after the passage of the mixing zone as the lens was being formed. In this way, the salinity of the freshwater and the thickness of the mixing zone are reduced. Freshwater Injection in a Coastal Aquifer Saltwater Intrusion It is commonly observed in coastal aquifers which discharge to the sea that a wedge of nearly static saltwater extends inland from the coast underneath the flowing freshwater. Actually, it was the study of this phenomenon which led W. Badon Ghyben and A. Herzberg to develop their hydrostatic theory for interfaces. Numerous subsequent investi gators have brought forth more refined methods for analyzing saltwater wedges (Glover, 1959; Henry, 1959; Bear and Dagan, 1964; Van Der Veer, 1977) and have found the GhybenHerzberg results to be quite accurate except for that part of the interface closest to the sea. The application of this principle, combined with the Dupuit assump tion, to a confined aquifer of constant thickness gives the following relationship between h and x. x = Kh2/(26Qx) (5.1) in which the symbols are as shown in Figure 5.1. The position of the toe of the interface is given by Kb2 x = 2 (5.2) t 26Q Injection Into a Coastal Aquifer If an injection well pumps freshwater at a constant rate into the coastal aquifer pictured in Figure 5.2 then the equilibrium position of the saline wedge will be modified and a steady freshwater lens will be formed in what was formerly a saline part of the aquifer. This is a case of injection into an aquifer in which only the freshwater is flow ing, so that the GhybenHerzberg principles can still be used. The flow pattern around a well in a uniformly flowing aquifer is customarily studied by the methods of potential theory in which the velocity poten tials for source flow and uniform flow are superimposed. However, in the presence of the saltwater wedge, even though the vertical velocity components are neglected, flow is not two dimensional because of the 71 o N I  L. re" I LA ' \LUl i tn / u I \ I Cz I u , O I I L I < Lu Lu I SICD . i (1 I y I X  I I C I II  N I ~U ; I L . u I e8 72 Sll* LL  LX / I LiiL N 0j LL 0 i ^ /, ^  '" ^ L .I.I I L.C) Lu "J I \ L I \  U u L\ / I   O ^ ' changing thickness of the flow region. Strack (1976) has found a modified form of the vertically integrated potentials originally pro posed by Girinskii (see Bear, 1972; p. 157) which can account for the presence of the interface. If the thickness of the flow region is a linear function of the piezometric head, 4, i.e., h = #i + 8 (5.3) then Strack gives the potential for the interface zone as Di = 1/2aK(4 + B/a)2 (5.4) For the inland part of the aquifer where there is no interface, the potential is defined by P = Kb) + c (5.5) where the constant, c, is used to make the potential function continuous at the toe of the interface. It is easily verified from equations (5.4) and (5.5) that the discharge components in the x and y directions are obtained from 0 in the usual way. Qx = a3/ax and Qy = D/ay (5.6) in which Q and Q represent the discharge through a unit width and the x y total thickness of the aquifer in the x and y directions, respectively. From equations (5.6) it follows that the continuity condition for the confined coastal aquifer is given by the Laplace equation a20/ax2 + a2o/ay2 = 0 (5.7) and 0 is a harmonic function of x and y. The effect of this is that through use of Girinskii's potential, variations in the vertical direc tion have been integrated out of the problem. The resulting potential function can be manipulated, using the theory of complex variables to describe a twodimensional flow field with the desired boundary condi tions. At the same time, the numerical value of the potential function at any point gives the depth to the interface at that point through the GhybenHerzberg relation and thereby supplies threedimensional information about the freshwater lens. Before this can be done, however, it is necessary to determine the constants a and B in equation (5.4) and c in equation (5.5). Confined coastal aquifer. From the geometry of the saltwater wedge, as shown in Figure 5.2, it is seen that s = b a (5.8) The thickness of the freshwater region, h, is related to the "push up", s, through h = 6s a (5.9) which is a consequence of the GhybenHerzberg principle. The combination of equations (5.8) and (5.9) results in h = 64 6b (1 + 6)a (5.10) But it has already been assumed that the relationship between h and < in the interface zone is given by equation (5.3). Consequently the values of a and B must be: a = 6 and 0 = 6b (1 + 6)a (5.11) Substituting these values into equation (5.4) gives the definition of the potential function in the interface zone for a confined coastal aquifer. i 6K[p b (1 + 1/6)a]2 (5.12) Now, in order that the potential function may be continuous throughout the aquifer, it is necessary to determine c in equation (5.5). Note that at the toe of the saltwater edge s = (a + b)/6 (5.13) which along with equation (5.8) gives : = (a + b) (1 + 1/6) (5.14) Using this value of p in both equations (5.5) and (5.12) and requiring that they be equal at the toe of the wedge results in c = 2 Kb(a + b) (1 + 1/6) (5.15) The definition of the potential function in the inland zone of a confined coastal aquifer then becomes P = Kb[b + 1/2 b/6 (a + b) (1 + 1/6)] (5.16) At the toe of the saltwater wedge this reduces to (t = 1/2 Kb2/6 (5.17) In order to construct the complex potential function for flow in the aquifer the boundary condition at the sea coast must be known. From Figure 5.2 it is seen that s = a/6 at the coast. Substituting this into equation (5.12) gives the boundary condition ( = 0 at x = 0 (5.18) This condition can be satisfied by locating an image well, with the same strength as the injection well but opposite sign, on the other side of the boundary at the same distance from it, x,, as the injection well. The superposition of the velocity potential for the two wells and the uniform seaward flow results in (x Xw) y22 = xQ + In ( )2 y2 (5.19) x 47r (x + x) 2+ y2 When particular values of 0 are computed from equation (5.12) and subs tituted into equation (5.19) which is then solved for the particular equipotential line, the shape of the freshwater lens produced by the injection well is outlined. Such a case is shown in Figure 5.3. From a practical point of view it is most important to locate the toe of the saltwater wedge around the periphery of the lens because this is the outer boundary of the usable portion of the freshwater. The equipotential line which delineates this outer boundary has the value given by equation (5.17). Solutions for this line can be made dimen sionless by defining the following dimensionless parameters: x' x (5.20a) t y, = Qx y (5.20b) It x Q x (5.20c) x w t w Q. = Q (5.20d) (t A series of such dimensionless solutions for various conditions of injection is shown in Figure 5.4. and Figure 5.5. .) z cr S V SI V 'I 0 C\J r. II Lfl SII i" C) LU (5 I LLI Z Im LL Lu C) Lu j cr: VL I4 2: LJ LII LO C, L u o Cn x($/Ix)= ,x Lfl LO II U) X(,/Xb)= ,x Phreatic coastal aquifer. The Girinskii potentials defined by Strack for interface flow can also be used in aquifers which have no upper confining layer because the thickness of the freshwater region is still a linear function of the piezometric head. The linear function is not the same as for a confined aquifer, however, so the potentials must be redefined to conform to the aquifer geometry shown in Figure 5.6. In the inland zone the specific discharge, Qx is given by K2 Qx K K2 (5.21) ax The potential for this region, then, should be defined as D = 1/2 KD2 + c (5.22) In the interface zone, note that s = ( b (5.23) and from the GhybenHerzberg principle the thickness of the flow region is h = s(l + 6) (5.24) which leads to h = (1 + 6)( (1 + 6)b (5.25) If this is to be in the form of equation (5.3) then the constants a and 6 must be a = (1 + 6) and B = (1 + 6)b (5.26) Substituting these values into equation (5.4) results in a definition 81 UJ rw LU C, 0  LLJ < /c SI LJ l LU I_ LU < r _j I u. SLLU I I I ULU LU 0 II F 4  r \ 0 D< LU I~ 4 0Z for the velocity potential in the interface zone Di = 1/2 (1 + 6) K(q b)2 (5.27) At the toe of the saltwater wedge the piezometric head is t = b(l + 1/6) (5.28) Substituting this into equations (5.22) and (5.27) and requiring that be single valued at this point leads to the value of the constant c in equation (5.22). c = 1/2 Kb2 (1 + 1/6) (5.29) The definition of the velocity potential in the inland zone of the phreatic aquifer .then becomes S= 1/2 K[P2 (1 + 1/6)b2] (5.30) At the toe of the saltwater wedge this reduces to Dt = 1/2 Kb2(l + 6)/62 (5.31) which is greater than the corresponding value for a confined aquifer by a factor of (1 + 1/6). Axisymmetric Flow in a Leaky Aquifer Consider a homogeneous saline aquifer with hydraulic conductivity, K, of uniform thickness, b, and unlimited area extent. It is confined at the bottom by a horizontal impervious layer and at the top by a horizontal aquitard of thickness, b', and hydraulic conductivity, K'. A fully penetrating well injects freshwater at a constant rate, Q, forming a freshwater lens, which remains centered on the well since there is no preexisting flow in the surrounding salt water. Injection continues until a steady state is eventually reached in which the injec tion rate is matched by the rate of leakage through the aquitard. The mixing zone at the boundary of the lens then becomes relatively thin and can be approximated by a sharp interface. This situation is depicted in Figure 5.7. The flow region is divided into two zones, each of which has its own governing equation. In zone 1, which is occupied entirely by freshwater, the customary analysis for well flow in leaky aquifers can be used. This includes the assumption of 90 degree streamline refrac tion at the aquitard so that flow in the aquifer is horizontal while in the aquitard it is vertical. The governing equation is (see De Weist, 1965, p. 264) d2s 1 ds 1 d2 + d )2 s = 0 (5.32) dr2 r dr where B2 = Kbb'/K' This is a modified Bessel equation of order zero, for which the general solution is given by s = C1I0(r/B) + C2KO(r/B) (5.33) where 10 and K0 are the modified zeroth order Bessel functions of the first and second kinds, respectively. The boundary condition at the 