Citation |

- Permanent Link:
- http://ufdc.ufl.edu/UF00084179/00001
## Material Information- Title:
- A study and evaluation of saltwater intrusion in the Floridan aquifer by means of a Hele-Shaw model
- Creator:
- Evans, Andrew Joseph, 1940-
- Publication Date:
- 1975
- Language:
- English
- Physical Description:
- xiv, 213 leaves : ill., diagrs., maps ; 28cm.
## Subjects- Subjects / Keywords:
- Aquifers ( jstor )
Cradles ( jstor ) Fresh water ( jstor ) Groundwater ( jstor ) Limestones ( jstor ) Modeling ( jstor ) Parametric models ( jstor ) Pumps ( jstor ) Saltwater ( jstor ) Velocity ( jstor ) Civil Engineering thesis Ph. D Dissertations, Academic -- Civil Engineering -- UF Saltwater encroachment -- Florida ( lcsh ) - Genre:
- bibliography ( marcgt )
non-fiction ( marcgt )
## Notes- Thesis:
- Thesis--University of Florida.
- Bibliography:
- Bibliography: leaves 210-212.
- General Note:
- Typescript.
- General Note:
- Vita.
- Statement of Responsibility:
- by Andrew Joseph Evans, Jr.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Resource Identifier:
- 000158804 ( ALEPH )
02609315 ( OCLC ) AAS5126 ( NOTIS )
## UFDC Membership |

Full Text |

A STUDY AND EVALUATION OF SALTWATER INTRUSION IN THE FLORIDAN AQUIFER BY MEANS OF A HELE-SHAW MODEL By ANDREW JOSEPH EVANS, JR. 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 1975 ACKNOWLEDGEMENTS To my supervisory committee, I wish to extend my sincere appreciation for their efforts on my behalf. A special thanks goes to Dr. B. A. Christensen and Dr. J. H. Schaub for their continuing support and advice. It's not an easy task to pull an oar in Hagar's longboat. The efforts of Mr. C. L. White, Mr. William Whitehead, and the staff of the Mechanical Engineering machine shop in the construction of the model are grate- fully acknowledged, as well as the assistance of Mr. Richard Sweet and Mr. Tom Costello'in the operation of the model, and the mathematical advice of Dr. Jonathan Lee. The participa- tion of Mr. Floyd L. Combs in all aspects of the project was especially valuable. Without the financial support of the Office of Water Resources Research, United States Department of the Interior, and the administrative assistance of Dr. W. H. Morgan and the staff of the Florida Water Resources Research Center, this project would not have been launched. Finally, to my mother and my many dear friends who have watched the "ins and outs" of my academic and profes- sional career, bless you for your interest and encouragement. I want to leave the reader with this closing thought, which I believe is paraphrased from Benjamin Franklin: "Waste not want not, you never miss the water 'till the well runs dry." iii TABLE OF CONTENTS Page ACKNOWLEDGEMENTS ...................................... ii LIST OF TABLES ........................................ vi LIST OF FIGURES .... ........................... ........ vii KEY TO SYMBOLS OR ABBREVIATIONS ........................ x ABSTRACT .............................................. xiii Chapter I. INTRODUCTION AND STATEMENT OF PROBLEM ......... 1 Topography ................................. 2 Western Highlands ....................... 3 Marianna Lowlands ........................ 6 Tallahassee Hills ....................... 7 Central Highlands ....................... 7 Coastal Lowlands ........................ 7 Climate .................................... 8 Geology .................................... 12 II. MODELS, NUMERICAL AND PHYSICAL ................ 15 A: NUMERICAL METHODS ............................. 15 Method of Finite Differences ................ 15 Method of Finite Elements ................... 17 Relaxation Methods ......................... 18 B: PHYSICAL MODELS ............................... 19 Sandbox Model .............................. 21 Hele-Shaw Analog ............... ............. 22 Electric Analog .............................. 23 Continuous Electric Analog .............. 23 Discrete Electric Analog ................ 24 Ion Motion Analog ....................... 24 Membrane Analog .............................. 25 Summary .................................... 26 TABLE OF CONTENTS-Continued Chapter Page III. THE HELE-SHAW MODEL ........................... 28 Viscous Flow Analog ........................ 29 Scaling .................................... 37 Time .................................... 39 Anisotropy .............................. 40 Leakage ................................. 50 Storativity ............................. 54 Discharge ............................... 55 Accretion ............................... 59 Volume .................................. 59 IV. SITE SELECTION AND PROTOTYPE GEOLOGY AND HYDROLOGY .............................. 61 Site Selection ............................. 61 Prototype Geology .......................... 67 Prototype Hydrology ......................... 69 V. DESIGN, CONSTRUCTION, AND OPERATION OF MODEL ................................... 71 Design ..................................... 71 Prototype ............................... 71 Model ................................... 71 Construction ............................... 86 Frame ................................... 86 Plexiglas Plates and Manifolds .......... 105 Saltwater System ........................ 125 Freshwater System, General .............. 125 Freshwater System, Accretion ............ 132 Freshwater System, Wells ................ 132 Freshwater System, Flow Meters .......... 132 Operation ........ .......................... 156 VI. RESULTS, CONCLUSIONS, AND RECOMMENDATIONS ..... 158 Results .............. ....................... 158 Conclusions and Recommendations ............ 191 Appendix A. UNIFORM FLOW THROUGH A CONDUIT OF RECTANGULAR CROSS SECTION .................. 193 B. FLOW METER CALIBRATION ........................ 195 BIBLIOGRAPHY .......................................... 210 BIOGRAPHICAL SKETCH ................................... 213 LIST OF TABLES Table Page 2.1 APPLICABILITY OF MODELS AND ANALOGS .......... 27 5.1 PROTOTYPE PARAMETERS ......................... 72 5.2 MODEL PARAMETERS ............................. 81 5.3 SIMILARITY RATIOS ............................ 82 6.1 DEPTH TO SALTWATER ........................... 186 LIST OF FIGURES Figure Page 1.1 Topographic divisions of Florida ................ 5 1.2 Mean annual precipitation ...................... 11 3.1 Free body flow diagram for Hele-Shaw model ..... 32 3.2 Section of anisotropic grooved zone in Hele-Shaw model ............................. 44 3.3 Head loss in grooved anisotropic zone ........... 46 3.4 Section of leaky zone in Hele-Shaw model ....... 53 3.5 Storage manifold for Hele-Shaw model ........... 57 4.1 Regional area of prototype ..................... 64 5.1 Ghyben-Herzberg interface model ................ 75 5.2 Cradle and cradle dolly ........................ 88 5.3 Cradle rotation ................................ 91 5.4 Frame and model setup .......................... 93 5.5 Air hose and valve arrangement ................. 96 5.6 Stub shaft and pillow block arrangement ........ 98 5.7 Back up air supply .............................. 100 5.8 Internal support and sealing system ............ 102 5.9 Model mounting system .......................... 104 5.10 Model back and front plates .................... 107 5.11 Detail of the model front and back plate ....... 109 5.12 Detail of the model back plate ................. 111 vii LIST OF FIGURES-Continued Figure Page 5.13 Front plate with accretion manifolds ........... 114 5.14 Accretion manifolds ............................. 116 5.15 Detail of accretion manifolds .................. 118 5.16 Connections between model and fluid supply system ............................... 120 5.17 Saltwater constant head tank ................... 122 5.18 Back and front plate clamp up .................. 124 5.19 Fluid supply network schematic ................. 127 5.20 Saltwater reservoir and pump ................... 129 5.21 Freshwater supply system ....................... 131 5.22 Freshwater reservoir and accretion pump ........ 134 5.23 Well supply manifold and pump .................. 136 5.24 Well supply manifold and accretion supply manifold ............................. 138 5.25 Opposite view of Figure 5.24 ................... 141 5.26 Flow meter bank ................................ 143 5.27 Flow meter detail .............................. 145 5.28 Flow meter to model connections ................ 148 5.29 Flow meter pressure sensing lines .............. 150, 5.30 Flow meter switching device and pressure transducers ........................ 152 5.31 Carrier demodulator and strip chart recorder .................................... 155 6.1 Interface location, tm = 0 min. ................. 160 6.2 Interface location, tm = 32 min. ............... 162 6.3 Interface location, tm = 48 min. ................ 164 viii LIST OF FIGURES-Continued Figure 6.4 6.5 6.6 6.7 6.8 6.9 6.10 6.11 6.12 6.13 location, location, location, location, location, location, location, location, location, location, tm = 62 min. t = 79 min. t = 92 min. m tm = 102 min tm = 107 min t = 117 min t = 132 min m tm = 147 mir tm = 162 mir tm = 172 mir L. L. t. L. 6.14 Wedge location at time zero ..... Page ............... 166 ............... 168 ............... 170 .............. 172 ........ ..... 174 ..... ......... 176 .............. 178 .............. 180 .............. 182 .............. 184 ............... 188 6.15 Time overlay of wedge location .. Interface Interface Interface Interface Interface Interface Interface Interface Interface Interface 190 KEY TO SYMBOLS OR ABBREVIATIONS Symbols A = area, Fourier coefficient b = width of model and prototype B = body force, Fourier coefficient f = subscript denoting freshwater, Floridan aquifer g = acceleration due to gravity h = height of water j = summation limit k = intrinsic permeability .K = hydraulic conductivity 1 = distance between storativity tubes, subscript denoting leaky layer L = distance center to center of groove in anisotropic zone, flow meter tube length m = subscript denoting model n = effective porosity p = pressure, subscript denoting prototype q = specific discharge Q = total flow r = subscript denoting ratio R = accretion R = Reynolds number KEY TO SYMBOLS OR ABBREVIATIONS-Continued Symbols S = specific storage t = time T = transmissivity U = volume V = velocity x = horizontal direction parallel to test section y = horizontal direction perpendicular to test section z = vertical direction 1,2 = subscripts denoting zone 1 and zone 2 a = width adjustment factor y = unit weight A = length adjustment factor p = absolute viscosity E = geometric parameter in anisotropic zone p = mass density v = kinematic viscosity = potentiometric head = velocity potential, potential Abbreviations D. substantial derivative Dt V2 = Laplace operator A = difference operator C = Centigrade KEY TO SYMBOLS OR ABBREVIATIONS-Continued Abbreviations cfs = cubic feet per second *F = Fahrenheit fps = feet per second g/cm3 = grams per cubic centimeter gpd = gallons per day Hg = Mercury I.D. = inner diameter MGD = million gallons per day msl = mean sea level O.D. = outer diameter pcf = pounds per cubic feet psid = pounds per square inch differential psig = pounds per square inch gage RPM = revolutions per minute xii 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 A STUDY AND EVALUATION OF SALTWATER INTRUSION IN THE FLORIDAN AQUIFER BY MEANS OF A HELE-SHAW MODEL By Andrew Joseph Evans, Jr. August, 1975 Chairman: B. A. Christensen Major Department: Civil Engineering Continuing development of the coastline zone in the middle Gulf area of Florida is increasing the demand for groundwater supplies and, in turn, increasing the probability of saltwater intrusion. Methods must be developed to make long-range predictions on the effects of increased demands on the Floridan aquifer. A Hele-Shaw model is a physical model which fits the requirements for long-range planning. It is well-suited for handling anisotropic aquifers, difficult boundary conditions and can simulate years of field conditions in minutes of model time. xiii The section selected for study lies in a line from the Gulf coast near Tarpon Springs to a point near Dade City and passes through the Eldridge-Wilde well field. The Eldridge-Wilde well field is the major water producer for Pinellas County. This region has experienced several years of dry weather, and pumping has lowered the water levels in the aquifer by a significant amount. This loss of fresh- water head is certain to induce saltwater intrusion. A Hele-Shaw model has been built for this area, and all pertinent geological and hydrological features of the area are included. Steady state characteristics of the aquifer system have been considered. In particular, the long-term effects due to pumping and artificial recharge were examined. xiv CHAPTER I INTRODUCTION AND STATEMENT OF PROBLEM Florida, with the possible exception of California, is the fastest growing state of the United States. The rapid influx of people since World War II has greatly increased the demands for land and water. In the past, there has been an almost total lack of wide range planning for the uses of these resources. Even fewer investigations have been made into the consequences of their rapid and unordered development. Recently, water supplies have had to be rationed in South Florida. Overall, land and wet- lands required for fish and wildlife have so diminished that, in some instances, there has been a marked decrease in their numbers. It seems reasonable to conclude that in some areas of the state, land and water resources cannot support much larger populations with current locally avail- able supplies without almost irrecoverable damage to the groundwater system in the form of saltwater encroachment. With the growing affluence of the American people, and the availability of economically priced air conditioning units, it can be expected that even more people will leave the colder northern climates for the southern and western states. Florida can expect to receive more than its share of the migration. Frequently now, environmental protection groups are making forecasts of impending doom. At worst, their predictions may come true and people are beginning to look at all growth with a jaundiced eye. It is doubtful, however, that growth and development can be stopped. The history of man indicates a continual effort to better his life-style, his private environment. There is little doubt that this has sometimes caused a degradation of other portions of his world. Unless the cessation of all growth and development is acceptable, new ways must be found of forecasting, or predicting, the results of all growth so as to combat possible undesirable results. Consequences of all growth must be known, even of those resulting when pragmatic short-term solutions are used. Hopefully, the remainder of this study will present a modeling method which will be useful in forecasting the results of pumpage and use of groundwater in our coastal zones so that we may better plan their usage. But first, a little background on Florida. Topography Florida lies between latitudes 24-40' and 31-00' north, and longitudes 80-00' and 87*-40' west, and is the most southerly unit of the continental United States. In its southernmost extension it is less than 10 of latitude north of the Tropic of Cancer. Florida is bounded on the east by the Atlantic Ocean; on the south by the Straits of Florida and the Gulf of Mexico; on the west by the Gulf of Mexico and the state of Alabama; and on the north by Alabama and Georgia. The shape of the state in relation to the remainder of the United States suggests two distinctive parts: the Floridan panhandle and the peninsula of Florida. The panhandle is a strip roughly 225 miles long that stretches in an east-west direction. The peninsula is a south-southeast extension at approximate bearing S 17 east. From the northern boundary of the state to the tip, not including the chain of keys, the peninsula is approximately 415 miles long and includes two-thirds of the land mass of the entire state. Its coastline, some 1350+ miles long, is the longest with the exception of Alaska. No place in the interior of Florida is more than 60 miles from either the Gulf or the Atlantic coast. Cooke (1939) divided the terrain of Florida into five sections, Figure 1; the Western Highlands, the Marianna Lowlands, the Tallahassee Hills, and a narrow band of Coastal Lowlands, which comprise the panhandle, and the Central High- lands and Coastal Lowlands, which comprise the peninsula. The topography of each is described briefly. Western Highlands Extending eastward from the Perdido River (the western boundary of Florida) to the Apalachicola River, the northern part of this region near the Alabama state line FIGURE 1.1 Topographic divisions of Florida. Western Highlands Marianna Lowlands Tallahassee Hills Central Highlands Coastal Lowlands 86 0 50 100 150 Scale of Miles After: Cooke (1939) is not much higher than 300 feet. It is considered to be hilly when compared to the broad gently rolling southern parts of this region which drop to 100 feet elevation as one approaches the coastal lowlands. The highest elevation in the state, 345 feet, is found in this region in the north- west corner of Walton County. The western highlands are underlain with the sand of the Pliocene Citronelle Formation. The steepness of the bankslopes at the headwaters of the many streams is the most unique physiographic characteristic of this section. Marianna Lowlands This roughly triangular-shaped region of Holmes, Jackson, and Washington counties, with somewhat smaller contributions from Bay and Calhoun counties, lies between the Tallahassee Hills and the Western Highlands. It is difficult to recognize this area of gently rolling hills as lowlands. Cooke (1945) attributes the lower elevations to the solubility and consequent degradation of the under- lying limestone. This area is one of the two in the state where the Ocala Formation is exposed to the surface and the only area of the state where the Marianna limestone, the soft white limestone of the Oligocene Age, is found exposed. The region is dotted with sinks, sinkhole lakes, and springs. Tallahassee Hills From the Apalachicola River east to the Withlacoochee River, the Tallahassee Hills extends along the Georgia- Florida border and is only 25 miles in width. The western section is a nearly level plateau some 300 feet above mean sea level. The remainder consists of rolling hills carved out of the Citronelle Formation. In addition to this, a red clayey sand and Fuller's earth of the Hawthorn Formation are found in this area. This is a fertile farming region. Central Highlands The Central Highlands forms the backbone of the Floridan peninsula and extends from the Georgia line between the Withlacoochee and St. Mary's rivers south-southeastward some 250 miles into Glades County west of Lake Okeechobee. This region is highly diversified. It includes high swampy plains, hills, and innumerous lakes. Soils are sandy. Many of them were derived from Pleistocene (Ice Age) marine ter- races. However, a distinguishable amount comes from the Miocene Hawthorn and Pliocene Citronelle Formations. The lakes and sinks which dot the entire area indicate the presence of limestone below the surface. Elevations of this region average just slightly more than 150 feet; however, they vary from less than 100 feet to approximately 300 feet. Coastal Lowlands The Coastal Lowlands, or Coastal Plains as it is sometimes called, borders the entire 1350 mile Florida coastline. Flanking on both sides of the Central Highlands, the Coastal Lowlands is widest just south of Lake Placid and narrowest between the western border and the Choctawhatchee Bay just south of the Western Highlands. The elevations everywhere within this region are less than 100 feet. The soil for the most part is sandy except in the Everglades and Big Cypress Swamp locales where Pliocene limestone, muck, and peat prevail near the surface. The keys, which extend some 100 miles into the Straits of Florida, are mostly sandy oolitic limestone like that of the mainland; however, some limestone with coral heads is found. The islands seldom reach 15 feet elevation. The entire region is generally flat, typical of recently deposited material with little or no erosion. Climate The sea surface temperatures east and west of Florida average, respectively, 780 and 77 Fahrenheit. Water tem- peratures range from 740 to 83 Fahrenheit in the east and 700 to 84* Fahrenheit in the west. The coldest month in both cases is February; the warmest month is likewise August. The relatively homogeneous distribution of sea temperature, the lack of high relief and the peninsula shape of Florida contribute to its climate. The temperature is everywhere subtropical. Mean annual average temperature in the north is 68 Fahrenheit, and in the southern tip 750 Fahrenheit. The tradewinds, which shift from northern Florida to southern Florida and back semiannually, bring a mildly mon- soon effect to Florida. In November, the tradewinds are at their southernmost extension and Florida's climate is con- trolled by frontal, or cyclonic, activity moving in from the continental United States. Rainfall during this period is of low intensity and long duration. Beginning in early May, the tradewinds move north, again bringing with them the moist warm air of the Atlantic. The cyclonic activity is greatly reduced over the state and convectional instability begins to become established. June through September is known as the rainy season in Florida. The thunderstorms of this period are intense and very specially varied. They usually occur during the hottest part of the day and only on rare occasions last longer than two hours. About 60 percent of the total average annual rainfall occurs during this period, Figure 2. The mean average rainfall of Florida is in the neighborhood of 53 inches. It varies from 38 to 40 inches in the lower keys to over 65 inches in the southeast corner of the peninsula and the western portion of the panhandle. Most of the interior, that is the central highlands, receives approxi- mately the mean annual average. FIGURE 1.2 Mean annual precipitation (in inches). I Over 66 S62 66 S58 62 I 54 58 50 54 S46 50 LZ Less than 46 0 50 100 150 Scale of Miles ^- <- - Note: Precipitation normals compiled from records published in Climatological Data: Florida Section, December, 1971. _I~_ _ ~ Geology The Floridan peninsula and the offshore submerged lands above 50 fathoms, which Vaughan (1910) called the Floridian Plateau, have existed for several million years. The region has not been subject to violent earth movement and, consequently, there has been a gentle doing resulting in the formation of an oval arch above the basement rock. The rock of the core underlying the plateau is hypothesized to be pre-Cambrian; however, no drill has penetrated the core. The oldest rocks penetrated, to date, are quartzites found at about 4500 feet below the surface in Marion County. The borehole encountered another 1680 feet of quartzite before drilling was suspended. This metamorphized rock, believed to be a continuation of the Piedmont region of Georgia, was assigned by Cooke (1945) to the Pennsylvanian period. The arch above the metamorphized basement, composed of almost pure porous limestone, is known as the Lake City, Avon Park, and Ocala Formations. Dated in the Eocene period, the Ocala Formation has an estimated maximum thickness of 360 feet. It is found at, or above, mean sea level through- out northeast and north central Florida and is this section's principal aquifer. In southern Florida, in the vicinity of the Everglades, the Ocala is found at depths approaching 1300 feet. The Lake City and Avon Park limestones found below the Ocala are the principal aquifer used by agricultural interests in central and south central Florida and are known locally as the Floridan aquifer. Above the Eocene series are the formations of the Oligocene epoch. These are represented by the Marianna limestone and the Byran limestones found and mined in the Marianna Lowlands of the northern part of the state, and the Suwanneelimestone found over the Ocala Formation as far south as Hillsborough County. The next higher formations are those of the Miocene epoch. These are well represented by the Tampa limestones of the early Miocene which are found above the Suwannee and Ocala limestone in south Florida, the Chipola and Shoal River Formations of the Alum Bluff group found in northwest and north central Florida, the Hawthorn Formation, and the Duplin marls. The latter three formations, Alum Bluffs, Hawthorn, and Duplin, chiefly are sands, clays, and marls that form a confining layer over the Eocene and Oligocene limestones. The Hawthorn, with the possible exception of the Ocala, is the most extensive formation within the state. It occurs at, or near, the surface in most of north Florida. It overlies the Tampa limestone formation in Hillsborough County, and is, itself, overlain by the Duplin marls and younger deposits in the south central and southern parts of the state. 14 The surface material of most of the coastal lowlands are of the Pliocene, Pleistocene, and Recent periods. The most widely distributed are the sands formed along the old shorelines of previous ocean levels. Cooke (1945) defines seven of these marine terraces. Some small deposits of coquina, oolite, coral reef limestone, and freshwater marls are found among these deposits. CHAPTER II MODELS, NUMERICAL AND PHYSICAL The purpose of this chapter is to enumerate some of the more widely used modeling techniques in groundwater flow, along with a brief description of each. The reader is referred to Bear (1972) for additional information and references. A: NUMERICAL METHODS Numerical methods are used in many cases where the partial differential equations governing flow through porous media cannot be solved exactly. Various techniques have been developed for obtaining numerical solutions. Method of Finite Differences The method of finite differences is one such technique. The first step is to replace the differential equations by algebraic finite difference equations. These difference equations are relationships among values of the dependent variable at neighboring points of the applicable coordinate space. The resulting series of simultaneous equations is solved numerically and gives values of the dependent variables at a predetermined number of discrete or "grid" points throughout the region of investigation. If the exact solution of the difference equations is called D, the exact solution of the differential equation is called S, and the numerical solution of the difference equation is called N, two quantities of interest may be defined. They are the truncation error, IS DI, and the round-off error, ID NI. In order for the solution to converge, it is necessary that IS DI 0 everywhere in the solution domain. The stability requirement is such that ID NJ ->- 0 everywhere in the solution domain. The general problem is to find N so that IS NJ is smaller than some predetermined error. Noting that (S N) = (S D) + (D N), it is seen that the total error is composed of the truncation error and the round-off error. The arbitrary form selected for the finite difference equation leads to the truncation error. This error is frequently the major part of the total error. The actual computation proceeds by one of two schemes. They are the explicit, or forward-in-time, scheme and the implicit, or back-in-time, scheme. The explicit scheme is simpler but more time consuming than the implicit scheme due to the stability constraint. The implicit scheme is more efficient but requires a more complicated program as compared to the explicit scheme. Method of Finite Elements The finite element technique employes a functional associated with the partial differential equation, as opposed to the finite difference method which is based on a finite difference analog of the partial differential equation. A correspondence which assigns a real number to each function or curve belonging to some class is termed a functional. The calculus of variations is employed to minimize the partial differential equation under consideration. This is done by satisfying a set of associated equations called the Euler equations. Thus, one seeks the functional for which the governing equations are the Euler equations and proceeds to solve the minimization problem directly, rather than solving the differential equation. The procedure is continued by partitioning the flow field into elements, formulating the variational functional within each element and taking derivatives with respect to the dependent variables at all nodes of the elements. The equations of all the elements are then collected. The boundary condition is expressed in terms of nodal values and incorporated into the equations. The equations are then solved. Relaxation Methods This method may be applied to steady state problems which are adequately described by the Laplace or Poisson equations. The process involves obtaining steadily improved approximations of the solution of simultaneous algebraic difference equations. The first step of the procedure is to replace the continuous flow domain under investigation by a square or rectangular grid system. The governing differential equation is also replaced by corresponding difference equations. Next, a residual, say R is defined corre- sponding to the point o on the grid. Ro represents the amount by which the equation is in error at that point. If all values of the equation are correct, R0 will be zero everywhere. In the initial step, values are assigned at all grid points and, in general, the initial residuals will not be zero everywhere. The process now consists in adjusting values at each point so that eventually all residuals approach zero, or to at least some required accuracy. The reduction of residuals is achieved by a "relaxation pattern" which is repeated at different grid points so as to gradually spread the residuals and reduce their value. B: PHYSICAL MODELS As implied in section A, direct analytical solutions are frequently inadequate or impractical for engineering application. In many cases, the analytical solutions which are found are difficult to interpret in a physical context. In an attempt to circumvent some of the shortcomings of a purely mathematical approach, model and analog methods are frequently employed. The analog may be considered as a single purpose computer which has been designed and built for a given problem. Modeling, then, is the technique of reproducing the behavior of a phenomenon on a different and more convenient scale. In modeling, two systems are considered: the proto- type, or system under investigation, and the analog system. These systems are analogus if the characteristic equations describing their dynamic and kinematic behavior are similar in form. This occurs only if there is a one-to-one corre- spondence between elements of the two systems. A direct analogy is a relationship between two systems in which corresponding elements are related to each other in a similar manner. A model is an analog which has the same dimensions as the prototype, and in which every prototype element is reproduced, differing only in size. An analog is based on the analogy between systems belonging to entirely different physical categories. Similarity is recognized in an analog by two characteristics: (1) for each dependent variable and its derivatives in the equations describing one system, there corresponds a variable with corresponding derivatives in the second system's equations, and (2) inde- pendent variables and associated derivatives are related to each other in the same manner in the two sets of equations. The analogy stems from the fact that the characteristic equations in both systems represent the same principles of conservation and transport that govern physical phenomena. It is possible to develop analogs without referring to the mathematical formulation; an approach which is particularly advantageous when the mathematical expressions are exces- sively complicated or are unknown. Analogs may be classed as either discrete or continuous with respect to space variables. In both cases, time remains a continuous independent variable. The need for complete information concerning the flow field of a prototype system is obvious, and no method of solution can bypass this requirement. However, in many practical cases involving complicated geology and boundary conditions, it is usually sufficient to base the initial construction of the analog on available data and on rough estimates of missing data. The analog is then calibrated by reproducing in it the known past history of the prototype. This is done by adjusting various analog components until a satisfactory fit is obtained between the analog's response and the response actually observed in the prototype. Once the analog reproduces past history reliably, and within a required range of accuracy, it may be used to predict the prototype's response to planned future operations. Sandbox Model A reduced scale representation of a natural porous medium domain is known as a sandbox model, or a seepage tank model. Inasmuch as both prototype and model involve flow through porous media, it is a true model. A sandbox model is composed of a rigid, watertight container, a porous matrix filler (sand, glass beads, or crushed glass), one or several fluids, a fluid supply system and measuring devices. The box geometry corresponds to that of the investigated flow domain, the most common shapes being rectangular, radial, and columnar. For one-dimensional flow problems, the sand column is the most common experimental tool. Transparent material is preferred for the box construc- tion, especially when more than one liquid may be present and a dye tracer is to be used. Porosity and permeability variations in the prototype may be simulated by varying the corresponding properties of the material used as a porous matrix in the model according to the appropriate scaling rules. The porous matrix may be anisotropic. In order to measure piezometric heads and underpressures, piezometers and tensiometers may be inserted into the flow domain of the model. Wall effects are often eliminated by gluing sand grains to the walls of the box. This effect can also be reduced by making the porous matrix sufficiently large in the direction normal to the wall. Inlets and outlets in the walls connected to fixed level reservoirs or to pumps are used to simulate the proper boundary and initial conditions of the prototype. Water is usually used in models which simulate ground- water aquifers, although liquids of a higher viscosity may be used to achieve a more suitable time scale. The sandbox model is used extensively because of its special features which permit studies of phenomena related to the microscopic structure of the medium such as hydro- dynamic dispersion, unsaturated flow, miscible and immiscible displacement, simultaneous flow of two or more liquids at different relative saturations, fingering, wettability, and capillary pressure. The capillary fringe in a sandbox model is disproportionately larger than the corresponding capillary rise in the prototype and, for this reason, the sandbox model is usually used to simulate flow under confined rather than phreatic conditions. Hele-Shaw Analog The Hele-Shaw or viscous flow analog is based on the similarity between the differential equations governing two- dimensional, saturated flow in a porous medium and those describing the flow of a viscous liquid in a narrow space between two parallel planes. In practice, the planes are transparent plates, and the plates are usually mounted in a vertical or horizontal orientation. The vertical Hele-Shaw analog was selected for this study because it is more appropriate for the prototype system under investigation. Also, it is not possible to model a free goundwater table or percolation in a horizontal model. A detailed description of this analog is presented in Chapter III of this study. Electric Analog Three types of electric analogs are powerful tools in the study of flow through porous media. They are the continuous electric analog, the discrete electric analog and the ion motion analog. Continuous Electric Analog This analog takes two forms: the electrolytic tank and the conducting paper analogs. The analogy rests on the similarity between the differential equations that govern the flow of a homogeneous fluid through a porous medium, and those governing the flow of electricity through conducting materials. In particular, Darcy's law for flow in a porous medium and Ohm's law for the flow of an electric current in a conductor may be compared. Also, the continuity equation for an incompressible fluid flowing through a rigid porous medium may be compared with the equation for the steady flow of electricity in a conductor. One concludes from this comparison that any problem of steady flow of an incompressible fluid having a potential may be simulated by the flow of electric current in an analog. Discrete Electric Analog This analog also takes two forms: the resistance network analog for steady flow, and the resistance-capacitance network for unsteady flow. In this analog, electric circuit elements are concentrated in the network's node points to simulate the properties of portions of the continuous prototype field around them. The unknown potentials are the solution of the problem, and they can only be obtained for those points which correspond to the nodes of the analog network. The discrete electric analog is based on the finite-difference approximation of the equations to be solved; therefore, the errors involved in the discrete representation are the same as those occurring in this approximation. The electric resistor corresponds to the resistance of soil to flow through it, and capacitors are used at the nodes to simulate storage capacity of the prototype. Ion Motion Analog This analogy uses the fact that the velocity of ions in an electrolytic solution under the action of a DC voltage gradient is analogous to the average velocity of fluid particles under imposed potential gradients in a porous medium. In this case, both electric and elastic storativi- ties are neglected. The primary advantage of the ion motion analogy is that, in addition to the usual potential distri- bution, it permits a direct visual observation of the movement of an interface separating two immiscible fluids. In ground- water interface problems where gravity is involved, this analog cannot be used. Scaling for the analog is based on the similarity between Darcy's law and Ohm's law governing the ion motion in an electrolytic solution. Physically the analog consists of an electrolytic tank having the same geometry as the investigated flow domain. Inflow and outflow boundaries are simulated by positive and negative electrodes, and two- and three- dimensional flow domains may be investigated. Membrane Analog The membrane analog consists of a thin rubber sheet, stretched uniformly in all directions and clamped to a flat plane frame. The achievement of.equilibrium of various forces and stresses in the membrane (caused by distorting the frame or transversal loads) leads to the Laplace equation and to the Poisson equation. The analogy is based on the similarity between these two equations and the corresponding equations that describe the flow in the prototype. This method is applicable mainly to cases of steady two-dimensional flow involving complicated boundary geometry and point sources and sinks within the flow field. Summary Following Bear (1972), Table 2.1 is presented as a summary of the models and analogs discussed in section B of this chapter. In section A, the numerical methods discussed are most likely to be carried out on a digital computer. It is important for the investigator to examine both the cost and the applicability of these various numerical and physical methods to his particular case. An analog is usually pre- ferred to a digital solution when the accuracy and/or amount of field data is small. In many simple cases, the analog is likely to be less expensive than a digital computer; whereas, for large regions or unsteady three-dimensional problems, the computer may be less expensive. The Hele-Shaw model also has definite advantages when demonstration of the saltwater intrusion phenomenon to a public body, or other laymen involved in political decision making, is considered. This type of model allows for direct observation of the phenomenon without the numerical interpre- tations used in the computer models. TABLE 2.1 APPLICABILITY OF MODELS AND ANALOGS Characteristic Dimensions of field Steady or unsteady flow Sandbox Model Hele-Shaw Analog Vertical Horizontal two both Simulation of phreatic surface yes' Simulation of capillary yes fringe and capillary pressure Simulation of elastic yes, storage di Simulation of anisotropic yes media' Simulation of medium yes inhomogeneity Simulation of leaky formations yes Simulation of accretion yes Flow of two liquids with appr an abrupt interface Simultaneous flow of two yes immiscible fluids Hydrodynamic dispersion yes Observation of streamlines yes, dime tran for sion or three two two both both yes* no yes no Electric Analogs Electrolytic RC Network Ion Motion two or steady yes1 no three two or three both no' no for two yes yes yes, for two yes mensions dimensions yes yes yes yes kx5kz kx ky yes5 yes5 yes yes oximately yes yes yes' yes yes yes yes, for two yes dimensions yes yes no6 no6 no no r for two I nsions, near parent walls three dimen- s no no yes no two(horizontal) steady no no no yes kxk y yes no no yes no no no I Subject to restrictions because of the presence of a capillary fringe. z By trial and error for steady flow. I By trial and error for steady flow, or, as an approximation, for relatively small phreatic surface fluctuations. By scale distortion in all cases, except for the RC network and sometimes the Hele-Shaw analog where the hydraulic conductivity of the analog can be made anisotropic. s With certain constraints. 6 For a stationary interface by trial and error. Membrane Analog two(horizontal) steady no no no yes kx y yes no yes no no no no CHAPTER III THE HELE-SHAW MODEL The viscous flow analog, more commonly referred to as a parallel-plate or Hele-Shaw model, was first used by H. S. Hele-Shaw (1897, 1898, 1899) to demonstrate two- dimensional potential flow of fluid around a ship's hull and other variously shaped objects. The analog is based on the similarity of the differential equations which describe two-dimensional laminar flow, or potential flow for that matter, of a viscous fluid between two closely spaced parallel plates; and those equations which describe the field of flow below the phreatic surface of groundwater, namely Darcy's law: qx = K (3.la) x x ax qz = K (3.1b) where qx, qz = Darcy velocity of specific discharge in the x-direction and z-direction, respectively K K = hydraulic conductivity in the x-direction Kx z and z-direction, respectively x = horizontal direction (major flow direction) z = vertical direction = potentiometric head and, by use of the conservation of mass principle, the Laplace equation 2 + 822 = 0 (3.2) ax2 az2 Viscous Flow Analog To demonstrate the analogy of model and prototype, the equations of motion and continuity for laminar flow of a viscous fluid between two closely spaced parallel plates will be developed and then compared to equations 3.1a, 3.1b, and 3.2. Consider a viscous incompressible fluid flowing ever so slowly between two parallel plates which are spaced such that the Reynolds' number, R based on the interspace width is less than 500 (Aravin and Numerov, 1965). In Cartesian coordinates, the general Navier-Stokes equations, i.e., the equations of motion, are DV B + 1 I- + V2Vx (3.3a) DVy = B + + V2V (3.3b) Dt y+ y ] Vz = B + 1 3P + v2V (3 .3c) Dt z p DZ where D = substantial derivative =- + V a Dt at x ax + V a + V z - V2 = Laplace operator = 2 + y2 + Vx, V Vz = velocity in the x-, y-, and z-directions, respectively B B Bz = body forces in the x-, y-, and xz z-directions, respectively p = pressure p = density of the fluid p = absolute viscosity of the fluid y = horizontal direction (minor flow direction) t = time Referring to the free-body diagram of the idealized flow regime shown in Figure 3.1, if no slip conditions (adherence to the walls) of the fluid particles are assumed in the molecules closest to the walls of the parallel plates, it is easily seen that the velocity gradient in the y-direction is much larger than the velocity gradient in either the x- or z-directions. Thus., the first and second order partial derivatives taken with respect to both x and z may be neglected when compared to those taken in the y-direction. Secondly, because of the very low velocities ("creeping" motion) the intertia terms, that is the terms FIGURE 3.1 Free body flow diagram for Hele-Shaw model. b m _~I _ _II~ _I __ on the left side of equations 3.3, are very small when compared to the viscous terms, those on the right side of equations 3.3, and may be neglected. Thirdly, because of the restriction to two dimensions, the velocity in the y-direction is taken to be zero; consequently, all rates of change of velocity in the y-direction.must be zero. Finally, the only noncancelable body force acting on the fluid is gravity which acts only in the vertical. Mathe- matically, Bx -x (gz) = 0; By = (gz) = 0; and Dz -~ a (gz) = -g = -32.17 ft./sec.2. Incorporating all of the above arguments and values into equations 3.3, the equation of motion becomes: 32V 3+ 0 (3.4a) 2= 0 (3.4b) 82V pg a- = 0 (3.4c) Defining the potentiometric head, or potential, i = z + p/y, where y equals the unit weight of the fluid, and taking the partial derivative with respect to x, y, and z, the following results are obtained after multiplying through by the unit weight of the fluid: Y M =_ ap (3.5a) Y ax ax Y -- = DP (3.5b) Tay ay Y (y y 'I Introducing these relationships into equations 3.4 and dividing through by the unit weight yields: 32V 3t = P2 x (3.6a) 3x y ay2 D| 0 (3.6b) ay 32V S_ 1 2 z (3.6c) 2z y ay2 It is evident from equation 3.6b that the potentiometric head is constant in the y-direction. It is possible then to inte- grate the first and third equations of equation 3.6 with respect to y. After separating variables and integrating 3V once, using the boundary condition y = 0, -x = 0, and 3V ay z = 0, the following equations are obtained: ay y a 1 x (3.7a) a x y 8y y = z (3.7b) Yz y ay Integrating, once again, using the second boundary condition y = b/2, Vx = 0, Vz = 0 (no slip) the above becomes, after solving for the respective velocities: Vx y2 ) (3.8a) Vz 2 (2 ] (3.8b) Note, where b is the spacing between the plates, that if a potential = y- b-4 is defined, equations 3.8 can be written: V (3.9a) x ax V (3.9b) P, the velocity potential, is dependent only on y. Inte- grating the velocity profiles established by equations 3.8 between the limits of b/2, and dividing by b, the directional specific discharges are obtained: q +b/ 2 y= b/2 qx b _/2 x y b = x N -b/2 =- b2 Y (3.10a) f+b/2 qz b j-b/2 vz dy 2-i z (3 j-b/2 b2 Y (3.10b) TZ T az It is obvious that for a model of constant spacing b, the quantity b2 does not vary in either the x- or z-direction. Defining the model hydraulic conductivity as K = K = b* ~ equations 3.10a and 3.10b become: qx = Kxm x (3.11a) qz = Km (3.11b) which, of course, is analogous to equations 3.1. Consider, now, the two-dimensional continuity equation for flow between parallel plates: aV 3V _x + z 0 (3.12) ax az The specific discharge, or Darcy velocity, is related to the velocity by the vector equation n q = V, where ne is the effective porosity of the flow media. In the model, ne equals 1. From the analogy Vx = q x; Vz = qz' and substituting the relationship obtained from equation 3.11 into equation 3.12: x- m a-x x- Kzm = z (3.13) or dividing by -Km and recalling that for a model Km = Kxm = Kzm 9 + 2i-- (3.14) ax2 az2 (3.14) which is clearly analogous to equation 3.2. The similarity of equations 3.1 and 3.11 and equations 3.2 and 3.14 establish the analogy. Scaling The two-dimensional equation along the free surface, or water table, of an anisotropic porous media given by Bear (1972) is: K p .2 + K 2 = n (3.15) xp xp zp z az ep at Sj p p where the subscript p denotes the prototype. For a Hele- Shaw model, using the subscript m, the equation can be written as: 2 Mjt 2 K Im + K =- n mm (3.16) xm jx-m- zm azm a m em atm Introducing the similitude ratios, denoted by the subscript r, of the corresponding parameter of model and prototype: K Kxr = K (3.17a) xp K Kzr = (3.17b) zp x Xr -x (3.17c) P z zr = z (3.17d) P r -I (3.17e) P n ner = n (3.17f) ep t tr = t-- (3.17g) P and substituting these relationships into equation 3.15, the following is obtained: Kxm [( m ) 2 Kzm [*a(m/r) 2 (m r) T Tv K+ _zz Kxr m/x Kzr Zm/Zr) (Zm/Z nem < tm/tr) (3.18) The ratios of model to prototype quantities are constant and can be removed from behind the differential; therefore, equation 3.18 can be rearranged in the following way: Kxr r2 xm 2 Z r m 2 K r 2 xm ax m zm K ^ 2 3z Zr m'I tr m Kzr r -zmJ nerzr nem tm (3.19) Comparing the equations 3.16 and 3.19, it is evident that, if the equations are identical, the following must be true: x 2 z 2 z t 1 K r r r r (3.20) xr r zr r zr'r er r Solving the third equality for zr, the following important relationship is found: zr =r (3.21) The second equality, after cross-multiplying, yields: S xr K (3.22) zr Kzr Recalling the definitions of Kxr and Kzr (equations 3.17a and 3.17b) and remembering that K = K in an isotropic xm zm model, the above equation can be rewritten: x r 2 K/Kxp K r xm xp zp (3.23) r K zm zp Kxp K The ratio of zK is called the ratio or degree of anisotropy xp of the prototype. Time Using the fourth equality of equation 3.20, the time ratio of the model and prototype is established: n z t er r (3.24a) r K zr n x 2 t er r (3.24b) r K rr Substituting the vertical ratio, zr, for the potentiometric head established by equation 3.20 and the similitude ratios of time, hydraulic conductivity and porosity into equation 3.24b results in tm nem X r (3.25) tp nep Kxm z The effective porosity of an isotropic model, nem, is unity. The hydraulic conductivity of the model was defined pre- viously as p?, thus the time scale for the model is finally written as: K x 2 t 12 vK xx r t (3.26) m g ep zr P Anisotropy The Hele-Shaw model is normally isotropic. This is because of the nonvariance of the. spacing of the parallel plates. There are, however, two methods for simulating anisotropy in a model. Equation 3.23 gives a clue as to the first possibility of simulating anisotropy: r xm x_ zp (3.23) r zm /Kzp xp Since Km = Kzm, the x or z ratio can be adjusted so that the model's hydraulic conductivities are kept equal. This is usually done by choosing a suitable horizontal ratio. Know- ing the prototype parameters, a vertical scale for the model is computed so that the aforementioned conductivities are kept equal, demonstrating: zr zm X m (3.27) p zp p Solving for zm: K x zzm p xm z (3.28) zp p Unfortunately, the geometric distortion method is adequate for modeling only one ratio of anisotropy. If there is a second aquifer, within the prototype which has a different vertical or horizontal hydraulic conductivity, the second aquifer cannot be correctly simulated unless, of course, the second aquifer's ratio of anisotropy is the same as the ratio of the first. This restriction would severely limit the use of the Hele-Shaw analog in modeling of regional groundwater problems unless another method were available to correct the ratio of anisotropy. Polubarinova-Kochina (1962) suggests using a grooved plate within the model to correct the ratio of anisotropy of the second flow zone. The plate may be grooved in any one of several methods. It matters little whether a grooved plate is sandwiched between the parallel plates, or if rectangular bars are attached to the front or back plate. The degree of anisotropy of the second aquifer and the amount of geometric distortion used to model the first flow zone determines the directions in which the grooves, or bars, are placed; however, the grooves are normally placed horizontally or vertically. Collins and Gelhar (1970) have developed the conductivity equations for the flow zone in which Polubarinova-Kochina's grooved plate is used. The analysis assumes one-dimensional flow and can be used equally well with either vertical or horizontal orientation of the grooves. Following Collins and Gelhar (1970), consider flow in a grooved portion of a model. For simplicity, assume Figure 3.2 is a section of the grooved zone. Assuming such, the horizontal direction then corresponds to the x-direction and the grooves, which are vertical, lie in the z-direction. Area 1 is associated with the wider spacing of length ab. Area 2 is associated with the narrower spacing of length b. Since flow area 1 is the much larger of the two areas, most of the frictional head loss occurring through the total length L is developed in flow area 2 which has length (1 M)L. Lambda, X, is a length correction factor.. Referring to Figure 3.3, the potentiometric gradient ax across area 2 is: FIGURE 3.2 Section of anisotropic grooved zone in Hele-Shaw model. Back Plate Rectangular Bar AL Front Plate L (1 A)L Plan View Perspective View FIGURE 3.3 Head loss in grooved anisotropic zone. A 1-) x or z XL 1 --- -- __ ax (1 A)L (3.29) For high values of a: ___ A = 2C(3 .30) ax L L but, from equation 3.29, 2 = (1 A) 3 so that: a = (1 ) ~2 (3.31) ax ax Applying Darcy's law to area 2: q= K (3.32) qx x2 ax and, substituting the previous expression for :2- ax K S- x2 1 (3.33) x (1 A) ax The effective hydraulic conductivity in the x-direction then is: K Kx b2g (3.34) xm (1 X) 12v 1 ) (334) Consider vertical flow through the grooved zone illustrated in Figure 3.2. In particular, consider flow downward through areas 1 and 2. The total discharge through these areas can be written as the sum of the dis- charge through each, that is, Qzm = Q1 + Q2. Applying Darcy's law for the total discharge, Q: 48 Q, = K (abAL) (3.35) Q2 = K ( ) Lb (3.36) z2 az Adding Qi and Qz: Qm= K (aX) + K (1 )]bL (3.37) zm z1 z2z z For flow area 2 it is not unreasonable to assume that the frictional forces in the fluid boundary to either side of area 2 are negligible. Therefore, the vertical conduc- tivity in this area is the same as defined by the earlier analysis, that is: K b2 = b2 (3.38) Z2 = l- 12 where v is the kinematic viscosity. Furthermore, if b/ab << 1, the flow in area 1 can be assumed roughly equivalent to flow through a rectangular hole. According to Rouse (1959), the equation of motion through a rectangular cross section of length XL and width ab is given by: 00 v x(x XL) + X sin JTx z 2i 3z jl AL x A. cosh 1J + B. sinh (3.39) S1AL j AL 49 where A. 2y(AL) 2 (cos j I) (3.40) and, cosh jab 1( B. = A. (3.41) sinh jab3 XL By integrating Vz over the area and dividing by the total area ALab, the mean velocity is given by: V= (a b)2 (3.42) z az where S= 192 b j (cos (j ) 1)2 _ = 2 -L j( 4j5 2cab (3.43) Since the terms of the infinite series decrease as j5s, only the first term of the series need be considered and retained, so that: 1 19 tanh fXL (3.44) Additional information may be found in Appendix A. From equation 3.42, the equivalent hydraulic conductivity in area 1 is given by: K Y (ab)2 (3.45) Kz P 12 Finally, introducing the values found for K and K into zE Z2 equation 3.37: Qzm =-] (b) (aX) + b (1 U) bL (3.46) or, Qzm = (3 + 1 A) bL (3.47) from which it is seen that the effective vertical hydraulic conductivity is given by: Kzm = (a 3X1 + 1 X) (3.48) after defining Qzm/bL = Vz, where Vz = qz is the effective vertical specific discharge. Equations 3.34 and 3.48 give the second method available to correct the hydraulic con- ductivity of a model so that it can simulate the true ratios of anisotropy found in the prototype. Leakage An aquiclude can be defined as a soil stratification in which the hydraulic conductivities are zero. In certain geohydrologic problems, it is convenient to assume such conditions. However, in reality few soil masses are truly impervious. The degree of perviousness in a stratum is referred to as leakance and it is generally assumed that the direction of flow is only vertical. There is no horizontal flow, that is, Kxp = 0 (3.49) Bear et al. (1968) suggest the use of vertical slots to model such a semipervious layer. To accomplish this, the spacing between the parallel plates of the Hele-Shaw analog is filled with a slotted middle plate. See Figure 3.4. The analysis to determine the effective vertical hydraulic conductivity of a model's leaky layer closely parallels that for an anisotropic grooved zone. Again, following Collins (1970), Darcy's law for flow through a vertical slot is: QZ = ALab (3.50) z = z 3z The effective specific discharge through the slot found by integrating the Rouse equation (equation 3.39) is the same as equation 3.42 from which is found the hydraulic conduc- tivity: Kz (b)2 (3.51) and, introducing the above into equation 3.50: Qz = (b)2 Lab (3.52) z U -2-- aL z FIGURE 3.4 Section of leaky zone in Hele-Shaw model. y Back Plate Flowspace Solid ab x l b F-I v PL (1-X)L t Front Plate L Plan View Perspective View Again, the effective specific discharge, or mean velocity, is equal to: z 3 b (3.53) Vz bL 12 az 353) so that the effective hydraulic conductivity of a leaky layer in the model is: Kzm a b2 (3.54) Storativity While the problem of storage has not been completely solved, it has, in general, been neglected by most researchers. Bear (1960) suggests that discrete tubes attached to either the front or back plate and connected to the aquifer be used to model the specific storage of a confined aquifer. For a nonisotropic aquifer, the right hand side of equation 3.13 is not zero, but, in fact, equals the specific storage, S , times the rate of change of the potentiometric head. Rewrit- ing the two-dimensional equation 3.13 for both model and prototype to include the above gives the following: K -2 + K --- = S -- (3.55) xp ax 2 zp aZ 2 op at p p p K 2m + K m m (3.56) xm Dxm2 zm 3Z%2 om atm m m m Defining a ratio of storativity: S Sor = (3.57) op it follows from inspection that, z K z K Sr- -- (3.58) xr x2 zr or t or that, K t S zr r om (359) or z S (3.59) r op Referring to Figure 3.5, the storage represented by the model in the discrete length 1m is equal to: A S om_ m (3.60) mo m b 1 z mmm where A is the cross-sectional area of the storativity tube. Introducing the above into equation 3.59 and solving for Am t A bl z Sop K r (3.61) m mmm op zr z 2 r Discharge The discharge scales are.obtained from Darcy's law. Written for both prototype and model with the usual sub- scripting, these are in the x-direction: Q K -- b z (3.62) xp = xp ax pp FIGURE 3.5 Storage manifold for Hele-Shaw model. Area, A bm m Front Plate Back Plate End Section I I Elevation and, (3.63) xm -xm 3x mm m Dividing equation 3.63 by equation 3.62 and recalling the definitions for the various parameters' ratios, it follows that: z z 2 S K r b r=K b r xr xr r r x xr rx Similarly, in the z-direction, x Q = K r b r K b x zr zr r r zr zr r r r (3.64) (3.65) Solving equation 3.22 for the hydraulic conductivity in the x-direction: Kx K xr zr S2 Xr r (3.66) and, substituting that: Qxr = Kzr this result into equation 3.64, it follows r b r K b x r r or, Qxr = Qzr = Q *xr -zr -r (3.67) (3.68) Accretion Accretion, R, is the rate at which a net quantity (precipitation and surface inflow minus evapotranspiration, runoff, etc.) of liquid is taken into the flow system at the phreatic surface. It is measured as a volume per unit horizontal area per unit time, that is: Q Rr (3.69) rr From equations 3.64 or 3.65, it follows that: K z 2 R x Kzr (3.70) r Volume On occasion, volume, U, is of some importance. The volume scale follows directly from continuity, that is: Ur = Q tr (3.71) Substituting the values found from equations 3.65 and 3.24a for Qr and tr, respectively, the above equation becomes: z r zr r r er K r er r r (3.72) zr As inferred earlier in this section's opening sentence, the volume scale is usually neglected; however, in the case of free surface water bodies, lakes, rivers, etc., if the volume exchange of liquid is of interest and has to be modeled, the 60 volume scale requires an additional restriction. In the following analysis, the bar above the width dimension indicates the free water surface of a river, lake, ocean, or such. In the portion of the model simulating the body of water, the spacing of the model is increased to maintain hydrostatic pressure distributions within the model. The narrower spacing of the model is, of course, a measure of the hydraulic conductivity of the aquifer. In the proto- type, however, the width of the open water and the aquifer are equal and this leads to the following (Bear, 1960) for the model and prototype, respectively: U = n b x z (3.73a) r er r r r U = n b x z (3.73b) r er r r r The same volume ratio must be applicable to both the narrow and the enlarged interspace; therefore, Ur Ur. It follows that: n er b = ner b (3.74) but, n n em 1 (3.75) er n ep so, br = ner br (3.76) Note that for an anisotropic media, nem does not necessarily equal one. CHAPTER IV SITE SELECTION AND PROTOTYPE GEOLOGY AND HYDROLOGY Site Selection The site selected for this study is the middle Gulf area of Florida. This region has a rapidly expanding popu- lation with a corresponding growth in water demand. The increased pumping to satisfy this demand also increases the likelihood of saltwater intrusion, and, in fact, a number of municipal supply wells in the coastal zone have been shut down in recent and past years due to chloride contamination. Black et al. (1953) list eight factors responsible for saltwater intrusion. They are: 1. Increased water demands by municipalities. 2. Increased water demands by agriculture. 3. Increased water demands by industry. 4. Excessive drainage. 5. Lack of protective works against tidewater in bayous, canals, and rivers. 6. Improper location of wells. 7. Highly variable annual rainfall with insufficient surface storage during droughts. 8. Uncapped wells and leakage. Of these eight factors, numbers 1, 2, 3, 6, and 7 would apply to this area. The city of St. Petersburg is an outstanding example of these problems. Their original water supply was from local artesian wells, but increased demands caused salt- water intrusion and forced the closing of these wells. In 1929 the present Cosme-Odessa field was located farther inland to escape this problem. One of the major water supply systems in this region is the Pinellas County Water System, and this study is centered around the Eldridge-Wilde well field of this system. The location of Eldridge-Wilde in relation to several of the population centers of this region is shown in Figure 4.1. It is about 8 miles east of the Gulf of Mexico and encom- passes an area in the northeast corner of Pinellas County, at the intersection of the boundaries of Pinellas, Hills- borough, and Pasco counties. This system was instituted in 1937 to supply the towns along the Gulf coast from Belleair Beach to Pass-a- Grille. Its original form was that of a raw water reservoir, and the first wells were not drilled until 1946 in the McKay Creek area. These wells were soon contaminated with salt water, and investigations were begun in 1951 to locate the well field at its present site. This well field has grown over the years, and in 1970 (Black, Crow and Eidsness, Inc., 1970), the waterworks facilities at Eldridge-Wilde included: sixty-one water wells, over 11 miles of raw water FIGURE 4.1 Regional area of prototype. il MEXIUI0GQ;:1;:.: , iiii -- / ko V -,^ / 10' . .. Dade City 60' / 40' , \ 70' 0 S80' ^ ./ I K N "- 50' S 30' Eldridge-Wilde Well Field Lake Tarpon 0O Tampa :: :1 :i:: Peter .sbug : ::? :::: ::::::l ::.. : : ... .. : :::::: ::::: : ...................i 9 collection piping, water treatment facilities consisting of aeration and chemical treatment, including chlorination and fluoridation, and high service pumping units. All wells are open hole and penetrate the Floridan aquifer at depths from 140 to 809 feet below ground surface, averaging 354 feet. The design capacity, of the field at the present time is 69 million gallons per day, although the maximum allowable pumpage has been set by the Southwest Florida Water Management.District at 28 million gallons per day on the average with a maximum day of 44 million gallons per day. In selecting the prototype location within the site area, two characteristics of the vertical Hele-Shaw analog must be considered. The first characteristic is that there can be no general flow normal to the parallel walls of the model. This means that the flow from one end of the model to the other is streamline flow. The second characteristic is that the ends of the model are finite. Therefore, the prototype must be along a streamline in the flow domain and have boundary conditions which are "infinite" reservoirs or water divides. The prototype selected meets the above requirements and includes the point of interest, i.e., Eldridge-Wilde well field. The center line of the prototype is shown in Figure 4.1 as the unbroken line passing through Eldridge- Wilde in a southwest to northeast direction. The dotted contours in the figure define the potentiometric surface of the Floridan aquifer in feet above mean sea level as of May, 1971. They were obtained from a map publication entitled "Potentiometric Surface of Floridan Aquifer Southwest Florida Water Management District, May, 1971" prepared by the U. S. Geological Survey in cooperation with the Southwest Florida Water Management District and the Bureau of Geology, Florida Department of Natural Resources. Now, in a flow field, streamlines are perpendicular to potentiometric lines. As can be seen from the figure, the prototype orientation reasonably satisfies the streamline requirement. The proto- type is terminated on the southwestern end at the 15 feet depth contour in the Gulf of Mexico, and it is assumed that this satisfies the infinite reservoir boundary condition. The northeastern terminus is located in the center of the 80 feet contour, southwest of Dade City. This location satisfies the water divide boundary condition. The area in the vicinity of the 80 feet contour is known as the Pasco High. The overall length of the prototype is 36 statute miles. The width of the prototype is taken to be 3.5 statute miles. This dimension is sufficient to include the cone of depression caused by pumping in Eldridge-Wilde well field, and is based on the results of a study by Mr. Evans (employing a numerical model) for Black, Crow and Eidsness, Inc. The land surface contours of the prototype were obtained from U. S. Geological Survey topographic maps. The bottom boundary of the prototype is taken to be the base of the Lake City Formation, with depths being determined from available well logs of wells in the prototype vicinity. The maximum depth from highest land surface to deepest point is 1340 feet. Prototype Geology Stewart (1968) identifies eight formationsas being of interest in terms of water production in the prototype area. They are in descending order, the Undifferentiated Deposits, Tampa Limestone, Suwannee Limestone, Crystal River Formation, Williston Formation, Inglis Formation, Avon Park Limestone, and Lake City Limestone. Underlying the Lake City Limestone is the Oldsmar Limestone which is not used as a source of water at present. The Undifferentiated Deposits are interbedded sand, silt, and clay of Post-Miocene age and range in thickness from zero near the Pasco High to 60 feet in the Eldridge- Wilde well field. The thickest deposits are in northeast Pinellas County around the north end of Lake Tarpon where sand dunes, as much as 40 feet high, overlie alternating layers of clay, thin limestone beds, and sand greater than 70 feet. thick. The Tampa Limestone is a hard, dense, sandy, white to light tan, or yellowish-tan fossiliferous limestone of Miocene age. This limestone is near the surface in the area of the Pasco High and about 80 feet below land surface at the Eldridge-Wilde well field. At Eldridge-Wilde, the thickness varies erractically from about 20 to 240 feet. The Tampa Limestone is a poor to fair producer of water. The Suwannee Limestone is a soft to hard, nodular or grandular, fossiliferous white to tan limestone of Oligocene age and is about 200 feet thick. The Suwannee and Tampa Limestones are the major water producers for wells in the area. The Crystal River, Williston, and Inglis Formations comprise the Ocala Group of late Eocene age. The Crystal River and Williston Formation are lithologically similar units of white to cream, porous, soft, coquinoid limestone and are generally poor producers of water. The Inglis Formation is a hard, cream to brown to gray fossiliferous limestone and is generally a good producer of water. The Avon Park and Lake City Limestones are litho- logically similar units of soft to hard, cream to brown, fossiliferous limestone with beds of dolomitic limestone and some gypsum. Both formations are good producers of poor quality water. The Oldsmar Limestone is a fragmental dolomitic limestone with lenses of chert, thin shale beds, and some gypsum. In this study, two formations are considered, the Undifferentiated Deposits and the Floridan aquifer. The Floridan aquifer is considered to contain all formations from the Tampa to, and including, the Lake City Limestone. The transmissivity of the Floridan aquifer ranges from about 165,000 to 550,000 gallons per day per foot, and the coefficient of storage ranges from about 0.0005 to 0.0015. The coefficient of leakage is approximately 0.0015 gallons per day per cubic foot. Based on groundwater discharge and water levels, the estimated recharge (leakance) to the Floridan aquifer was computed (Stewart, 1968) to be about 103 million gallons per day. Based on aquifer test data, the estimated recharge for a 250 square mile area was 90 million gallons per day. The Undifferentiated Deposits act as a confining layer, and the Floridan aquifer is thus under artesian conditions. Prototype Hydrology The surface waters of the area consist of many lakes and few streams. Because of the flat topography, little water runs off into streams, and swampy wetlands are numerous. Most rainfall evaporates or is transpired by plants. The Floridan aquifer is recharged through the Undifferentiated Deposits by surface and groundwater derived from local rainfall. Many millions of gallons of water are also admitted to the aquifer by numerous sink- holes in the region. Water levels in the Floridan aquifer respond to rainfall since this is the recharge source. This response is not immediate, but usually fluctuates with the wet and dry seasons. Water levels in wells which are not directly affected by local pumping show yearly lows in the dry season, April and May, and yearly highs during the wet season, late summer or early fall (Black, Crow and Eidsness, Inc., 1970). The aquifer recharge has been estimated (Black, Crow and Eidsness, Inc., 1970) from available data and the use of the following formula: Aquifer Recharge = P + SWI + GWI ET R GWO The basin area is 575 square miles, changes in storage are assumed zero, and evapotranspiration is assumed to be 75 percent of the precipitation. The applicable values are listed below in million gallons per day: P = Precipitation SWI = Surface Water Inflow GWI = Groundwater Inflow ET = Evapotranspiration R = Runoff GWO = Groundwater Outflow Aquifer Recharge This value is in reasonable reported values (Stewart, 1968). = +1492 = + 0 = + 0 = -1119 = 218 = 37 = + 118 agreement with previous CHAPTER V DESIGN, CONSTRUCTION, AND OPERATION OF MODEL Design Prototype The selection of the prototype area was discussed in Chapter IV. Table 5.1 is a summary of the prototype characteristics. The leaky layer is synonymous with the undifferentiated deposits. The top and bottom of the Floridan aquifer were determined by straight-line extrapo- lation from available well logs. The hydraulic data are within the reported range of values and are the result of a trial and error process to stay within the range and still-produce a reasonable model. Model The purpose of the Hele-Shaw analog in this study is to model saltwater intrusion." Before discussing the model design, it seems appropriate at this point to pro- vide some background about saltwater intrusion. Water, in general, whether it be surface water or groundwater, is continually migrating towards the sea, where an equilibrium, or moving freshwater/saltwater interface, TABLE 5.1 PROTOTYPE PARAMETERS Floridan Parameters Aquifer Leaky Layer Geometric x (ft.) z (ft.) yp = bp (ft.) Hydraulic Txp (gpd/ft.) Tzp (gpd/ft.) Kxp /Kzp Leakance (gpd/ft.3) Kzp (gpd/ft.2) Kxp (gpd/ft.2) S v @ 77* F (ft.2/sec.) gp = 32.2 ft./sec.2 190,080 1,340 18,480 225,000 184,426 1.218 161.4 196.7 0.00158 0.965 (10-5) 190,080 55 18,480 0 0.0015 0.09 0 0.965 (10-5) is established. The two fluids are miscible, but because of the difference in densities and the very low velocities, the interface is formed. Across the interface, the salinity varies from that of the fresh groundwater to that of the ocean. The transition zone, as it is called, is due to hydrodynamic dispersion and, although it is anything but abrupt, it is usually assumed to be. The interface then is generally selected to occur at some measured electric conductivity or salt (chloride) concentration. The earliest investigations of saltwater encroch- ment were made by Badon-Ghyben (1888) in Holland and Herzberg (1901) in Germany. Working independently, both investigated the equilibrium relationships between the shape and position of the freshwater/saltwater interface. Figure 5.1 shows a coastal phreatic aquifer and the Ghyben-Herzberg interface model. Badon-Ghyben and Herz- berg assumed static equilibrium and a hydrostatic pressure distribution in the fresh groundwater and stationary saline groundwater near the interface. Considering a point P on the interface, and choosing mean sea level as the datum, the pressure at point P is: pp = h 5 (5.1) where h = vertical distance from mean sea level to s point P s = unit weight of sea water FIGURE 5.1 Ghyben-Herzberg interface model. 75 Water Table Vertical h h pdcA Interface Y^ \pdA Y This pressure may, also, be expressed by: p = {hf + hs] f (5.2) in which hf equals the vertical distance from mean sea level to the phreatic line at the location of p and yf equals the unit weight of freshwater. Equating the, preceding two equations: hSys = hfyf + hsf (5.3) and, rearranging, hsYs hf = hff (5.4) Solving for hs : yf hs Y= Yf hf (5.5) Introducing y = pg, where p is the density, factoring and canceling out the gravity term, the Ghyben-Herzberg relation is found: f h = -f hf (5.6) S P-S f for a saltwater density of 1.981 lb sec.2/ft.4 and a fresh- water density of 1.933 lb sec.2/ft.4, @ 250 C, the quantity pf/ps pf) = 40. The implications of equation 5.6 are rather dramatic. For instance for every foot of freshwater above the datum, there is 40 feet of freshwater below the datum. More importantly however, consider the effects of lowering the phreatic surface. For every one foot drop of the water table, the interface raises 40 feet. It must be remembered that the above analysis assumes static conditions. This, in fact, is not always the case. The position of the interface is a function of dynamic conditions rather than static. Even so, in cases where flow is quasi- horizontal, i.e., the equipotential lines are nearly vertical, equation 5.6 is valid. Many investigators have incorporated dynamic forces into the analysis of the stationary interface. Hubbert (1940) was able to ascertain a more accurate determination of the shape of the interface near the coast line. He assumed that at the interface the tangential velocity was zero in the saltwater, but increases with horizontal dis- tance in the freshwater as the coast line is approached. This then is the cause for the interface to tilt upwards as the sea is approached and the greater depths found than those estimated by the Ghyben-Herzberg relationship. Hubbert showed that the Ghyben-Herzberg equation holds between points on the water table-and the interface along an equipotential line, rather than along a vertical plane. R. E. Glover (1959) modeled an infinitely deep coastal aquifer by assuming no flow in the saltwater region, a horizontal water table and a horizontal seepage face located seaward of the coast line. He found an exact solution for the shape of the wedge, giving the following relationship: z2 2qx + Pf 2 (5.7) K s f K2 S f- Pf Pf where x and z are the horizontal and vertical directions, respectively, q is the seepage rate per unit width and K is the hydraulic conductivity. De Wiest (1962) using complex variables and a velocity potential of D = Kx . p Pf)/pf derived the same equation. Bear and Dagan (1964b), using the Dupuit assumptions and the Ghyben-Herzberg equation, developed the approximate shape of the interface for a shallow aquifer of constant depth. All of the above investigated the equilibrium position of the saltwater/freshwater interface. If there is a change in the freshwater flow regime, a transition period is caused during which the interface moves to a new point of equilibrium. The nonlinear boundary conditions along the interface make the solution for the shape and position of the transient interface all but impossible except for the simplest geometries. Bear and Dagan (1964a), as well as other investigators, have used the Dupuit assumptions to approximate the rate of movement of an interface in a confined aquifer. Following Polubarinova-Kochina's (1962) suggestion, they assumed quasi-steady flow and were able to approximate the inter- face shape and position for both a receding motion and landward motion of the interface. Characteristically, the solutions obtained by investigators to date have all had simple,geometries and involved simplifying assumptions, some of which have had little resemblance to actual conditions. Therein lies the advantages of a Hele-Shaw analog, complex geometries and boundary conditions can be modeled with relative ease. In order to satisfy additional similitude require- ments for the flow of two liquids with an abrupt interface, as in this study, and to provide a suitable time ratio, two liquid silicone fluids were chosen to be used in the model. Dow Corning Corporation Series 200 silicone fluid was used to model fresh water. Series 510 silicone fluid from the same company was used to simulate salt water. The 200 Series fluid and the 510 Series fluid have densi- ties of 0.977 gm/cm3 and 1.00 gm/cm3, respectively, at 250 C. Both fluids have kinematic viscosities of 500 centistokes at 250 C. Dow Corning 200 fluid is a clear dimethyl siloxane which is characterized by oxidation resistance, a relatively flat viscosity-temperature slope and low vapor pressure. Dow Corning 510 fluid is a clear phenylmethyl polysiloxane which also has a relatively flat viscosity-temperature slope. In order to locate and follow the interface movement, the denser fluid was dyed blue and the lighter fluid dyed orange. The dimensions of the model were selected so that a unit of reasonable size would be produced. These dimensions are 11.75 feet long, 0.5 inch wide inside, and 2 feet deep. The parameters xr, Zr, and br are therefore set, and the result is a distorted model. This distortion requires the use of a slatted inner zone as discussed in Chapter III. Tables 5.2 and 5.3 list the applicable parameters. The following analysis is shown for the design of the model Floridan aquifer, leaky layer and storage coefficient: A. Slatted Anisotropic Zone for the Floridan Aquifer. 1. Compute k xm/kzm from equation 3.22 x 2 Kxr r _= xr (3.22) Sr Kzr Noting that k = K v (5.8) g Then 3.22 becomes, x '2 k k zr XM Z kz (5.9) r ixp zm and finally, p [ -6.182 10-0|2 xMm = (1.22) 0 [1 = 0.00209 zm Kzp zr1.493 x 1 TABLE 5.2 MODEL PARAMETERS Floridan Parameters Aquifer Leaky Layer Geometric x (ft.) z (ft.) ym (ft.) a A b (in.) L (in.) ab (in.) XL (in.) (1 A)L (in.) Hydraulic kxm (in.2) kzm (in.2) xm /kzm S (ft.-1) vm @ 77* F (ft.2/sec.) gm = 32.2 ft./sec.2 11.75 2.0 0.0417 16.50 0.8097 0.030 1.235 0.6905 0.495 1.00 0.235 3.94. x 10-4 .1884 .00209 0.0369 5.382 x 10-3 11.75 0.082 0.0417 0.25 0.3846 0.500 1.235 0.8349 0.125 0.475 0.760 0 1.05 x 10-4 0 5.382 x 10-3 82 TABLE 5.3 SIMILARITY RATIOS This is the value which the slatted zone must produce. 2. Select dimensions for the slatted zone as shown in Figure 3.2, as follows: hold ab = 0.5 inch select b find a select (1 X)L select XL find L find x 3. Compute kxm b2 ( 1 (3.34) 4. Compute kzm = ((3AE + 1 A) (3.48) 5. Compute kxm/kzm and compare to the results of step (1). 6. Repeat the process until the result of step (5) equals the result of step (1). B. Slotted Leaky Layer. Refer to Figure 3.4 for this section. kxlm = 0 (5.10) b2 (3.54) kzlm = i-Z a (3.54) 84 Also, continuity of flow between the leaky layer and the Floridan aquifer requires that: (5.11) Qlr = Qfr Also, Q = Qzr = Qxr = K bzrb x "r zr "xr zr rr Therefore, K b x =K b x zlr rr zfr rr so, Kzlr = Kzfr or, kzlr = kzfr (3.67) & (3.68) (5.12) (5.13) (5.14) Therefore, kzlm = kzlp kzfr Now rewrite equation 3.54 as, a3 k 12 zlm U_ 1. Set b = 0.5 inch. 2. a3A = constant, since b is set and kzlm can be computed from previous information. (5.15) 3. Select ab and find a. 4. Select XL and L and find X. 5. Compute E. 6. Compute a3A. 7. Compare the result of step (6) to the result of step (2). 8. Repeat the process until the result of step (6) equals the result of step (2). 9. Make sure the physical size of this layer is compatible to the slatted layer, especially in terms of slot spacing, i.e., blockage. C. Storage Coefficient Manifolds. Refer to Figure 3.5 for this section (bm = 0.5"). 1. Take average storage coefficient and average depth in prototype to compute S 2. Select convenient time ratios, in this case 1 minute = 1 year, and compute Sor from equation 3.59. 3. Compute S om 4. Equation 3.60 is now used to compute A for m various 1 's with z averaged over 1 . m m m In this case the model was apportioned into five zones with one manifold per zone. Construction Because the Hele-Shaw analog is capable of modeling complex geometries and boundary conditions, it is desirable that it be as adaptable to as many different prototype geometries and hydraulic parameters as possible. This would facilitate model construction and investigations of many different areas in the state of Florida where saltwater intrusion is, or in the future might be, a problem. A reduction in cost of investigation would also be achieved if many of the parts were reusable. A list of general specifications would then be as follows: 1. The Hele-Shaw model should be housed in a frame in which it can be easily installed and removed. 2. The front and back plates with interior model parts should not be permanently sealed together. 3. The front and back plates should be as adaptable as possible to different situations. 4. The model should have as few opaque parts as possible. 5. The model should be mobile. Frame As shown in Figure 5.2, the frame is composed of two assemblies; the cradle and the cradle dolly. The function of the cradle is to support and orient the Plexiglas plates. It also contains the inflatable neoprene hose which seals the plates. It is fabricated of 2-1/2" x 2" x 3/8" steel angles which are welded into a channel |