Two- DIMENSIONAL FINITE ELEMENT ANALYSIS DEV. TRANSIENT FLOW AND TRACER NDYENIENT IN CONFINED. 2, i=1,- 5 , ‘ ’9‘ ' AND PHREATIC AQUIFERS ‘ ~ DIssertatIon for the Degree of Ph D MICHIGAN STATE UNIVERSITY SIROUS HAJI- DJAFARI 1976 TSHH}S§T§§:}7’TI1;gfigf‘girigtgj‘g‘gz'azgakhn':vii-7‘31. I.. - ,_,- , 1‘.41.I‘..',‘1\-.:. .'.=1-; “If“ _A _ .w, . . _ H V ..... 12mm AH .. 9 m IIIIIIIIININIIIIIIIIIINNNNIIIIIIIIIIIIIII W 4 Unifvm. % TWO-DIMENSIONAL FINITE ELEMENT ANALYSIS OF TRANSIENT FLOW AND TRACER MOVEMENT IN CONFINED AND PHREATIC AQUIFERS presented by Sirous Haji-Djafari has been accepted towards fulfillment of the requirements for Ph.D. degree in Civil & Sanitary Engineering Major 0 ssor /> C fl/CZ/flQZS/ p{ r Date W6 0-7639 IN” ABSTRACT TWO-DIMENSIONAL FINITE ELEMENT ANALYSIS OF TRANSIENT FLOW AND TRACER MOVEMENT IN CONFINED AND PHREATIC AQUIFERS BY Sirous Haji—Djafari In this study the movement of a tracer in an aquifer with transient flow conditions is investigated, both on a regional as well as local scale. For the regional scale the two-dimensional horizontal plane is considered, while for the local scale a vertical cross section of a site is chosen. Special emphasis is placed upon solving the flow and mass transfer phenomena in a phreatic aquifer with a time variable boundary. Finite element formulation of the flow and convective-dispersion equations leads to a set of first order partial differential relations. In addition, with use of the finite element concept, higher order time approximations for the system of equations are derived. In order to obtain continuous flow across elements and at the nodes, the Galerkin formulation of the Darcy law is constructed and velocity vectors are calculated simul— taneously at the nodes. These transient velocities are Sirous Hafi-Djafari subsequently used in shifting the phreatic surface and in computing dispersion coefficients and convective terms of the mass-transport equation. A procedure is adapted which locates the phreatic boundary of an aquifer without repositioning the nodal coordinates of the elements. The validity of the proposed techniques is estab- lished by first comparing the numerical flow results with existing analytical, experimental, and field data. Upon verification of the solution of the flow equation, the prediction of the movement of a tracer in an unconfined aquifer with a time—variable phreatic boundary, and in a confined aquifer with a transient flow condition, is con- ducted. Numerical examples are presented to demonstrate the capability of the proposed techniques. It is shown that the Galerkin finite element method can be used to solve the flow and convective— dispersion equations, both in confined and phreatic aquifers under time-dependent flow situations. TWO-DIMENSIONAL FINITE ELEMENT ANALYSIS OF TRANSIENT FLOW AND TRACER MOVEMENT IN CONFINED AND PHREATIC AQUIFERS BY Sirous Haji-Djafari A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Civil and Sanitary Engineering 1976 To my mother and in memory of my father ii ACKNOWLEDGMENTS Sincere appreciation is expressed to the wri— ter's major professor, Dr. David C. Wiggert, for his encouragement and guidance throughout this study. Dis- cussions with the writer's committee members, Drs. J. C. Beck, J. F. Foss, and L. J. Segerlind, were very inspir- ing and constructive, and their cooperation and assistance gratefully appreciated. Acknowledgment is also extended to Drs. J. Bear, P. W. France, G. F. Pinder, and J. T. Oden for their technical contributions. Gratitude is expressed to the Division of Engineering Research for fellowship funds. The investigation was supported in Ipart by funds provided by the United States Department of the Interior, Office of Water Research Technology, as authorized under the Water Resources Research Act of 1964, and by a grant from the Rockefeller Foundation to the Water Quality Project, Institute of Water Research, Michigan State University. Special thanks are due to Mrs. K. C. Cornelius for her effort in typing the manu- script and careful drafting of most of the figures. The writer is grateful to his wife, Azizeh, and sons, Valla and Sina, for their sacrifice, understanding, and encouragement throughout the pursuit of this thesis. TABLE OF CONTENTS LIST OF TABLES . . . . . . . . . . . LIST OF FIGURES . . . . . . . . . . . NOMENCLATURE . . . . . . . . . . . . Chapter I. INTRODUCTION . . . . . . . . . 1.1 Motivation . . . . . . 1.2 Objective and Plan . . . . . l 3 Scope of Study . . . . . . . II. LITERATURE REVIEW . . . . . . . III. 2.1 Introduction . . . . . . . 2.2 Mechanics of Flow . . . . 2.2.1 Confined Aquifer . . . 2.2.2 Unconfined Aquifer With Phreatic Surface . . 2.3 Mechanics of Convective- Dispersion Phenomena . . . . . . . . MATHEMATICAL REPRESENTATION OF PHYSICAL SYSTEM . . . . . 3.1 Background . . . . . . . . 3.2 Darcy Equation . . . . 3.3 Regional Groundwater Flow . . . 3. 3.1 Basic Assumptions . 3.3.2 Two- —Dimensional Horizontal Flow . . 3.3.3 Initial and Boundary Conditions . . 3.4 Phreatic Aquifer With Accretion . 3.4.1 Basic Assumptions 3.4.2 Governing Differential Equations . 3.4.3 Boundary and Initial Conditions . . . . iv Page viii ix DUNE-4 mmfl 12 l7 17 21 21 22 22 26 28 29 Chapter Page 3.5 Convective-Dispersion Phenomena . . 31 3.5.1 Basic Assumptions . . . . 32 3.5.2 Hydrodynamic Dispersion Coefficient in a Porous Medium . . . . 33 3.5.3 Convective- Dispersion Equation in Cartesian Coordinates . . 35 3.5.4 Initial and Boundary Condi- tions . . . . 36 3.6 Closure . . . . . . . . . . 38 IV. FINITE ELEMENT FORMULATION . . . . . 39 4.1 The Galerkin Finite Element Method . 40 4.2 Finite Element Formulation of Flow Equation . . . . . 42 4.3 Finite Element Formulation of Convective- Dispersion Equation . . 46 4.4 Finite Element Computation of Velocity Vectors . . . . . 50 4.4.1 Direct Calculation . . . 51 4.4.2 Simultaneous Calculation of Velocity Vectors . . . . . 53 V. FINITE ELEMENT FORMULATION OF TIME DERIVATIVE _ , , _ , _ , . . . . 59 5.1 Introduction . . . . . 59 5.2 Finite Element Formulation . . . . 62 5.3 First Order Approximation . . . . 62 5.4 Second Order Approximation . . . . 67 5 5 Third Order Approximation . . . . 72 5.6 Summary . . . . . . . . . . 77 VI. SOLUTION OF SYSTEM EQUATIONS CONCERNING THE FLOW IN CONFINED AND UNCONFINED AQUIFERS . . . . . . . . . . . 79 6.1 Regional Groundwater Flow . . . 79 6.1.1 Computation of Piezometric. Heads . . . 80 6.1.2 Introducing the Dirichlet Boundary Condition . . . . 83 6.2 Solution of Flow Vectors . . . . 85 6.2.1 Direct Calculation . . 85 6.2.2 Simultaneous Calculation of Velocity Vectors . . . . . 86 Chapter VII. VIII. 6.3 Solution of System Equation for Flow in Unconfined Aquifer With Transient Phreatic Surface . . . 6.3.1 Background . . . 6.3.2 Locating the Phreatic Surface by Modifying the Elements . 6.3.3 Location of Phreatic Surface Using Fixed Elements . . 6.3.4 Steady- State Condition . . 6.3.5 Reasonable Time Step . 6.3.6 Iteration Within Time Step 6.3.7 Comparison Between the Two Methods of Locating the Phreatic Surface . . . . 6.3.8 Summary . . . . . . . SOLUTION OF SYSTEM EQUATIONS FOR CONVECTIVE-DISPERSION PHENOMENA . . . 7.1 Introduction . . . . . 7.2 Calculation of Dispersion Coefficients 7.3 Computation of Tracer Concentration 7.4 Discretization of Time Derivatives and Solution of System Equation . 7.5 Stability and Convergence Criteria NUMERICAL RESULTS FOR SIMULATION OF GROUNDWATER FLOW . . . . . . . . 8.1 Flow in a Confined Aquifer . . 8.1.1 Pumping in a Single Well Field . . . . . 8.1.2 Effects of Various Time Approximations on the Accu- racy of the Results . Comparison of Two Different Methods of Solution of the Velocity Vectors Flow in a Phreatic Aquifer . 8. 3.1 Transient Buildup of a Mound Due to Accretion . . . 8.3.2 Numerical Modeling of a Field Problem . . . . 3 Steady- State Solution . . .4 Two— Layer Aquifer . 5 Change of the Effective Porosity . 8.3.6 Effect of the Depth of the Water Table . . vi 91 100 114 114 115 117 120 121 121 122 123 126 127 129 129 129 134 139 143 143 148 153 156 161 161 Chapter 7 Page 8.3.7 Decay of a Groundwater Mound . 164 8.3 8 Maximum Applicable Time Interval . . . . . . . 164 8.4 Summary . . . . . . . . . . 170 IX. NUMERICAL RESULTS FOR PREDICTION OF CONCENTRATION OF A TRACER . . . . . . 171 9.1 Longitudinal Dispersion in Steady Uniform One-Dimensional Flow . . . 172 9.2 Sensitivity Analysis for Time Step and Grid Size . . . 177 9.3 Two— Dimensional Dispersion With Uniform Flow . . . 181 9.4 Point Source With Uniform Flow . . 186 9.5 Two- -Dimensional Dispersion With Transient Flow . . . . 188 9.6 Dispersion in a Phreatic Aquifer With Accretion . . . 196 9.7 Reasonable Time Step for Calculating Tracer Concentration in a Phreatic Aquifer . . . . . . . . . . 202 9.8 Summary . . . . . . . . . . 204 X. SUMMARY AND CONCLUSIONS . . . . . . 205 XI. RECOMMENDATIONS FOR FUTURE STUDIES . . . 208 LIST OF REFERENCES . . . . . . . . . . . 211 APPENDICES . . . . . . . . . . . . . 220 I. FINITE ELEMENT DEVELOPMENT . . . . . 221 II. INTEGRATED ELEMENT MATRICES FOR ONE- DIMENSIONAL QUADRATIC AND TWO- DIMENSIONAL TRIANGULAR ELEMENTS . . . . 236 III. DEVELOPED COMPUTER PROGRAMS . . . . . 245 LIST OF TABLES Boundary Conditions . . . . . . . Values of a, B and y for Second Order Approximation in Terms of e . . . Values of the Coefficients of Equation (5.1.5) for the First and Second Order Time Approximation . . . . Values of a, B and y for Third Order Approximation in Terms of 9 . . . Coefficients of Equation (5.1.5) for Third Order Time Approximation for Different 6 . . . . . . . . Different Combinations of Hydraulic Conductivities for Aquifer Shown in Figure 8-18 and Maximum Specified Time Step . . . . . . . . . viii Page 25 69 72 75 78 159 Figure 3-1 3-2 LIST OF FIGURES Definition of Piezometric Head in a Confined Aquifer . . . . . . . Vertical Cross Section of a Phreatic Aquifer . . . . . . . . . . Saturation Curve and Pressure Distribu- tion in the Capillary Zone . . . . Phreatic Boundary With Accretion . . Domain Divided Into Finite Elements . A Typical Finite Element . . . Interelement Zone Depicting How Continu— ous Function ¢ May Have Discontinuous Gradients as x1 + 0 . . . . . . One—Dimensional Finite Element Variation of d and 8 With 6 for First Order Approximation . . . . . . One—Dimensional Quadratic Element . . Variation of a and 8 With 9 for Second Order Approximation . . . . . . One—Dimensional Cubic Element . . Variation of a and 8 With 6 for Third Order Approximation . . . . . Dividing the Actual Discharge Curve Into a Series of Step Functions . Velocity Vectors in the x1, x2 Plane Movement of the Phreatic Surface Simple Linear Quadrilateral Isoparametric Elements . . . . . . . ix Page 20 20 27 30 39 52 54 63 65 67 70 73 76 82 87 93 95 Figure 6-5 Location of the Temporary and Actual Nodes . . . . . . . . . . . . Typical Grid System for Solving Phreatic Aquifer Problem . . . . . . . . Phreatic Surface Passing Through the Elements . . . . . . . . . . . One-Dimensional Linear Element . . . . One-Dimensional Quadratic Element . . . Location of the Phreatic Surface and Definition of the Terms . . . . . . Scheme for Solving the Piezometric Heads in Unconfined Aquifers With the Movable Node Technique . . . . . . . . . Scheme for Solving the Piezometric Heads in Unconfined Aquifers With the Fixed Node Technique . . . . . . . . . Schematic Representation for Calculating the Tracer Concentration in a Confined or Unconfined Aquifer . . . . . Nonleaky Artesian Aquifer With a Fully Penetrating Well . . . . . . Grid System Which Is Used in Simulating a Single Well Field . . . . . Comparison of the Numerical and Analyti- cal Solutions for Drawdown, 1697 and 2163 Meters Away From the Well . . . Finite Element Representation of a Con— fined Aquifer With Linear Isoparametric Quadrilateral Elements . . . . . Calculated Piezometric Heads Versus Time at Point W Using Different 6 With At = 0.05 Days . . . . . . . . Calculated Piezometric Heads Versus Time at Point W Using Different 6 With At = 1.0 Day Page 97 99 101 104 106 107 118 119 125 130 132 133 135 137 138 Figure 8-7 8-8 Computed Values of V at the Nodes by Direct and Simultaneous Methods . . Computed Values of V2 at the Nodes by Direct and Simultaneous Methods . . Simulation of Transient Buildup of a Mound Due to Accretion . . . . . Prediction of Transient Buildup of a Mound Due to Accretion With Dif- ferent Techniques . . . . . . . Portion of a Grid System . . . . . Prediction of Transient Buildup of a Mound Due to Accretion . . . . . Sketch of the Numerical Model for a Field Problem and Location of Water Table for Different Times . . Variation of the Relative Conductivity With Moisture Content for Two Soils . Representative Sketch of a Vertical Cross Section of a Phreatic Aquifer and Grid System . . . . . . . Location of the Phreatic Surface at the Steady-State Condition and Rise of Water Table With Respect to Time at the Center Line - . . - . . Equipotential Lines and Specific Dis— charge Vectors Beneath and in the Vicinity of the Recharge Zone at the Steady—State Condition . . . Vertical Cross Section of a Two-Layer Phreatic Aquifer . . . . . . Effects of Two-Layer Aquifer on the Water Table Profile After 5.1 Days . Equipotential Lines and Specific Dis— charge Vectors forzaTwo-Layer Aquifer Locathmiofthe PhreaticSurfaceAfter 5.1 DaysforlfifferentEffectivePorosities xi Page 140 141 144 146 146 147 '150 152 154 155 157 160 162 163 Figure Page I 8-22 Rise of Water Table for Different Ini- I tial Depths of the Free Surface ' After 10.0 Days . . . . . . . . 165 I 8-23 Rise and Fall of Water Table in a Phreatic Aquifer . . . . . . . . 166 I I 8-24 Rise and Fall of Water Table at the Center Line of a Recharge Zone . . . 167 8-25 Oscillation of Water Table Due to a Large Time Step . . . . . . . . 169 9-1 Cross Section of a Homogeneous and Iso- tropic Porous Medium With a Plane Source Maintained at x1 = 0 . . . . 172 9-2 Finite Element Model to Simulate One- Dimensional Longitudinal Dispersion . . 174 9-3 Numerical Solution of the Longitudinal Dispersion at t = 30 Sec. and t = 80 Sec. Using Second Order Time Approxi- mation and Comparison With the Analytical Results . . . . . . . 175 9-4 Numerical Solution of the Longitudinal Dispersion at t = 30 Sec. and t = 80 Sec. Using First Order Time Approxi— mation and Comparison With the Analytical Results . . . . . . . 176 9-5 Variation of ZT With Respect to At for Different e . . . . . . . . . . 179 9-6 Effects of Different e on the Accuracy of the Concentration Distribution at x = 2 Cm. With At = 7 Sec. . . . . . 180 9-7 Variation of ZX With Respect to Ax at t = 80 Sec. . . . . . . . . . . 182 9-8 Two—Dimensional Dispersion With One- Dimensional Flow . . . . . . . . 183 9-9 Steady—State Solution of Two—Dimensional Dispersion With One-Dimensional Flow . 185 9-10 Two—Dimensional Dispersion With a Point Source . . . . . . . . . . . 187 xii Figure 9-11 Page Representative Domain of a Porous Medium With Time Variable Velocity and Dis- persion . . . . . . . . . . 189 Variation of the Piezometric Head, Velocity, and Concentration With Time at Point A . . . . . . . . 190 ISO-Concentration Lines for a Medium With a Transient Flow Field After 100 Seconds . . . . . . . . . - . 191 Domain of Two-Dimensional Dispersion With a Transient Flow Field . . . . . . 193 Variation of Tracer Concentration With Time at Points A and B . . . . . . 194 Iso-Concentration Lines After 10 and 20 Days . . . . . . . . . . . 195 Representative Model Used in Simulating Tracer Movement in a Phreatic Aquifer . 197 Discharge Vectors in a Phreatic Aquifer With Recharge After 2 and 8 Days . . . 199 Iso—Concentration Lines in a Phreatic Aquifer After 4, 8 and 12 Days . . . 200 Variation of Concentration at Point B With Time Using Different Time Steps . 203 xiii aII m D: O O O O m D" ai-O * d)ij NOMENCLATURE longitudinal dispersivity transversal dispersivity coefficients in Equation (5.1.5) thickness of aquifer dimensionless length mass concentration unknown variable (in Chapter V) initial value of C analytical solution of C shifting distances of phreatic surface in the x1 and x3 directions shifting distance along the normal line shifting distance along the normal to the free surface average shifted distance domain of interest maximum specified time step molecular diffusivity mechanical dispersion coefficient hydrodynamic dispersion coefficient longitudinal dispersion coefficient transversal dispersion coefficient coefficient of molecular diffusion xiv M/L M/L Lz/T Lz/T L2/T L2/T L2/T L2/T FNT W W W) W U- 2° NELS superscript representing an element Fixed Node Technique acceleration due to gravity L/T2 subscript superscript representing iteration unit normal along xl accretion L/T magnitude of accretion along the L/T normal to the phreatic surface subscript subscript unit normal along x3 intrinsic permeability L2 hydraulic conductivity L/T directional cosine length L length of a segment L an integer number of nodes in one element order of approximation upper bandwidth of matrix plus one Movable Node Technique subscript effective porosity shape function number of elements in whole grid system XV NNDS NNFS NT NZ number of nodes in whole grid system number of nodes on the phreatic surface number of available data number of available data shape function for a specific point pressure discharge per unit area discharge per unit area at a specific point specific discharge seepage flux along the unit normal mass flux of source or sink known flux along boundary subscript distance from pumped well to obser- vation point drawdown storage coefficient surface boundary elastic specific storage specific yield time transmissivity medium's tortuosity velocity of propagation of the free surface propagation of the phreatic surface surface along the unit normal xvi M/LT L3/L2T L3/L2T L/T L/T M/L3T L/T L/T L/T ij At AX magnitude of velocity L/T seepage velocity L/T simultaneous solution of the L/T seepage velocity vectors weighting function Cartesian coordinates L location of the phreatic surface L at time t + At temporary location of the phreatic L surface at time t + At elevation head L averaged sum oftflmasquare of residu- als using different time steps averaged sum of the square of residuals using different Ax known functions for boundary condi- tions of dispersion coefficients which are used to approxi- mate unknown variable (Chapter V) an angle equal to (N/2) — w + 9 known functions for boundary condition of flow coefficients which are used to approxi- mate first derivative of unknown variable (Chapter V) coefficients which are used to approxi- mate second derivative of unknown variable (Chapter V) Dirac delta function Kronecker delta time step T grid size L xvii dimensionless curvilinear and local coordinates in x2— or x3—direction dimensionless time parameter angle that the tangent to the phre- atic surface makes with the positive xl-direction independent variable kinematic viscosity of fluid dimensionless curvilinear and local coordinates along xl-direction density of the fluid piezometric head average piezometric head known piezometric head angle that the tangent to the phre— atic surface makes with positive xl—direction a superscript representing the numerical approximation of unknown variable such as ¢ or C xviii L /T M/L CHAPTER I INTRODUCTION 1.1 Motivation In recent years the demand for fresh water has increased drastically due to population increase, indus- trial growth, and agricultural expansion. However, fresh water sources are limited and the bulk of this vital resource lies underground in the form of ground- water. It is now more or less generally accepted that man's future will depend on his ability to conserve this valuable resource for human consumption. Almost every pollutant that leaves a sewage plant or is part of runoff from rural lands or urban areas will reach some water basin such as a lake or an aquifer. Besides the customary disposal of treated effluent in nearby streams, rivers and lakes, irrigation of farmlands with secondary and tertiary wastewater is being sought as one of the possible ways to replenish groundwater sources while maintaining their high quality. In one type of system, treated sewage effluent is passed through a series of ponds before transferal to farmlands for use as irrigation water [see, e.g., Bahr 1974]. The tertiary wastewater slowly infiltrates into the soil while pro- viding intake water and nutrients for the plants. Dur- ing this passage the constituents are mixed, dispersed, and diffused through the flowing mass, and some of them such as phosphate, sulfate, and certain heavy metals are adsorbed by the soil [Ellis 1973]. After the water leaves the biologically active zone of the soil it enters the groundwater reservoir and eventually appears in wells and springs in the region. Although studies of nutrient intake by plants and studies of pollutants deposited along the course of the water's flow deserve great attention, it is important to know the concentration of the dissolved chemicals which travel easily with the water through the porous medium at different stages of the movement. The analysis of the effects of possibly contaminated recharged effluent on the quality and quantity of the transient groundwater is extremely worthwhile and beneficial. 1.2 Objective and Plan The objective of this thesis is to investigate the effects of treated wastewater recharge on the quality and quantity of the groundwater resources. Consideration is given to regional problems, but primary emphasis is focused upon the local region, that is, the region beneath and in the vicinity of the recharge site. This will be done by considering a two—dimensional horizontal plane for the regional scale (on the order of one square kilometer) and a vertical cross section for the local scale (one hundred square meters). For the purpose of this study the pollutants of concern are dissolved chemical substances such as chloride and chromium, which remain unaltered during the transport process. The dis— persion and convection of a tracer through a confined or unconfined aquifer with transient conditions are simula— ted by a Galerkin formulated finite element method. Calculated transient velocity vectors are used to obtain the dispersion coefficient, and they are introduced into the convective—dispersion equation. The specific objec- tives include: 1. Calculating the location of a transient phreatic surface in unconfined aquifers. 2. Investigating the tracer movement in a confined and an unconfined aquifer with a time-dependent flow domain. 1.3 Scope of Study The scope of the study presented in this thesis includes a detailed description of the mathematical equations for both flow and convective-dispersion in confined and unconfined aquif.rs, solutions of the above equations by the finite element method, comparison of the numerical results of this investiation with existing data, and finally the application of the employed tech- niques for more complex problems. The mathematical equation describing the magni- tude of the piezometric head in an unconfined aquifer is nonlinear, because of the existing phreatic boundary. The difficulties of solving related equations will increase when the boundary is transient due to accretion or other events. Furthermore, the mass transport equation is also a nonlinear equation because of the dependence of the convective terms and dispersion coefficients on the velocity COmponents. Both the dispersion and flow equa— tions have to be solved simultaneously or consecutively for each time step in order to predict mass distribution in a porous medium. It is customary to modify the numerical grid systems such that the movement of the phreatic surface can be handled. However, because of the changing location of the nodes, it is difficult to obtain the values of velocity components within the grid system so that they might be used in computing the convective terms and dispersion coefficients. In this study this difficulty has been overcome by adapting a technique pro- posed by France [1971] which can locate the phreatic sur- face within the grid system without repositioning the nodal coordinates. The method provides a feasible way to solve the convective-dispersion equation for uncon- fined phreatic aquifers. The velocity vectors play a dominant role in the accuracy of the predictions governed by the dispersion model. A simplified procedure for solving velocity com- ponents is introduced; this procedure provides continuous velocity values at the nodes. The technique also enables one to predict the tracer movement with a transient flow condition. The finite element formulation of the field prob- lem leads to a set of ordinary differential equations of the form [A]{C(t)} + [H] {g—E} = {F(t)} where C(t) is an unknown variable such as concentration or piezometric head. Discretization of the time deriva— tive of this equation is one of the major concerns for many investigators. With the finite element concept, recurrence formulae for the above equation for three different orders of time approximation are derived. In this process a simple procedure to obtain the finite dif- ference relation for a variable and its first and second derivatives is also shown. In this study the validity of the proposed tech— niques is established by first comparing the numerical flow results with existing analytical, experimental, and field data. Upon verification of the solution of the flow equation, the prediction of the movement of a tracer in an unconfined aquifer with a transient phreatic boundary is conducted. Numerical examples are presented to demonstrate the capability of the proposed techniques. Sensitivity analyses are made to explore the effects of time steps and element size on the accuracy of the numerical results. CHAPTER II LITERATURE REVIEW 2.1 Introduction The major objective of this study is to investi- gate the movement of a tracer in a confined or an uncon— fined aquifer which is experiencing a transient flow regime. The piezometric head and the velocity of ground— water flow must be known in order to predict the rate and direction of movement of dispersive substances. The simulation of mass transport through porous media in two- dimensional horizontal flow in a confined aquifer has received considerable attention, but little attention has been given to the solution of the dispersion equa- tion for an unconfined aquifer with a phreatic boundary. The difficulties associated with calculation of the convective term and dispersion coefficient with a time- dependent phreatic boundary might be considered the major obstacles in the way of progress in this area. There is a vast number of publications available concern- ing the solutions of the flow equation and the convective- dispersion equation; only those publications relevant to this study are referenced herein. 2.2 Mechanics of Flow 2.2.1 Confined Aquifer The finite element solution of the differential equation describing flow in a two-dimensional horizontal plane for a nonhomogeneous and anisotropic confined aqui- fer is well-established and is available in the literature. Zienkiewicz et a1. [1966] employed the method to obtain a steady-state solution for the heterogeneous and aniso— tropic seepage problem. Javandel and Witherspoon [1968] used the Rayleigh—Ritz procedure to solve groundwater problems with linear triangular elements. Pinder and Frind [1972] utilized Galerkin's technique to analyze groundwater problems with isoparametric elements. The numerical model developed in this study to calculate the piezometric head in a confined aquifer is similar to the procedure presented by Pinder and Frind. 2.2.2 Unconfined Aquifer With Phreatic Surface The equations governing boundary and initial value problems of the phreatic surface for liquid flow through porous media are known [Polubarinova-Kochina 1962, Bear 1972]. The exact analytical solutions of these equations are extremely complex. To simplify the treatment of such problems, Dupuit in 1863 assumed that the gradient of the phreatic surface in the vertical plane away from the wells and mounds is very small, thus the groundwater flow is essentially horizontal and can be considered as a uniform flow. This assumption led to the well-known Boussinesq equation. Because of the nonlinearity of the Boussinesq equation, only a small number of analytical solutions are known to date [e.g., Polubarinova-Kochina 1962, and Bear 1972]. An approach to overcome the problem is to linearize either the par— tial differential equation describing the phreatic sur- face boundary or the Boussinesq equation. Bear [1972] has outlined several linearization techniques with rela— ted references. Marino [1967], following Hantush [1963, 1967], employed a linearization method to solve the problem of the rise and decay of a groundwater mound below a spreading site, and justified the solution with experimental study. All analytical solutions are limited to flow systems in which the boundary conditions are simple, the porous medium is relatively uniform, and the Dupuit approximation is valid, i.e., the vertical gradi- ents throughout are not too large and hence are negligible. Models and analogues are tools for achieving the solutions of problems where the direct analytical solu— tion is not possible because of the complexity of the system. The Hele-Shaw or viscous flow analog and resist- ance network analogues are most commonly used. Bear [1960] discussed the scales of viscous analog models for groundwater studies. Bear [1972] presents extensive 10 bibliographies of past applications of the Hele-Shaw model to studies of groundwater flow. Marine [1967] used the viscous analog model to study the growth and decay of groundwater ridges. Tinsley and Ragan [1968] employed the Hele-Shaw model to investigate the response of an unconfined aquifer to localized recharge. Herbert [1968] used a resistance network analog to study the time-variable movement of the water table in unconfined saturated strata. His work is based on the assumption that for any time-variant system of water, the flow can be approximated as a series of steady-state solutions, each of slightly varying shape and satisfying the Laplace equation. In his technique the location of the water table is known at a given time, and the new position of 'the free surface is predicted for a chosen finite time interval. This process is repeated for several time intervals until the maximum required time has been reached or the system has reached the steady state. In the solution of groundwater flow problems, the digital computer offers enormous advantages and has become a dominant computational tool. There are several numeri- cal techniques available for solving the governing equa- tion for groundwater flow. Of these, the finite difference and the finite element techniques are most commonly used. Todsen [1971] used the finite difference method for solving the free surface flow problems. Amar 11 [1975] investigated the two—dimensional hydrodynamic behavior of recharge of an unconfined aquifer with the finite difference technique. The finite difference technique is simple to program, but manipulation of curved boundary conditions which most likely appear in nature is difficult. The finite element technique eliminates this problem and the computer model can be developed in such a way that it can be used for any type of boundary without modifi— cation of the program. Taylor and Brown [1967] presented finite element solutions of steady seepage through dams using a network of triangular elements. In their technique, the location of the phreatic surface was guessed and subsequently adjusted until the free surface boundary conditions were satisfied. Neuman and Witherspoon [1970, 1971] improved and extended their technique to problems of steady—state and transient seepage with a free surface using linear triangular elements. Desai [1972] used the finite ele— ment procedure with isoparametric elements to analyze transient unconfined seepage under drawdown conditions in porous media. Sandhu et a1. [1974] introduced the variable time step analysis of unconfined seepage. France et a1. [1971], following Herbert [1968], used the finite element method with isoparametric elements to analyze 12 the free surface seepage problem. France [1974] has extended this work to three-dimensional problems. In all numerical solutions it is necessary to modify the elements to accommodate the movement of the phreatic surface. Based on the characteristics of curved isoparametric elements, France [1971] introduced a new method which permits one to locate the phreatic surface without altering the position of the nodes for each element. This procedure will be called "location of the phreatic surface by use of fixed nodes." The technique is adopted and modified in this study to determine the location of the free surface with transient recharge. This method will provide a tool to solve the convective—dispersion equation in unconfined aquifers with a time-dependent phreatic boundary. 2.3 Mechanics of Convective- DispersIon Phenomena The analytical solution of the convective— dispersion equation, except for a small number of simple one- and two-dimensional cases, is not easy to determine. Ogata and Banks [1961] used Laplace transforms to obtain the solution of the one—dimensional longitudinal disper- sion equation. Harleman and Rumer [1963] gave a steady— state solution for two—dimensional dispersion. Bruch and Street [1967] formulated the analytical solution for unsteady dispersion in an idealized study of 13 one-dimensional seepage flow through an isotropic porous medium. Hoopes and Harleman [1967] introduced an expres- sion for the distribution of dissolved concentration substances which were added to the steady-state flow between a recharging and pumping well in a homogeneous, isotropic aquifer of infinite, horizontal extent. Marino [l974d] gave a mathematical solution to predict the distribution of concentration in saturated porous media I resulting from a variable source concentration. Much experimental work has been attempted to investigate the behavior of dispersion coefficients and ‘ their relation to the seepage velocity, porous structure, and concentration gradient; literature concerning this subject is given by Bear [1972]. Numerous investigators such as Bear [1961], and de Josselin de Jong and Bossen [1961], showed that the dispersion coefficient is a function of true velocity and medium properties. Rumer [1962] experimentally determined the longitudinal dis- persion coefficient for one—dimensional transient flow within a certain range of values of the Reynolds number. Harleman and Rumer [1963] investigated the dependency of the dispersion coefficient upon the Reynolds number and porous structure. They conducted laboratory experiments to study the convection and dispersion of salt water in a two—dimensional confined aquifer. Shamir and Harleman [1967] developed an analytical solution for two problems 14 of dispersion in layered porous media and verified their results experimentally. Fattah [1974] investigated and verified a model of the dispersion coefficient tensor in flow through anisotropic and homogeneous porous media. The conclusion that can be made from the results of these investigations is that the dispersion coeffi- cient is a second-rank tensor and is a function of the true velocity vectors, porous media properties, and the Peclet number. However, there is still no universal agreement regarding the degree of dependency of the dis- persion coefficient on these parameters. In the simulation of the movement of a tracer in a porous medium, the flow and convective-dispersion equations are solved simultaneously or consecutively. The following paragraphs concern only those techniques which are used to solve the partial differential equa- tion describing the dispersion of a dissolved—chemical constituent in a medium. The finite difference method is the most com— monly used scheme in attempting the numerical solutions of the mass transport equation. Douglas et a1. [1959] employed an alternating direction—implicit procedure to solve a two-dimensional, two-phase, incompressible flow model. A similar technique was used by Peaceman and Rachford [1962] in calculating the multidimensional miscible displacement. The Crank—Nicolson approximation 15 is frequently used in the area of mass displacement. Fried and Combarnous [1971] summarize some of the ana- lytical and numerical methods of resolution of the convective-dispersion equation. Pinder and Cooper [1970], and Reddell and Sunada [1970], applied the characteristics approach to solve salt water intrusion including the effect of dis— persion. Bredehoeft and Pinder [1973] used similar concepts to investigate the groundwater contamination at Brunswick, Georgia. Konikow and Bredehoeft [1974] used the method of characteristics to investigate the chemical quality changes in an irrigated stream-aquifer system. Robertson [1974] used the same approach to model the transport of radioactive and chemical waste in the Snake River Plain Aquifer. The method of character- istics involves placing several moving particles in each cell of the finite difference grid. The location and concentration associated with each particle varies with time. Although this method gives good results compared to the analytical solution and is simple in concept, it is tedious to program and is suitable only for Specific situations commonly encountered in the field [Pinder 1973]. In an effort to circumvent difficulties associ- ated with the method of characteristics, Price et a1. [1968] introduced a Galerkin-based variational method to 16 approximate the solution of the dispersion equation, in which various different base functions were used. The finite element technique has been_used recently as a numerical tool to solve the convective— dispersion equation. Guymon et a1. [1970] and Nalluswami et a1. [1972] employed a finite element integration scheme with triangular elements. Cheng [1973] solved the convective—dispersion equation based on the Galerkin procedure, using a family of triangular elements or ; quadrilateral isoparametric elements. Pinder [1973] also used the Galerkin-finite element formulation to simulate the groundwater contamination in Long Island, New York. Wang and Cheng [1975], using Dupuit's assump— tion, solved the convective-dispersion equation by quadratic isoparametric elements for homogeneous and isotropic media with uniform horizontal flow and constant dispersion coefficients. Segol et a1. [1975] realized that the velocity vectors used in the mass transport equation should be continuous across elements. Thus, they solved three equations (two components of Darcy and one of mass conservation) simultaneously, and investigated the distribution of salt concentration with steady flow. In the present work the technique for solving the convective-dispersion equation in unconfined aquifers with a transient phreatic surface is presented, and higher order approximations of time—dependent variables are introduced. CHAPTER III MATHEMATICAL REPRESENTATION OF PHYSICAL SYSTEM 3.1 Background In this study the dispersion of a tracer on a regional as well as a local scale is investigated. For the regional scale the two-dimensional horizontal plane is considered, while for the local scale a vertical cross section of a site is chosen. To predict the move- ment of a tracer in a porous medium, the flow regime and its behavior in the medium should be well understood. In this chapter the system is defined, the mathematical descriptions of flow in the aquifers with and without a phreatic boundary are given, the convective—dispersion equation is presented, and the initial and boundary con— ditions are discussed. 3.2 Darcy Equation Darcy established a linear relationship between the seepage velocity and the gradient of the piezometric head. This law, which is a consequence of the equation of motion neglecting inertia effects, can be generalized 17 '9. 18 for either a two— or a three—dimensional situation. It is given by: qi Ki. a¢ . . Vi = H; = ‘ 5;; 5;; 1,J=l.2.3 (3.2.1) where Vi is the seepage velocity, qi is the specific discharge, Kij is a second order tensor whose elements are called the hydraulic conductivities, ne is the effec- tive porosity of the aquifer, ¢ is the piezometric head, and xj are the Cartesian coordinates. The hydraulic conductivity is a scalar coefficient which depends on both solid matrix and fluid properties. It is defined as = .11.. - -= Kij V 1,3 1,2,3 (3.2.2) where kij is the intrinsic permeability of the porous matrix and depends solely on the properties of the solid matrix, g is the acceleration due to gravity, and v is the kinematic viscosity of the fluid. From purely physical considerations, it would seem that the hydraulic conductivity tensor must be sym- metric [Eagleson 1970], in which case K = K21, K13 = 12 and K = K and its components reduce to six. K31' 23 32' Since the principal axes of the symmetric permeability tensor will be orthogonal, it is possible to orient the 19 coordinate axes (x1, x2, x3) parallel to the principal axes so that only the three orthogonal terms remain. Thus _ Kll 3¢ Vl _ - n 3x e 1 K22 3¢ V2 = ' n— Dx— e 2 K 8¢ V3 = - Iii QT (3.2.3) e 3 For incompressible fluid, the piezometric head is defined (see Figures 3—1 and 3-2) ¢ =—P—+ x (3.2.4) where p is the pressure deviation from atmospheric pres— sure, p is the fluid density, and X3 is the elevation above datum. In Equation (3.2.4) the term p/pg is called the pressure head, and X3 is known as the elevation head. For confined aquifers the piezometric surface is an imaginary surface to which water rises in a tapped well (Figure 3-1). In unconfined aquifers the piezo— metric surface coincides with the upper surface of the zone of saturation, called the water table or phreatic surface,where the pressure is atmospheric (Figure 3—2). 20 .__ Ground 775 7F ’5 7‘: :1: IE I“ surface 2? Y Piezometric T surface P/Y ’ ’ , I . Wlmpermeable X3 (I) b X2 ' Impermeable I Datum Figure 3 1 —-Definition of Piezometric Head in a Conf1ned Aquifer. accretion) I I I I I I arm—WWW“ surface 2: Phreatic ' surface X3 X2 X1 L I4 Impervious boundary and datum level —J I‘ I] Figure 3-2.--Vertica1 Cross Section of a Phreatic Aquifer. 21 3.3 Regional Groundwater Flow For regional problems, two-dimensional horizontal flow is considered. The governing equations are well established [e.g., see Bear 1972, Pinder and Frind 1972]. 3.3.1 Basic Assumptions The following assumptions are valid for regional groundwater flow: (a) The flow is essentially horizontal in a two— _dimensional plane. This assumption is valid when the variation of thickness of the aquifer is much smaller than the thickness itself. This approximation fails in regions where the flow has a vertical component. (b) The fluid is homogeneous and slightly compressible. (c) The aquifer is elastic and generally non— homogeneous and anisotropic. The consolidating medium deforms during flow due to changes in effective stress with only vertical compressibility being considered. (d) For the two—dimensional horizontal flow assumption, an average piezometric head is used where the average is taken along a vertical line extending from the bottom to the top of the aquifer, i.e., b _ l where b is the thickness of the aquifer. 3.3.2 Two-Dimensional Horizontal Flow The combined equation of motion and continuity for flow in a two-dimensional horizontal plane may be written 3 23¢ _ _ 23¢ . ._ where Tij is the transmissivity tensor equal to the aquifer thickness multiplied by the hydraulic conduc- ., S is the storage coefficient, t is time, 3 I is the vertical recharge or infiltration into the t1v1ty Ki aquifer, and P is strength of a sink (or source) function [Pinder and Frind 1972] defined by M P = mgl Pw[(xl)m,(x2)m]6[xl-(xl)m][x2-(x2)m] where Pw is the discharge (or recharge) from the aquifer, M is the number of nodes in one element (details are given in Chapter IV), and 6 is the Dirac delta function. 3.3.3 Initial and Boundary Conditions 3.3.3.1 Boundary conditions.—-In order to solve a partial differential equation describing a physical phenomenon, it is necessary to choose certain additional conditions imposed by the physical situation at the 23 boundaries (S) for the domain (D) under consideration. In general the equation for the boundary condition can be written BlT.. 39— A. + 32¢ + 33 = o i,j=l,2 (3.3.2) 13 xj 1 where 9i are the directional cosines, and 81, 82, and B3 are given functions of position and possibly time. For flow through an aquifer, three different boundary conditions are applicable: (a) Dirichlet or prescribed potential: In this case the potential is specified for all points along the boundary ¢=_I5—; 82#0 (b) Neumann or prescribed flux: Along a boundary of this type, the flux normal to the boundary surface is prescribed for all points of the boundary as a function of position and time A special case of the Neumann condition is the impervious boundary where the flux vanishes everywhere on the boundary, i.e., 83:0 24 (c) Cauchy boundary: This problem occurs when the potential and its normal derivative are prescribed on the boundary in the combined_form, and the entire Equa- tion (3.3.2) is used. Different forms of Equation (3.3.2) for three types of boundary conditions are summarized in Table 3-1. In general, for a flow problem one will have mixed boundary conditions in which the Dirichlet condi— tion will apply over a part of the boundary and the Neumann condition will be specified for the remaining portion [Bear 1972]. 3.3.3.2 Initial conditions.-—At the initial time, either the piezometric heads are known in the entire domain (D) or the hydrologic stresses (such as pumping and recharge) are specified and boundary condi- tions are known. For the second case the system has reached the steady state, so the solution of the equation 3 8¢ _ _ . ._ W [Ti]. 3X] P + I — O l,j—l,2 (3.3.3) 1 3 will yield piezometric heads for the initial time. The procedure of solving Equation (3.3.3) is discussed in Chapter VI. H H x ma Tm - I 9 Im + .. II@ ..a o A o R aaosmo m mm mm em .ommu ma mummcson map am a mxm ma so 30am may .0 n ma DH mm I N .a em ..a o o A acmssoz m .A we AD .wi N.MI I N e o N o DoaQOHHHQ a 5 m 2 Mnmfimm QOHumsqm mm Hm oEmz meme .OEHD can mommm mo mcoflpoqsm :3ocx who mm paw mm.am “mocflmoo Hmcoflwomuflp mum Ha muons m . ma 0 n mm + 9mm + as Iwm ..efim copufluz on cmo COHuHccoo mHMUQDOQ map How coflpmoqo Hmuocom one .mcofluflccoo mnmccsomuu.aum mqmae 26 3.4 Phreatic Aquifer With Accretion A typical cross section of a phreatic aquifer with accretion is illustrated in Figure 3—2. The govern- ing equations are discussed in literature [e.g., Bear 1972, France 1974]. 3.4.1 Basic Assumptions 1. Usually immediately above the water table (p = 0) there is a zone that is saturated or nearly so. This nearly saturated zone above the phreatic surface is called the capillary fringe or capillary rise, where the pressure is negative. In Figure 3—3 a typical saturation curve and pressure distribution in the capillary fringe at equilibrium areshown. The capillary rise might range from 2-5 cm for coarse sand up to greater than 200 cm for clay [Bearl972]. In this study it is assumed the aquifer is fully saturated and the capillary fringe can be ignored. The resulting idealized diagram is given in Figure 3-3b. 2. When the saturated soil is being drained, the free surface gradually descends and some water is removed from the soil profile. In practice, the amount of water removed per unit volume of soil depends upon the water level, rate of drawdown, temperature, and atmospheric pressure; but for theoretical analysis it is usually taken as a constant and equal to the specific 27 Saturation curve Assumed top of capillary fringe —~-—._-—--- — —¢—-—I Water table p<0 C _‘ o,./,Ix3=0.p=o) p - Sw (saturation) l Capillary fringe Pressure distribution (a)[After Bear 1972] Idealized ’/// saturation curve Idealized pressure distribution (b) Figure 3-3.--Saturation Curve and Pressure Distribution in the Capillary Zone. (a) Actual form, (b) idealized form. 28 yield. The specific yield Sy is thus defined to be the volume of water drained over the gross volume of the porous medium. Quantitative information on specific yield is given by Todd [1959, pp. 23-26]. Because of the assumptions made in the above paragraph, in this study the magnitude of the specific yield is assumed constant and equal to the effective porosity which is defined to be the volume of water drained by gravity Wfrom a unit volume of saturated soil. 3. In unconfined aquifers the amount of water released from storage is usually small compared to the water available from normal movement through the aquifer and accretion, thus elastic specific storage (SS) can be ignored [Herbert 1968; Neuman and Witherspoon 1971]. 3.4.2 Governing Differ— ential Equations Consider an unconfined aquifer with its phreatic surface depicted schematically in Figure 3-2. The governing differential equation can be written 3 23¢ _ 23¢ . ._ — [KN —_] — sS E 1,3—1,3 (3.4.1) where SS is the elastic specific storage, and the other terms are defined previously. As discussed above, if Ss can be neglected then Equation (3.4.1) will have the form: 3 29¢ _ . ._ 'Wi [Kij E] — 0 l,j—l,3 (3.4.2) 3.4.3 Boundary and Initial Conditions Let D represent the flow region, which in general may possess up to four kinds of boundary conditions: in addition to the Dirichlet and Neumann conditions, a phreatic (free surface) boundary and a seepage face. Referring to Figure 3—2, the following boundary condi— tions can be written: 1. At x = 0 and x = L, the piezometric heads 1 l are known functions of time (Dirichlet) ¢(OI X3; t) = (I) ¢(L, x3, t) = CI) (3.4.3) 2. On the impervious boundary B, the normal flow is zero 8¢(Xl,0,t) . . Kij ————§§f——— ti = 0 1,j=3 nglgL (3.4.4) ‘ J 3. In the concept of successive changes of steady—state values, it is assumed the flow at each instant is steady but its boundary condition is time variable [Polubarinova—Kochina 1962]. Therefore, the flow rate for a small time increment is equal to the 30 change of volume filled with fluid divided by the time interval. Consider Figure 3-4 where the position of the phreatic surface at times t and t + At is shown. In represents the rate of accretion normal to the phreatic surface. By taking a control volume in the direction of the unit normal between two successive positions of the boundary at times t and t + At, and writing the continuity relation, one arrives at I (3.4.5) Un is the propagation of the phreatic surface, qn is the seepage flux, both along the unit normal, and ne is the —zs'—ps“’1RF"‘Z“"7§7"'757_~ Ground surface X1 Figure 3-4.--Phreatic Boundary With Accretion [Adapted from Todsen 1971]. 31 effective porosity. Equation (3.4.5) is a nonlinear boundary condition for Equation (3.4.2) because it con- tains an unknown dependent variable, i.e., ¢, in the flux term. As will be discussed in detail in Chapter VI, by assuming that at the beginning of each time increment the piezometric heads are known at the phreatic surface, Equation (3.4.5) is linearized. Equation (3.4.5) after multiplying by At can be written as U.9.. At = U At = A—t— (-K.. E— 9.. - I-IL) (3.4.6) 3 n 1 1 3 1 n . e 3 i,j=l,3 where I is the accretion (positive downward), and Uj is the velocity of propagation of the free surface at the point of consideration on which the pressure is maintained atmOSpheric. Todsen [1971] has also derived Equation (3.4.5). 4. On the seepage face, Initially, the surface configuration and the boundary conditions are known. 3.5 Convective-Dispersion Phenomena In this study the movement of a solute in a saturated flow through a porous medium is considered. 32 This solute will be referred to as a "tracer." The symbol C will be used to denote the concentration of a tracer, i.e., mass of tracer per unit volume of solution. The term tracer will be used to represent any species of interest in a solution. 3.5.1 Basic Assumptions 1. It is assumed that no chemical reactions occur between the water and the aquifer or soil material that affect the tracer concentration. 2. The porous medium is homogeneous and iso- tropic with respect to dispersivity. 3. The flow regime is laminar. 4. In general, variations in tracer concentra- tion cause changes in the density and viscosity of the liquid. These in turn affect the flow regime (i.e., velocity distribution). At relatively low concentrations it is assumed that the concentration does not affect the liquid properties [Bear 1972]. This assumption leads to the following conclusions: a. the viscosity is constant, b. the concentration does not affect the velocity distribution. 3.5.2 Hydrodynamic Disper— siOn Coefficient in a Porous Medium Hydrodynamic dispersion is the macroscopic out— come of the actual movements of individual tracer parti- cles through the pores and includes two processes [Bear 1972]. One mechanism is mechanical dispersion, which depends on both the flow of the fluid and the character— istics of the porous medium through which the flow takes place. The second process is molecular diffusion which basically results from variations in tracer concentration within the liquid phase, and is more significant at low velocities (e.g., less than 1 cm/hr). Thus the coeffi- cient of hydrodynamic dispersion Dij includes the effect of both the mechanical (or convective) dispersion Dij and * molecular diffusion (Dd)ij' Hence * ' = Dij Dij + (Dd)ij (3.5.1) - * * n In Equat1on (3.5.1), (Dd)ij — Tide’ where Dd 15 the molecular diffusivity and Tij is the medium's tortuosity. . is 3 approximately equal to 2/3 [Bear 1972, pp. 109-112]. * For homogeneous and isotropic media the value of Ti For most situations the contribution of molecular diffu- sion to hydrodynamic dispersion is negligible when compared to the mechanical dispersion coefficient. For example, for a gravel with seepage velocity ranging from 34 0.1 to 0.45 cm/sec, the magnitude of the dispersion coefficient varies from 0.01 to 0.08 cmz/sec [Rumer 1962]. The molecular diffusivity for solutes in water is very small and in the range of 0.5 to 4.0 X 10"5 cmz/sec [Welty 1969, p. 461]. Many investigators have attempted to model the dependence of the hydrodynamic dispersion coefficient on media, fluid properties, and flow characteristics, in order to understand the dispersion process in flow through porous media. A comprehensive discussion of the factors affecting the dispersion coefficient can be found in Bear [1972, pp. 605-616]. The mechanical dispersion coefficient for an iso— tropic medium in Cartesian coordinates can be written [Bear 1972] as: <3 .v+ (a -aII) vi Vj/V 1,3=1,2,3 (3.5.2) Dij =aII 13 I In Equation (3.5.2) aI and aII are the longitudinal and transversal dispersivities of the medium, respectively, Vi and Vj are components of the seepage velocity in the i and j directions, V is the magnitude of the velocity, and Gij is the Kronecker delta. Its value is one when i = j and is zero, otherwise. Equation (3.5.2) is com— monly used by investigators to calculate the mechanical 35 dispersion coefficient and hence is utilized in this study. It includes the major parameters causing the mechanical dispersion, and for any practical study it is assumed adequate. 3.5.3 Convective-Dispersion Equation in Cartesian Coordinates The equation describing the mass transport and dispersion of dissolved chemical constituents in a saturated porous medium may be written as 3c 2) 8 . 3C ' - [D—t + 5;.— “’1' 1. [fix— [Di]. 5?: + qc]— 0 (3.5.3) 1 l j (l) (2) (3) (4) i,j=l,2,3 where C is the mass concentration of the tracer; Dij is the coefficient of hydrodynamic dispersion, discussed in Section 3.5.2; Vi is the component of seepage velocity; go is the mass flux of source or sink; and xi is the Cartesian coordinate. The theoretical basis and the derivation of the diffusion—convection equation are dis— cussed in detail by Reddell and Sunada [1970], Bear [1972], and Bredehoeft and Pinder [1973]. In Equation (3.5.3) the first term represents the time rate of change of the tracer concentration. The second term describes the convective transport of C in the xi- direction, which is proportional to the seepage velocity. 36 The third term is the transport (redistribution) of C due to dispersion and molecular diffusion. Finally, the last term represents the time rate of production or decay of C. The convective-dispersion equation is a nonlinear partial differential equation of parabolic type. The relation is nonlinear because of the convective term, and because of the transport coefficient which is a function of the dependent variable V. The convective term “3/3XB(ViC)) is nonsymmetric and has been a principal source of difficulty in the numerical solution of the convective-dispersion equation [Guymon et a1. 1970]. 3.5.4 Initial and Boundary Conditions 3.5.4.1 Boundary conditions.--The general equa— tion of the boundary conditions for the mass transfer equation is similar to the flow equation. As discussed in Section 3.3.3, it can be written: I BC + _ alDij 5;; 2i d2C + d3 — 0 (3.5.4) where al, a2, and a3 are known functions. Three differ- ent boundary conditions are: (a) Dirichlet or prescribed concentration boundary condition: 37 C=-——,0L=0;0I27‘0 (b) Neumann or prescribed flux: '99.. =__§_ = Dij ax. £1 a ' 0‘1 I 0 ' 0‘2 0 j l for G3 = 0, one has the no-flow boundary. (c) Cauchy boundary: , ml and d2 # 0 Again, as in the flow situation, usually along the boundary one has mixed boundary conditions, i.e., the Dirichlet condition applies over a part of the boundary and the Neumann condition applies over the remaining part. 3.5.4.2 Initial conditions.--As an initial con- dition, the concentration distribution at some initial time t = 0 at all points of the flow domain must be specified: C(xi,0) = fl(xi) where fl is a known function of Xi' 38 3.6 Closure In this chapter the mathematical description of flow and dispersion phenomena in porous media is pre- sented, and basic assumptions are introduced. Using the concept of successive changes of steady-state values, the nonlinear boundary condition for an unconfined aquifer with a phreatic surface is linearized. For practical purposes the hydrodynamic dispersion coefficient can be replaced by the mechanical dispersion coefficient in predicting the tracer movement. Finally, the flow and convective—dispersion equations can be solved consecu— tively. In the following chapter, the finite element formulation of the flow and mass transfer equations is given. CHAPTER IV FINITE ELEMENT FORMULATION The finite element method is a numerical tech— nique which is used to approximate a continuous partial differential equation in a given domain D with specified boundary conditions along boundaries S. The key features of the finite element concept are [Norrie and de Vries 1973]: 1. The domain is divided into subdomains or finite elements, usually of the same form. 2. The trial solution is prescribed (func- tionally) over the domain in a piecewise fashion, element by element. 4 4 Domain De S \4 / (element) Nodes Domain D Figure 4-l.-—Domain Divided Into Finite Elements. 39 40 A detailed formulation of the finite element method is given by Zienkiewicz [1971], Norrie and de Vries [1973]. This technique has been utilized by several investiga— tors [Javandel and Witherspoon 1968; Pinder and Frind 1972; Neuman and Witherspoon 1971; Desai 1972; Cheng and Li 1973; and France 1971, 1974] to solve transient flow problems in a confined or unconfined aquifer. Recently the finite element procedure was also used to solve the convective-dispersion equation [Cheng 1973, Pinder 1973, and Segol et a1. 1975]. In this chapter a brief discussion of the Galerkin based finite element technique is given and the method is used to discretize the space derivatives of the flow and dispersion equations. The simultaneous solution of velocity vectors is also described, that is, the Galerkin formulation of the Darcy law is constructed and velocity components are calculated at the nodes. 4.1 The Galerkin Finite Element Method While the approximate minimization of a func— tional is the most widely accepted means of arriving at a finite element representation, it is by no means the only possible approach. The Galerkin method offers an alternative way to formulate a problem for the finite element solution without using variational principles. In the finite element technique the domain D is divided into subdomains De which are called elements. 41 Each element is designated by nodes. In this chapter NELS represents the number of elements, M is the number of nodes in each element, and NNDS stands for the total number of nodes in domain D (Figure 4—1). Consider a problem of solving approximately a set of differential equationsixlwhich the unknown function {C} has to be satisfied in the domain D with the boundary conditions specified along S. The governing equation can be written f({c}) = 0 A Let the trial solution for this equation be C A M c = [N]{c} = 2 N (2 (4.1.1) n=1 n n where [N] = [N(xi)] are shape functions (prescribed func- tions of coordinates) and {C} = {C(t)} is a set of M unknown parameters. In general, the equation of residual (or error) is formed in the following way: R = f e({Cl) - f each = -f eucn ;é 0 (4.1.2) D D D The best solution will be one in which the residual R has the least value at all points in the domain De. An obvious way to achieve this [Zienkiewicz 1971] is to make use of the fact that if R is identically zero else- where, then 42 f WRdD= 0 (4.1.3) De where W is any function of the coordinates. If the num- ber of unknown parameters {c} is NNDS and NELS linearly independent functions W are chosen, one can write a k suitable number of simultaneous equations as fDe wk R dD = fDe wk f([N]{C})dD = {0} (4.1.4) where Wk is called the weighting function. If the shape function Nk is to be chosen as the weighting function, the process is termed the Galerkin procedure. The ele- ment equations can be assembled by NELS [I e w R dD = 0 (4.1.5) D e=l to yield the global relations for domain D. 4.2 Finite Element Formulation of Flow Equation The residual equation for flow in a confined horizontal aquifer (Equation 3.3.1) can be written as A : i + — __a__ 3d.) . .2 R s p axi [Tij ij] 1,3 1,2 (4.2.1) 43 The symbol A represents the numerical approximation of ¢. Substituting Equation (4.2.1) into Equation (4.1.4), one obtains ¢ 3 _ Ll:5—+P’?T[uax[hk®-0 k4“.M (L22) By use of the Green theorem, the third term can be modi— fied + I N T.. 351 A. as (4.2.3) Se k 1] 8x 1 The last term in Equation (4.2.3) is nonzero only for elements which contain the Neumann flux boundary condition a$ _ 52, = Ise Nk Tij 3Xj 1 d8 [Se Nk Q2 dS (4.2.4) where O2 is known flux along the boundary. Substituting Equation (4.2.3) and Equation (4.1.1) into Equation (4.2.2) and rearranging the terms, one obtains 3Nk BNn IDe ¢n T13 Di? 3x. a” + I e S Nk Nn DE dD + I e Nk P d0 1 j D D + [s9 Nk QW d5 = 0 (4.2.5) 44 Since ¢n and its time derivatives are independent of the coordinates, they can be taken out of the integrals. Equation (4.2.5) can be written in matrix form [B]e{¢}e + [H]e §% {D}e = {F}e (4.2.6) where e e 8Nk BNn [B] =Bkn= e Ti.5—X— 6):— 6D (4.2.7a) D 3 i 3' k,n=1,...M e_ e _ [H] -—Hkn—!DeSNandD > (4.2.7b) i,j=1,2 {F}e=F}eW=—[eNkQ2dS-[ePdeD (4.2.7c) s D ) It is assumed that the storage coefficient is constant throughout the element and that the element coordinate axes coincide with the principal direction of the transmissivity tensor: the transmissivity can be defined either at the nodes or at each element. Evalu— ation of Equation (4.2.7) for different types of elements is discussed in Appendix II. Upon evaluation of Equation (4.2.7) for all elements and transformation to a global coordinate system, they are assembled by virtue of Equa- tion (4.1.5) into a global relationship IBIIAI + [Hugh = {F} (4.2.8) 45 The parameter {¢}, matrices [B] and [H], and force vector {F} are the summation of the corresponding terms in Equation (4.2.7) over all the elements in the domain D. [B] and [H] are banded symmetric matrices. Equation (4.2.8) is a set of first order linear differential equa— tions with unknowns {¢}amuican be solved simultaneouslyeu: the given nodes in the space domain. The solution of Equation (4.2.8) and a similar equation which is genera- ted from the finite element formulation of the convective- dispersion equation is presented in Chapter V. The governing equation for the unconfined aquifer with phreatic surface is 351] = o i,j=l,3 (3.4.2) 3 3 DR? [K13 3x. 1 The finite element formulation has the form [B']e{}e = {F}e (4.2.9) where 8N 8N i,j=l,3 I e _ k___Il “3 I - IeKij Wax. k,n=l,...M (4.2.10.4) D 1 3 e {F} __Ise Nk Q2 dS (4.2.10b) Assembling the element matrices leads to a system of equa- tions in the form [B]{¢} = {F} (4.2.11) 46 4.3 Finite Element Formulation of Convective-Dispersion Equation The residual of Equation (3.5.3) for each ele- ment can be written __3C 3 +_.__(cvi) -—8—.[D' 3C] ~83 i,j=1,2 (4.3.1) R-——— ..——— 3t axi Bxi 13 ij c Again, C represents the numerical approximation of C. In this development it is assumed that at every small time step the velocity vectors Vi and the dispersion coefficient Dij are known functions which are either determined independently or are replaced by values of Vi and Dij from the previous time step. Using Equation (4.1.1), C and qc can be written A M c = nil Nn cn = [Nn]{Cn} (4.3.2) 4 M . qC = nil N (qc)n = [Nn]{(qc)n} (4.3.3) In order that C be an exact solution of Equation (3.5.3), Equation (3.5.3) must be identically zero when C is sub- stituted into it. To minimize the errors of residuals as discussed in Section 4.1, the orthogonality condition requires that IDe R(C) Nk dD = [De R([Nn]{Cn})Nk dD= 0 (4.3.4) 47 where De represents the integral over the domain of an element. Substituting Equation (4.3.1) into Equation (4.3.4) yields A A ac -AIE‘if _ D? 3t+ 7(c vi ) axi [Dij axj qC Nk dD — 0 (4.3.5) By use of the Green theorem, A 3N A _8-[D. £INdD='I D. ._k_a.<'.dD [De {axi 1] axj k De 1] sxi axj . BC“: + [Se Nk Dij axj Ni ds (4.3.6) where Ni is the directional cosine of the boundary at the node under consideration. Substituting Equation (4.3.6) into Equation (4.3.5), one arrives at A A DN A 8C a - k 3C —— + ——— (C vi ) — q JN dD + DI. ——— ——— dD [De 8t axi c k De 1j axi axj —[ N D3. 3C 4. as (4.3.7) Se k 1 8x.1 j 1 Expanding Equation (4.3.7) and substituting Equations BNn + ViCnNk 5;: dB (4.3.2) and (4.3.3) yields BCn avi N N ——— dD + C N N ——— De n k at De n n k 8X 1 48 SN 3N A + D' k ——E C dD = N DI. 39— R. as e . e X . D S J + JDe (qc)n NnNk dD (4.3.8) As discussed in Section 3.5.4, the general equation for the boundary condition can be written as 3C ' —— = alDij ij 2i + NZC + d3 0 (3.5.4) If ul = O, C = a; = f(xiL provided d2 # O This boundary is the "Dirichlet" boundary. For nonzero al, the Galerkin formulation of Equation (3.5.4) will be A G’— e Nk Di. 3: 2i as = e Nk - c - W3. as = s 3 3 s 1 1 A 0‘2 0‘3 - e Nk C 5‘ as - e Nk 5‘ as (4.3.9) 5 1 s 1 QIQ N Substituting Equation (4.3.9) into Equation (4.3.8) yields _ DCn avi 3Nn e Nn Nk 8t dD + e Cn Nn Nk 8x. + viCnNk 8x. D D D 1 1 49 a3 . - Se Nk —— dS + De (qc)n Nn Nk dD (4.3.10) Equation (4.3.10) may be written in matrix form [Cheng 1973] T IHIe{%E—} + ([K]e+[S]e){C} + [Ele IIIIC} = [H]e{éc} + {F}e (4.3.11) where [H]e = H: n = I e Nk Nn dD ' (4.3.12a) I D e e 8Nk 3Nn [K] = Kk n = e 1311.? 5.)?— dD (4.3.12b) ’ D 3 i ' e e avi 3Nn [S] = Sk n = e NnNk §§_ + ViNk ax dD (4.3.12c) ’ D i i {E}e = Ee = 33-N dS d # O (4 3 12d) k Se a1 k 1 ‘ ° { }e _ e _ 0L3 F _— Fk — — e d— Nk d8 d1 # 0 (4.3.12e) S 1 If d1 = 0, then Ek = Fk = O. T [I] is the identity matrix, and [E]e is the transpose of the vector {E}e. F1".- ' 50 It is assumed that the velocity vectors are known at the nodes and that the dispersion coefficients are calculated for each integration point using Equation (3.5.1). Calculation of Equation (4.3.12) for different types of elements is given in Appendix II. Upon evaluation of Equation (4.3.11) for all elements, they are assembled to obtain the global relations [Hug-E} + ([K]+[S]){c} + [E]T[I]{C} = [H]{<§1c}+ {F} (4.3.13) [H] and [K] are banded symmetric matrices, while [S] is a nonsymmetric matrix. 4.4 Finite Element Computation of Velocity Vectors Mathematical equations of velocity vectors are discussed in Section 3.2. When the coordinate axes coin- cide with the principal direction of the hydraulic con— ductivity matrix, the Darcy equation can be written: Kij 23¢ Vi = - H—_ 5;? 1,3=l,2,3 ; Kij=0 when 1%] e 3 (4.4.1) In two-dimensional Cartesian coordinates, K _ _1_1 ELL V1 — n 8x (4.4.2a) .Velocity components play an important role in predicting the tracer movement in a porous medium, because they appear in the convective term and are used in the calcu— lation of dispersion coefficients. Special attention thus has to be given for the evaluation of the velocity vector. Two techniques are discussed below: 1. Direct calculation in this study is defined as the technique in which the velocity vectors are cal— culated using only the gradient of the shape functions multiplied by the corresponding piezometric heads, ele- ment by element, and is outlined in Section 4.4.1. 2. Simultaneous calculation is defined to be a procedure where the continuity of the velocity vectors is maintained and the gradient of the piezometric heads is multiplied by a weighted coefficient, as it will be shown in Section 4.4.2. The conjugate function concept is also used when referring to a similar method. 4.4.1 Direct Calculation One of the most common methods to calculate the velocity vectors is to substitute Equation (4.1.1) into Equation (4.4.2). The resulting equation will be _ ___ ¢ (4.4.3b) 2 ne 8x2 where ¢n's are piezometric heads at the nodes and are known. In Equation (4.4.3), BNn/axl and BNn/sz are first derivatives of the shape functions and are evaluated at the point of interest. Subscripts for x represent the coordinates, and subscripts for shape functions represent the node numbers. For example, for the element in Figure 4—2 with four nodes, V1 and V2 at point A are 4 3 X2 2 Xll Figure 4—2.—-A Typical Finite Element. K EN EN EN EN V1l = ' l a l ¢1'*a 2 d’2 a $3 a ¢4 A ne Xl X Xl X V =_KE_NI¢+W_2¢+E“_3¢+3N4¢ 2‘A 113 8x2 1 8x2 2 8x 3 8x 4 53 with the BNn/Bxl and BNn/ax2 computedau:point A. This procedure is easy to apply and provides accurate results at the centroid of the element, but the velocity compo- nents become discontinuous along the element edge as will now be demonstrated. Consider linear rectangular elements as shown in Figure 4—3a. The piezometric head ¢ and its first deriva- tive with respect max along lines AB and BZC are plotted l 1 in Figure 4-3. Figure 4—3b shows that the piezometric heads are continuous between two elements, and Figure 4-3C demonstrates how Sch/3xl is discontinuous between the elements. The gradient of ¢ takes different values depending on which element is used to compute its value. The discontinuity of the first or higher derivatives is also discussed by Zienkiewicz [1971]. To overcome the difficulties associated with discontinuity of the first derivative at the nodes, different techniques such as the conjugate function concept [Gallagher 1975] and Hermitian shape function [Zienkiewicz 1971] are used. In the next section the technique of simultaneous calculation of the velocity vectors at the nodes is pre— sented. 4.4.2 Simultaneous Calcula- tion of Velocity Vectors In order to provide continuous velocity functions between elements at the nodes, a smoothing technique, 54 9‘2 Axl 41'.— B B C X1 A l 2 (a) II II CD II [I II II I' $Xl I' (b) II II 23¢ IL “1 I,II III II lI l1 , ==x " (c) 1 ll Figure 4-3.-—Interelement Zone Depicting How Continuous Function ¢ May Have Discontinuous Gradients as Axl + O. 55 namely, the Galerkin-based finite element formulation of Equations (4.4.2), is developed. The concept is similar to that presented by Zienkiewicz [1971, pp. 44- 46] and also used by Segol et a1. [1975]. In the pro- posed technique which is outlined below, the computer storage and the complication of the solution process are reduced considerably. Let Vi, VS represent the simultaneous calculation of V1 and V2. The piezometric head and velocity vectors can be written in terms of nodal parameters. M 4 = [m1 {4} = nil Nn 4n (4.4.4a) GO - [N {vc} - § N (VC) (4 4 4b) 1‘ 11“_ n 1n '° n—l 9C - [N1{vc} — § N (vc) (4 4 4 1 2" 2‘n=1n 2n "0 Equation (4.4.4) will assure that continuity of the vari- ables is maintained along the element boundaries, regard- less of the type of elements used (by definition of the shape function). The technique of calculating the velocity vector along the x -coordinate is presented 1 below. The residual of Equation (4.4.2a) is K A _ “c 11 So 1) ‘ V1 + n ax e 1 (4.4.5) 56 The definition of V is given by Equation (4.4.2a) l where K and ¢ are known parameters. Substi- 11' ne’ tuting Equation (4.4.5) into Equation (4.1.4) yields { N 6C + El— 3$ dD = o (4 4 6) De k l ne 5x1 A K A c 11 3¢ _ D D e 1 Introducing Equation (4.4.4) into Equation (4.4.7) gives Kll BNn C ""— —_.._ {De Nan(Vl)ndD — [De Nk n 3x ¢ dD (4.4.8) e l n Since (V$)n is independent of the space coordinate, it can be taken out of the integral c K11 8Nn {(vl)n} De NandD = — 5;_ De Nk SE“ ¢n dD (4.4.9) The right-hand side of Equation (4.4.9) is known and becomes a column matrix. Thus [H]e{\3)(>\->\4) 1 3 (ll-12)(Al-l3)(ll-l4) -6(Al) N = (A—Al)(A-13)(A-X4) : (A-Al)(A—l3)(A-A4) 2 3 (42-41)(12-13)(12-14) 2(A1) N = (A-Al)(A-12)(A-A4) : (A-Al)(A-12)(A-A4) 3 3 (AB—ll)(l3-12)(A3-A4) —2(Al) N = (l-Al)(1-12)(A-l3) = (l-Al)(A-12)(A-l4) 4 3 (A4-Al)(A4-A2)(A4—A3) 6(AA) Define the variable 9: A - A 9 = 3 l4 - 13 (5.5.1) (5.5.2) 74 such that 9 = O for A = A3 and 0 g 6 IA [.1 The shape functions in terms of 6 can be written _ 1 N1 — g-[9(9+1)(l'9)] N = i [e(e+2)(e-1)] 2 2 _ 1 N3 — f [(e+1)(e+2)(1-e)] N4 = % [e(e+1)(e+2)] (5.5.3) The value of C and its first and second derivatives with reSpect to A can be written in the following form: C = alC4 + d2C3 + d3C2 + a4Cl (5.5.4) 39 = 3L [8 c + B c + B c + B C J (5 5 5) 3A AA 1 4 2 3 3 2 4 1 ° ° 2 3—9 = l [y C -+y C -+y C -+y C ] (5 5 6) 3A2 (AA)2 1 4 2 3 3 2 4 1 ° ° The values of a, B, and y in terms of 8 are given in Table 5—3. 75 TABLE 5-3.--Va1ues of a, B and Y for Third Order Approximation in Terms of 8. n a B Y n n n 1 % e(e+1)(1—e) % (1-392) -e 2 % e(e+2)(e-1) % (392+2e-2) (3e+1) 3 % (9+1)(e+2)(1-e) -% (362+46—1) -(3e+2) 4 % e(e+1)(e+2) % (382+68+2) (9+1) The variations of a and B with respect to 6 are given in Figure 5-6. For a time-dependent variable, Equations (5.5.4) through (5.5.6) can be written C = alC(t+At) + a2C(t) + a3C(t—At) + a4C(t-2At) (5.5.7) 8C __ 1 _ 5E - A? [BlC(t+At) + 82C(t) + B3C(t At) + B4C(t-2At)] (5.5.8) and 82c 1 ——— = [Y C(t+At) + y C(t) + y C(t—At) atz (At)2 1 2 3 + y4C(t-2At)] (5.5.9) 76 V (D Figure 5-6.--Variation of d and 8 With 6 for Third Order Approximation. 77 To obtain a third order recurrence formula for a set of first order partial differential equations, Equations (5.5.7) and (5.5.8) need to be substituted in Equation (5.1.1). The coefficients of Equation (5.1.5) for dif— ferent 6 for third order time approximation are given in Table 5-4. 5.6 Summary In this chapter the recurrence formula for the first, second, and third order approximations for a system of first order differential equations is derived, and the procedure of obtaining the finite difference relations using the finite element concept is presented. In this study the first order approximation was used in the solution of the flow equation, while the first and second orders were employed for solving the convective- dispersion equation. The third order was not examined in this work but for sake of completeness it was intro- duced. It will be shown in Chapters VIII and IX that 8 = 2/3 provides less oscillatory and more reasonable results for the flow and mass transport equations. However, choosing the value of 8 depends on the nature of the equation, the numerical technique, and the required use of previous known values. Giving any specific value for 8 at this stage would be premature. More work needs to be done in this area. 78 m m- m) u H o o m o o m o H HH H H m .E B pm- E- B... .E E E E. .s. E E Mus m H m em m «N om hm ow ow mm ow N all. .m- E. E- m E me 4| Ms. B .4) me Mums m H m mH m mH mw HN mv mH mm mH H Nam Ham mmm mmm Hmm mmm mmm Hmm mam Nam Ham Ame 801562 .m ncmummwfle How GOHDMEHxOHmm< wEHB Hmpuo UHHQB “Ow Am.H.mv COHumsqm mo mpcmHOHmmooonl.vlm mqm(t + 446)} = {6} (6.1.1) where [E] - [B] + a” [H] (6 1 2) " a11 —A_t_ ° ° and {G} - a [B]+EEE[H] + a {F(t+At)} ‘ 21 At 13 + a23{F(t)} (6.1.3) The values of the coefficients are given in Table 5-2. Typical boundary conditions are discussed in Sec- tion 3.3.3. From these conditions the lateral recharge Q2 as well as the Dirichlet boundary are assumed to be time invariant. For many applications, time variable pumping must be considered at various nodes throughout 81 the flow domain. In the computation it is convenient to handle the pumping as a series of step functions dis- cretized with respect to time as shown in Figure 6-1. Thus at any time interval recharge or discharge will be constant. There are different numerical techniques which can be employed to solve the system of Equations (6.1.1). Since [E] is a banded symmetric matrix with nonzero terms in the diagonal, it is possible to use Cholesky's square root procedure to decompose the matrix [8] (upper band) and solve with the companion subroutine for the unknown I{¢} as outlined by Weaver [1967], and employed by Pinder and Frind [1972] for groundwater flow. The required storage space for [E] in the computer will reduce to NNDS X MAXBW, where NNDS is the number of unknowns {4} and MAXBW is equal to time upper bandwidth plus one. The parameters such as transmissibility, storage coefficient, and time step participate in the construction of the [8] matrix. Thus if one desires to change any of the above mentioned parameters, it will be necessary to regenerate the [E] matrix again. Equation (6.1.1) can also be considered to be a "steady-state" problem if the coefficients for the steady— state condition from Table 5—2 are used. The resulting equation for this state is 82 ' '0 Representative discharge curve / ‘4 \\_;’/,/’7 | I ' I I l I l m l I U1 l l H I g l l l o l Step’J I -2 ' r‘Changel I I Q l l l l | I II . , I I I I I I l J L 1 :7 t Time Figure 6—l.--Dividing the Actual Discharge Curve Into a Series of Step Functions. 83 [B]{¢} = {F} (6.1.4) 6.1.2 Introducing the fiIrichlet Boundary Condition The placement of the Dirichlet boundary condition into the global equations can be accomplished by the deletion of rows and columns [Norrie and de Vries 1973]. Assume in Equation (6.1.1) that the piezometric head at node k is known. Equation (6.1.1) in expanded form can be written K \ \ K \ 13 E 15 rd) G Ekl . . . Bkk . . . Bkn J ¢k >=< Gk f (6.1.la) EDI ml CHI '9- Q nk nn n n Since ¢k is known, all coefficients in the [E] matrix at the kth row can be set equal to zero except the diagonal terms (i.e., gkk) which will remain unchanged. The Gk is replaced by Ekk¢k’ Equation (6.1.la) will have the following terms: 84 o . . . Bkk . . . o <¢k >=< Bkk¢k >(6.1.5) Bnl ° ° Bnk ' ' Bnn ¢n Gn Usually Ekk has a positive value. In the computer pro- gram a check has to be made to be certain that the value of E k does not fall below a certain small value, e.g., k one. Otherwise Ekk can be changed to any large value. In the computer program, Ekk is replaced by the average value of the diagonal terms. Equation (6.1.5) can handle the Dirichlet boundary condition, but the matrix is then banded and nonsymmetric. Since Cholesky's method is used to solve the system of equations, the matrix of Equation (6.1.5) must be symmetric. This can be accomplished by subtracting Eik¢k7 k, i= 1. . . n, i # k from both sides of Equation (6.1. ). The final form becomes 85 K \ /\ / \ B11 ' ° ' 0 ° ° ‘ B1n ¢1 G1‘B1k¢k o . . . Bkk . . . o < ¢k >==< Bkk¢k >(6.1.6) Bnl . . . o . . . Bnn 4n Gn—Bnk¢k \ —--—-J \/ J\ / This procedure is carried out for all Specified piezo- metric heads prior to decomposition of the [E] matrix. 6.2 Solution of Flow Vectors The mathematical equation of the velocity vectors is discussed in Section 3.2 and the methods of solution are described in Section 4.4. Once the piezometric heads have been determined, flow (or velocity) vectors can be computed. 6.2.1 Direct Calculation In the direct method, Equation (4.4.3), i.e., Kij aNn Viz-T725..— (bi l,j=l,2 (4.4.3) e 3 is employed to compute the velocity vectors. This tech— nique provides accurate results at the centroid of each element, but the velocity components calculated by this 86 procedure are discontinuous across element boundaries. The magnitude and direction of the velocity vectors (Figure 6-2) can be obtained by < ll j 1 2 9 ll arctan (VZ/vl) 6.2.2 Simultaneous Calculation of Velocity Vectors IV I = /V2 + V2 j=1,2 (6.2.1) Simultaneous calculation provides continuous velocity components at the nodes. The finite element formulation of this technique is presented in Section C 4.4.2. For example, for V1 the element equations are K 8N =[ 111.4,... e xn D The global matrix is [HHV‘l’I = {F } 1 Integrated forms of Equation (4.4.11) and feaN D (4.4.11) (4.4.12) (4.4.13) (aNn/axl)dD for different types of elements are given in Appendix II. A similar global relation is obtained for V 87 ‘ r—xl Figure 6—2.--Velocity Vectors in the x x Plane. 1' 2 C [H]{V2} = {F2} (4.4.14) The [H] matrix is a banded symmetric matrix. Cholesky's square root procedure (see Section 6.1.1) is used to decompose the matrix. As long as the nodal coordinates are fixed, the [H] matrix can be decomposed and stored with no need to recalculate it. If piezometric heads and velocity vectors were to be calculated at the same time in a computing scheme, the size of both the [H] and [B] matrices would increase up to 3NNDS X 3NNDS with many zero terms, NNDS being the number of nodes in the system. This procedure requires a great amount of computer 88 memory, and the transient solution of the problem is tedious. In this method, besides the boundary conditions for the piezometric heads the conditions for the velocity vectors need to be specified. But the procedure out- lined in the preceding paragraph, based on the simple realistic assumption that the piezometric heads can be calculated independently of the velocity vectors, pro— vides great simplification. With the construction of the [H] matrix which requires only MAXBW x NNDS storage core, the velocity vectors can be calculated at the nodes. Equation (4.4.12) can be calculated for each element and added to {F1}, a column matrix with NNDS rows. Similarly, '{FZ} will be evaluated. The [H] matrix for Vi and v: is identical. Thus the solution of Equation (4.4.13) and Equation (4.4.14) will provide the velocity vectors at each node. Equation (6.2.1) can be used to obtain the magnitude of the velocity and its direction at each node. Known velocities at the node can be handled by deletion of rows and columns, as outlined in Section 6.1.2. In Section 8.2 the two methods for calculation of velocity vectors are compared. 6.3 Solution of System Equation for Flow in Unconfined Aquifer With Transient Phreatic Surface In this study the flow in porous media on a regional as well as local scale is studied. For a 89 regional scale the two-dimensional horizontal flow is considered and the solution of related equations is given in Section 6.1. For a localized scale a vertical cross section of an unconfined aquifer is chosen and the solu— tion of the system of equations is presented in this section. A technique of locating the phreatic surface with fixed nodes is also shown. 6.3.1 Background Based on the assumption made in ChapterIELI,if the specific storage coefficient can be neglected, Equation (3.4.1) will be reduced to 8 8¢ 332‘; ij 53?; = 0 i,j=l,3 (6.3.1) The finite element formulation of Equation (6.3.1) leads to the global matrix of the form [B']{¢} = {F} (6.3.2) with its element components defined as e BNk aNn [B 1 : Kij ST WdD (6.3.3a) De 1 3 and e {F} — ~[e Nk Q2 dS (6.3.3b) 90 Equation (6.3.2) is solved with the related boundary conditions discussed in Section 3.4.3. In addition to the known potential and known flux boundary, there is the phreatic boundary condition which is described as At Uj-zi = UnAt.= ——- -K j2£_2. - I 23] i,j=l,3 (3.4.6) ne ij ij 1 Equation (3.4.6) contains a time derivative of the free surface and can therefore be used to determine the height at the later time when the other terms in the equation are known [France 1971]. An iterative technique is used to replace the original transient problem by a discrete number of steady-state problems, based on the assumption that the flow at each instant is steady but the boundary of the flow is time variable [Poluborinova- Kochina 1962, p. 572]. At the beginning of each time interval the position of the free surface and boundary conditions are known. Using Equation (3.4.6), the phre- atic surface is then shifted along its normal to a new position. Two possible methods are considered below. One technique most commonly used (e.g., Desai 1972 and France 1974) requires the modification of the elements such that the free surface always is the upper boundary of the grid system. Another method, presented by France [1971] and modified and improved herein, accomplishes the movement of the phreatic 91 surface within the grid system without repositioning the nodal coordinates of the elements. 6.3.2 Locating the Phreatic Surface by Modifying the Elements The movement of the phreatic surface is deter- mined in the following manner: Step 1. At the beginning of each time interval the surface configuration and boundary conditions are known. Equation (6.3.2) is used to find the value of the piezometric heads for all nodal points except those at the phreatic surface. If there is a known flux from boundaries other than the phreatic surface, the values of {F} in Equation (6.3.2) are calculated; otherwise {F} will be zero. At this step the nodal points at the phreatic surface act as the Dirichlet boundary. In Sec- tion 6.1.3 it is shown how to introduce the known potential into the global equations. Equation (6.3.2) is similar to Equation (6.1.4) and is solved the same way, as outlined in Section 6.1.1. Step 2. At the beginning of each time the piezo- metric heads are known in the system, and from Equation (3.4.6) it is possible to compute the location of the phreatic surface at a time t + At. The distance a point on the phreatic surface will propagate in the direction of the normal to this surface at that point 92 is equal to UnAt. If P(xl,x3) in Figure 6-3 is a point of the phreatic surface at time t and P'(xi,x5) is the location of that point at time t + At, then the shifting distance along the normal to the free surface dn’ PP' in Figure 6-3, is _— —— — I: __ — —_ — . + . . A where i and k are unit vectors along x and x3, respec- 1 tively; and £1 and £3 are directional cosines. I is the accretion term, positive downward. If 6 is the angle that the tangent to the phreatic surface makes with the positive x —direction, then 1 2o ll 1 sin 6 & 2 = cos 8 dn = At (Vl Sln 8 + V3 cos 8 - I cos e/ne) (6.3.5) Thus dn can be evaluated by Equation (6.3.5). It is Convenient to calculate the shifting distance along the nodal lijxs (d2). Let mi be the angle that the nodal line (i) makes with the positive x -direction, and 1 define B = N/2 - w + e [Desai 1972]. Then d n dg = cos 8 (6.3.6) 93 .oomwusm OHDNmHSm map mo ucoEm>ozcl.muo musmHm N lmxc.mlnxc Himxv.slaxv Hx I FE:C\\\\:\::\:\:\::C:\§\\\:: N3 H3 mwcHH Hmwoz mommHsm \ UHDNmHne \x.\ U. <+ U QEHU. My“ CD. HMEOZ \\\ mommnsm 0Hummucm///V. me.MMYm I m. .1 \\||\.\\ aw mm. D wEHD pm n H a wommusm UHumecm me. xvm U m unwom pm ucwmcme 94 d = d cos w (6.3.7a) D; II 3 d2 Sln w (6.3.7b) In Equation (6.3.7), d and d are the shifting distances l 3 in the x1 and x3 directions, respectively. Summarizing, (V' sin 6 +rV' cos 6 - I cos G/h.)At cos w _ l 3 e d1 - (6.3.8a) ast (Vi sin 8 + V: cos 6 - I cos e/ne)At sin w d3 = (6.3.8b) cmsB Thus the location of the phreatic surface at time t + At will be x = x + d (6.3.9a) X = X3 + (13 (6.3.9b) In order to complete step 2, one has to know the velocity components and the angle 6 at each phreatic node. As discussed in Section 6.2, direct calculation of the velocity vectors usually does not provide continuous results at the nodes. Most authors, e.g.,France et a1. [1971] and Desai [1973], have realized this deficiency and have used the average velocity components at the node calculated from two adjacent elements. Although this procedure reduces the errors significantly, because of 95 discontinuity of the velocity functions it is still not the best way to evaluate the velocity vectors at the nodes (see Section 4.4.1 for details). The proposed simultaneous calculation of the velocity components pro— vides continuity at the nodes and is used in this study. The slope of the element lines along the phreatic boundary is usually discontinuous at the nodes. For example, at node 8 the slope of line 4-8 in Figure 6-4 differs from line 8—12. In order to obtain a better estimate of 6 at the phreatic nodes, a polynomial of degree n is passed through the phreatic nodes of two elements while the node under consideration is almost Phreatic surface (9 8X " 3\Q o 14 7 11\ GD .15 2—~——~_____~______d 6 10““14 .3 <9 69 o l 5 9 13 X1 Figure 6-4.-—Simp1e Linear Quadrilateral Isoparametric Elements. 96 the middle point. When the polynomial is passed through the nodes the slope of the free surface is evaluated at the node. The value of n is dependent on the type of element. For linear and quadratic quadrilateral elements n = 2 will be sufficient, while for a cubic element n = 3 is recommended. For corner nodes (e.g., nodes 4 and 16 in Figure 6-4) the slope of the phreatic line is used. Because the phreatic nodes might not be equally spaced, Newton‘s divided-difference method [see Carnahan et a1. 1969, pp. 9—26] is used to pass a polynomial for the desired points. Step 3. Equation (6.3.8) assures that the shifted points always lie along the nodal line. In this case, Equation (6.3.9) represents the location of the new nodes at a new time. In some instances the nodal lines cannot be straight and the value of w changes along the line, such as for quadratic or cubic elements. Then * * the new point xl,x will be calculated by 3 * - xl — xl - dn Sln e (6.3.10a) * x3 — x3 — dn cos 6 (6.3.10b) * * where x1 and x3 represent the temporary location of the nodes at time t + At, and dn is defined in Figure 6-3. To find the location of the actual nodes, it is necessary to fit a polynomial to these temporary nodes and then . ._h 97 find the location of new points with the given xl-values (Figure 6-5). Polynomial curve-fitting procedures such.as Newton's divided—difference method (Carnahan et a1. 1969], or other methods such as the Newton-Raphson iteration technique or the least square curve-fitting [Pennington 1970, pp. 408-417], can be used. In this study the first method is employed, and the second technique is used by France et al. [1971]. Step 4. At the start of this step the location of the phreatic surface is known and can act as the Dirichlet boundary condition for solving Equation (6.3.2) to compute the piezometric heads in domain D. But since the free surface has moved, the nodes on the phreatic Phreatic surface at time t + At Phreatic surface at time t 0 Temporary nodes 0 Actual nodes Figure 6-5.——Location of the Temporary and Actual Nodes [After France 1971]. 98 surface do not coincide with the nodes of elements. There are two possible ways to handle this problem: (a) modify the elements such that the phre- atic surface becomes the upper boundary of the grid system, or (b define a new Dirichlet boundary condi- tion for the nodes above the phreatic surface such that it is not necessary to modify the elements. v Procedure (a) is discussed below, and details con— cerning procedure (b) are postponed until Section 6.3.3. The location of the phreatic surface is known, so the upper nodes of the phreatic elements are transformed to coincide with the phreatic surface. Since the coor- dinates of some nodes have changed, it is necessary to re-evaluate the [B'] matrix (Equation 6.3.3a). In order to reduce computation time the whole grid system can be divided into two groups (Figure 6-6): (a) fixed elements--their nodal coordinates will not alter during the entire calcu- lation, moveable elements—-their nodal coordi— nates will be affected by changing the phreatic surface. (b v Another important point is that the rise or fall of the phreatic surface is not uniform, and in some nodes the phreatic surface moves several times more than other nodes. In order to keep the elements in reasonable shape it is recommended that all moveable elements be modified rather than only shifting the nodes of phreatic 99 .EwHQOHm HOMHDw¢ OHDMOHSm mcH>Hom MOM Empmmm pfluw HMOHmmBII.mIm musmflm mopoc oHQm>OS O mmpoc ~5me . ® 3 :4 44 .- .. . .. ® ® ® O $33 3me $3 3.01 Z 3 m N G O G G O mm Saw 30 3G 5 O m J O G ® ® ® 4N om emu «we m.Lv v v 1:23:23. Acohumnnafluqflc H 100 elements. When the nodal points are shifted the program is ready to start a new time step. Steps 1 to 4 are repeated until the system reaches the maximum specified time or the steady state. 6.3.3 Location of Phreatic Surface Using Fixed Elements As discussed in Section 6.3.2, it is possible to handle the movement of the phreatic surface without altering the nodal coordinates of the elements. In the following method the phreatic surface travels within a fixed finite element grid system. 6.3.3.1 Computation of piezometric heads at the nodes above phreatic surface.——In order to illustrate the method, a plane linear element, Figure 6—7, is employed. The first part of this development is adopted from France [1971]. Referring to Figure 6—7 which simu- lates a small portion of the flow domain, the phreatic surface represented by the broken line cuts through ele- ments numbered (1), (3), and (4) at points a, b, c, and d. For element number (1) the piezometric head ¢ at any point within or along its boundaries is given by: 4 = N141 + N242 + N5¢5 + N4¢4 (6.3.11) 101 1 2 3 Q) n *__a_____\. o 1A- \‘\H T— ‘TDKLB: 5 ~\ \‘\.ET1:H1 \\ 4 5 C .\ 6 ‘~d 1; ” :1 E 6 x3 (:2 (:2 9 Figure 6-7.--PhreatL:Surface Passing Through the Elements [After France 1971]. The shape function for a plane linear element can be written [Zienkiewiez 1971, p. 109] 1 4 ) 2 ll (1+ 60) (l +n0 where 50 E Ej and no = n n. J For element number (1) the shape functions are N1 = % (1-6)(1+n) N2 = % (l+€)(l+n) N—ll- 1 N-ll+ 1 4 — 4 ( €)( -n> 5 — 3 ( €>< -n> (6.3.12) 102 At point (a) along line 1-4 (5 = -l), the piezometric head is _.1 4a — Z-{2(1+A)¢l + (0)(1+A)¢2 + 2(1-.A)¢4 + (0)(1HA)¢5} (6.13.13) where A is the value of n at point (a), and —l g A g 1. Finally, -1 _ 4a — 7 {(1+A)l + (1 A)¢4} (6.3.14) Similarly for point (b) along line 2-5, 4b = % {(1+B)¢2 + (l-B)¢5}, -1 g B 5 1 (6.3.15) For element number (3) the piezometric head is given as ¢ = N2¢2 + N3¢3 + N6¢6 + N5¢5 (6.3.16) At point (c) on line 5-6 the head is 40 = % {(1+C)6 + (l-C)¢5}, -1 g c g 1 (6.3.17) For element number (4) the head is given by ¢ = N5¢5 + N6¢6 + N9¢9 + N8¢8 (6.3.18) At point (d) on line 6-9, ¢d = % {(1+D)q>6 + (1—D)¢9}, -1 g D g 1 (6.3.19) 103 The objective is to determine the piezometric head distribution on and below the phreatic surface. How- ever, the finite element method yields values only at the nodal points, including those points above the phre- atic surface. It is therefore necessary to define the values of ¢l’ ¢2, ¢3, and ¢6 in terms of 0a, ¢b’ and ¢d‘ Since ¢3 does not affect values along 2—5 and 5—6, it can be determined arbitrarily as will be discussed later. The piezometric head at node 6, ¢6’ can be Speci- fied either in terms of ¢c or ¢d' From Equatfbns (6.3.14), (6.3.15), and (6.3.19); ¢l' ¢2, and $6 can be derived: 41 I24a - (l-A)¢4]/(1+A) e- M II [24b — (l-B)¢5]/(1+B) (6.3.20) 46 = [24d - (1—D)491/(1+D) It is necessary to know the values of ¢a’ ¢b’ and ¢d‘ This presents no difficulty since on the phreatic surface the piezometric head equals the elevation head. Thus Equation (6.3.20) can be used to interpolate the values of the piezometric heads at the nodes above the phreatic Surface, based on the location of the free surface and known piezometric head at the nodes beneath the surface. Equation (6.3.20) which has been derived by France [1971] 104 is simply a linear interpolation between three points where the distance between the nodes and the values of two nodes are known. It is desired to find the unknown value at the third point. using a one—dimensional, two—node element. This can easily be shown 01 Figure 6-8.--One—Dimensional Linear Element. In Figure 6—8, let $1 and ¢A be known with the intent to find ¢2. The distances between the nodes are defined in the above figure. Shape functions for nodes 1 and 2 are _ 12"LA LA N1 _ L ’ N _—_ 12 12 By definition, ¢A = Nl¢l + N2¢2 (6.3.21a) or 105 (>2 = (4A - N1l)/N2 (6.3.21b) Substituting the shape functions into Equation (6.3.21) yields L ¢ ' (L -L )4 _ 12 A 12 A 1 $2 — L (6.3.22) A . 2LA Define B = E—_ — 1 such that —l g B 5 l. 12 Substituting A in Equation (6.3.22) and simplifying, (1) _ 2¢A - (l - B)¢l , 2 - W (6.3.23) Equation (6.3.23) is identical to Equation (6.3.20) which has been derived using the properties of an iso- parametric element. Usually, ¢ represents the piezo— A metric head at the phreatic surface, $1 is a known head within the system, and calculated $2 is a Dirichlet boundary condition. The problem associated with Equation (6.3.21b) or Equation (6.3.23) is that it is singular when N2 + 0 or B + -l, i.e., when the point A is close to point 1. In practice, when the phreatic surface rises due to infiltration or other hydraulic stresses, Equation (6.3.23) does not provide a reasonable estimate for $2. 106 One way to obtain a better estimate for phreatic nodes (Figure 6-10) is to use two interior nodes besides ¢A to evaluate unknown ¢. Let ¢l and ¢2 be known and ¢A represent the value of the piezometric head at the phreatic surface. It is desired to find 03 (Figure 6—9). Using the properties of a one—dimensional quadratic ele— ment (see Section 5.4), ¢A ' N1¢1 ' N2‘1’2 ¢3 = ——-—-—-fi;—-————— ; N3 # 0 (6.3.24) where N (LA )(LA L13) 1 (L12 (L13) (L )(LA - L13) N2“ () ) 12 23 L13 ’3 L13IIA I .2 L 1L .1_ 12 A 1.. I Figure 6—9.--One-Dimensional Quadratic Element. 107 I (accretion) 9 * H++i /Nodes above the 8 ,‘ phreatic nodes I \) 4 Phreatic nodes —\ g~H~¥1\\~‘ r/Phreatic surface Y Interior nodes / I Nodal line Figure 6-10.——Location of the Phreatic Surface and Definition of the Terms. 108 3 (L13) ”‘23) Although Equation (6.3.24) provides justifiable results it still becomes singular when N3 + 0. A possi- ble way to compensate for the singularity of N2 in Equation (6.3.21) or N3 in Equation (6.3.24) is to use the node below point 1 of Figure 6-8, when point A is close to point 1. The scheme will be similar to Figure 6+9 but ¢ - N'¢ 43 = A l 1 (6.3.25) NI 2 where ._ L13 - LA ._ LA N1 ’ L ' N2 ‘ L 13 13 Equation (6.3.25) will never be singular, but since it is using $1 rather than ¢2 it underestimates ¢3 because usually $1 < $2- Returning to Figure 6-8, there is a third way to compute $2 and it is based on Darcy's Law: V3="ni3’—a:¢ e 3 or V a.) __ 3 n (42 - 4,) v3 “‘3’2 ' (X3)A K33 9 Since (x3)A E 0A, then V n _ _ 3 e _ (6.3.26) ¢2”¢A K (”3’2 1’21) 33 V3 and K33 are known values at point A. In Equation (6.3.26), $2 is calculated based on the location of the phreatic surface. If the values of are evaluated correctly, then it is believed V and K 3 33 that Equation (6.3.26) will provide accurate results for ¢2 without under or overestimating its values. When the phreatic surface is horizontal, V3 is zero and thus Equation (6.3.26) is not applicable. In this case, for the first time interval Equation (6.3.25) is employed to compute the value of the piezometric head at the points above the phreatic surface. For the second time interval and later times, one will instead use Equation (6.3.26). It can be concluded that among the Equations (6.3.23) through (6.3.26), Equation (6.3.25) gives a reasonable value at the first time interval and Equation (6.3.26) provides a better estimate at later times. 6.3.3.2 Computation of velocity vectors at phreatic surface.--When calculating the shifting distance for the phreatic surface (Equation 6.3.8), itis necessary 110 to know the velocity components on that surface. Since the phreatic surface moves within fixed nodes, in most cases it will intersect the element line as shown in Figure 6-7. The velocity components are known at the nodes. Consider again Figure 6-8 to calculate V1 and V3 at point A between nodes 1 and 2. Using an equation similar to Equation (6.3.21a) to compute velocity vectors, one will have (l—B) (v + (1+B) (v1)2 ) _ 11 (V1)A — 2 (6.3.27) where 2L B=L—A—l 12 LA /I(xl)l—(xl)A12 + [(x3)l-(x3)A]2 2 2 L12 = /1(xl)2—(xl)l] + [(x3)2-(x3)11 Similar equations can be written for (V3)A. 6.3.3.3 The procedure for solving piezometric heads with fixed nodes.--The procedure for solving the piezometric heads in an unconfined aquifer with this method is as follows: Step 1. The phreatic surface is initially chosen to coincide with the element boundaries. It is not lll necessary that this boundary be the upper limit of the grid system, i.e., the phreatic surface might be speci- fied within the system. This is especially important when a rise of the water table is expected. For example, in Figure 6-6 it is expected that the phreatic surface will rise, and line 2-6—10-14-18-22 can represent the free surface at the initial time. A finite element solution is performed, velocity components are calcula- ted, and the phreatic surface is shifted as in the fixed element method. Step 2. The piezometric head values at points such as a, b, and d of Figure 6-7 are stored. Their values are simply the elevation heads at the respective points. The piezometric head values at nodes 1, 2, and 6 are calculated using Equation (6.3.25), and an arbi— trary value is assigned to node 3 as will be discussed later. Step 3. Since the element stiffness matrices are unaltered, all that is necessary is to insert these new Dirichlet boundary values in the governing set of simultaneous equations and solve for a new set of piezo- metric heads. Step 4. The velocity components are evaluated and the phreatic surface is shifted as mentioned in Sec- tion 6.3.3.2, using Equation (6.3.5). 112 Steps 2, 3, and 4 are repeated until the maximum required time is achieved or until the system has reached the steady-state condition. For the second time interval or later times Equation (6.3.26) is used in step 2. 6.3.3.4 Assigning the value of piezometric heads above phreatic nodes.-—Assigning an arbitrary value of piezometric head to node 3 (Figure 6—7) or to similar nodes is an important task. To clarify the dis— cussion consider Figure 6-10 where the position of the water table is shown, and the phreatic nodes and the nodes above the phreatic nodes are defined. The piezo— metric heads at the phreatic nodes, i.e., nodes 3, 7, 11, and 15, are calculated using Equation (6.3.26) based on the location of the phreatic surface. The task now is to assign the values of the piezometric heads to nodes above the phreatic nodes, i.e., nodes 4, 8, 12, and 16. In the calculation of the piezometric heads within the system, these values do not play any major role. But they have a great importance particularly when the velocity vectors are calculated simultane— ously. It is assumed that there is no capillary fringe above the water table, that an abrupt interface exists between the saturated and unsaturated soil, and the .pressure above the phreatic surface is atmospheric. Due 113 to the physical characteristics of the porous media there is little or no flow movement above the phreatic surface, except the infiltration due to recharge which is assumed to travel directly to the water table. This concept, together with the above assumptions, directs one to assign the piezometric heads above the nodes such that the physical conditions are satisfied. This is done by equating the piezometric heads in the nodes above the phreatic nodes to the piezometric heads of the phreatic nodes along each nodal line. For example, $16 = ¢15 and ¢12 = ¢ll’ etc. This technique will provide the desired condition for velocity vectors along x3, i.e., V3 above the phreatic node will be zero, but there will be a gradient along x This situation will not cause any 1' major problem since the slope of the phreatic surface is small and its effect will be minimal in calculating the shifting distance (Equation 6.3.8). However, to improve the technique it is possible to equate the velocity com— ponents at and above the phreatic nodes (Figure 6—10) to zero, which is done in this study. Summarizing, the piezometric head at and above the phreatic nodes will be the same for each nodal line; refer to Figure 6—10. After the velocity components are calculated simultane- ously, the velocity vectors along x at and above the l phreatic nodes will be set to zero. 114 6.3.4 Steady-State Condition When the system has reached the steady-state condition the phreatic surface will not move considerably, i.e., the normal velocity will approach zero. From Equa- tion (3.4.6), -Kij gg? 2i - I 23 = o This means that at steady state, normal discharge is equal to the normal component.0faccreti0n- In the com— puter program Equation (6.3.5) is used. The average shifting distance is defined by NNFS k§1(dn)k d = NNFS where NNFS is the number of nodes in the phreatic surface. When a is less than or equal to a small value the steady— state condition is reached and the program is halted. 6.3.5 Reasonable Time Step There is no general agreement for choosing the time step. Rushton and Herbert [see France 1974] have suggested that if the time interval is chosen so that the velocity at which the phreatic surface moves changes by less than 30 percent between successive steps, then it iS not necessary to iterate within the time step. 115 Based on the Lipschitz criteria [Isaacson and Keller 1966, pp. 86-91], Sandhu et al. [1974] have intro- duced the variable time step. They state that in choosing At the following criteria should be satisfied: i+l i max 1(x3)k - (x3)k) = A(At) < 1 (6.3.28) i i-l where k = l, 2, ..., NNFS, i represents the number of iterations (see Section 6.3.6), and the value of A is dependent upon At. If A exceeds one, the time step will be reduced by some factor until the time interval satis- fies Equation (6.3.28). If, on the other hand, A is found to be very small compared with one, At will be increased by the factor. In the computer program (x3)£ is estimated, the value of the previous time is used, and only two computations of x3 are required to obtain an estimate for A and to choose the proper time step. 6.3.6 Iteration Within Time Step Usually, if one chooses a small time step there will not be any need to iterate within a time step, especially when the rise or fall of the phreatic surface is small compared to the size of the system, such as with f 116 actual field problems. However, Sandhu et a1. [1974], by using the mean value theorem and assuming a smooth change of geometry in the time domain, have introduced the following equations for iteration within the time step: i _ 1 + dl[xl(t+At), x3(t+At)]i—l}£ (6.3.29) [x3(t+At)]: = [x3(t)]k + %~{d2[xl(t), X3(t)] + d2[xl(t+At), x3(t+At)]i-1}L i and k are defined in Section 6.3.5. The inside bracket means that dl and d2 [see Equation (6.3.8)] are calcula— ted at [xl(t), x (t)] and [xl(t+At), X3(t+At)]. Since 3 the change of x3 is greater than the xl-value after every iteration, a check is made where i i—l [(x3)k - (x3)k ]/|(X3)k| is less than or equal to a (say, 10-5). If this condi- tion is satisfied the iteration is complete. Usually One or two iterations will be sufficient. The same tech- rnique is used in this study. 117 6.3.7 Comparison Between the Two Methods of Locating the Phreatic Surface Figures 6hll and 6-12 depict the schemes for solving the piezometric heads in an unconfined aquifer with the method of modifying elements (method one) and the method of fixed nodes (method two). It is extremely difficult to compare the two methods comprehensively. However, each technique has specific advantages as out- lined below. (a) In the first method the locations of nodal points are changing and the system is shrinking or expanding with time, in the second procedure the grid system is fixed and the coordinates of nodes are time invariant. The second method enables one to introduce convective-dispersion equations into the system, a major advantage of this technique. (b) In the first method it is necessary to reconstruct the global matrices and decompose them for every time step, while it is not required in the second method. Thus the fixed node procedure needs less com— puter time than the method of modifying elements. (c) Another advantage of using the fixed grid technique is that the position of the nodal point is fixed and the shape of each element does not change with time. Thus it is possible to calculate the location of the phreatic surface in any anisotropic and heterogeneous ll8 Input data [Build up global matricesl [Insert Dirichlet boundary condition] Solve for ¢ (with Dirichlet boundary on upper nodes along phreatic surface) [Print piezometric head of the systeml i ICalculate velocity vectors at the nodes] Shift phreatic surface along its normal and find new coordinates of the free surface [Frint the location of the phreatic surface] Check if system has reached steady state or required maximum time, stop Based on new location of free surface, modify elements such that phreatic surface coincides with upper element boundaries Figure 6—ll.——Scheme for Solving the Piezometric Heads in Unconfined Aquifers With the Movable Node Technique. 119 Input data [Build up global matrices] [Insert Dirichlet boundary condition[ above phreatic surface) Solve for ¢ (with Dirichlet boundary] [Print piezometric heads of the system Calculate velocity vectors at the phreatic surface Shift phreatic surface along its normal and find new coordinates of the free surface [Print the location of the phreatic surface] Check if system has reached steady state or required maximum time, stop Based on new location of the free surface, assign piezometric head values at the nodes above phreatic surface Figure 6-12.——Scheme for Solving the Piezometric Heads in Unconfined Aquifers With the Fixed Node Technique. 120 porous media. But in the modifiable element method, nonhomogeneity can be handled only for fixed node sec— tions or between the columns of nodal lines. (d) In the first method the boundary condition is evaluated once to yield the location of the phreatic surface. In the second method the Dirichlet boundary condition at the nodes above the free surface is computed baSed on the position of the phreatic surface. As dis- cussed in Section 6.3.3.1 there is no concrete method to evaluate the piezometric heads above that surface. Among the equations which are presented, Equation (6.3.26) gives reasonable results. Still, interpolation of the values of V3 and K and assigning the value of piezo- 33’ metric heads for nodes above the phreatic nodes, need to be investigated more. 6.3.8 Summary The calculation of Unapiezometric headstlconfined and unconfined aquifers is discussed. A modified tech- nique of calculating the phreatic surface with fixed nodes is presented, and it is shown that the velocity vectors can be calculated simultaneously to provide continuous functions at the nodes without using large computer memory. The feasibility of modeling the dispersion phenomena in an unconfined aquifer with a transient phreatic surface, which has become of universal interest, is apparent. CHAPTER VII SOLUTION OF SYSTEM EQUATIONS FOR CONVECTIVE-DISPERSION PHENOMENA 7.1 Introduction In this chapter the procedure of describing the convective-dispersion phenomena is discussed. Related mathematical equations are described in Section 3.5 and the finite element formulation is given in Section 4.3. For most practical purposes at relatively low concentra— tions, it can be assumed that the concentration of a tracer does not affect the velocity distribution [Bear 1972]. Hence the solution of a dispersion problem is made up of two independent subproblems. First, the velocity distribution is determined for all points of the flow domain. Second, the resulting velocity distri- bution is inserted into the dispersion equation, which in turn is solved to yield the concentration distribution in the flow domain. Velocity vectors play a dominant role in the convective-dispersion equation. They appear in the convective term and dispersion coefficients, and hence realistic and accurate evaluation of the velocity components is very important. Proposed simultaneous solution of the velocity vectors from piezometric heads 121 122 (Section 4.4.2) provides continuous velocities at each node. The concept is also presented by Segol et a1. [1975], where the pressure and velocity components are computed simultaneously. For an almost homogeneous liquid, the proposed method (in Section 4.4) has some advantages over their procedure: (1) the required storage space in the computer program is significantly reduced; (2) the transient solution of the velocity vec- tors is feasible; and (3) the boundary conditions are simple to apply. 7.2 Calculation<1fDispersion Coefficients The hydrodynamic dispersion coefficient defined in Section 3.5.2 regulates the degree of Spreading of the contact zone between two miscible fluids. For most practical purposes the molecular diffusion coefficient is negligible compared with the mechanical dispersion coefficient. The dispersion coefficient is given by Equation (3.5.1). The longitudinal and transversal dis- persivities are the two components of the mechanical dispersion coefficient, and are considered as porous media properties. In this study they are assumed to be constant over the entire domain. The reported values of a and a both range between 4 and 135 meters I II [Pinder 1973, Robertson 1974]. 123 Dispersion coefficients can be either evaluated at a node or assumed to be constant over each element at a given instant. For a transient flow where velocity components are also time variant, it is necessary to compute the dispersion coefficient at each node for every time step and reconstruct the global matrices (see Section 4.3). Consequently, the time dependency of the dispersion coefficient will increase the computation time considerably. 7.3 Computation of Tracer Concentration The finite element formulation of the convective— dispersion equation leads to a system of ordinary differ— ential equations of the form [HM—3%} + ([K]+[S]){C} + [E]T[I]{C} = {Huge} + {F} (4.3.13) The two most common boundary conditions, namely, the Dirichlet and the Neumann boundary conditions, are used in solving tracer concentration. The Dirichlet boundary condition is handled by specifying the known con— centration at the appropriate boundary nodes and factoring out rows and columns in the coefficient matrix associated with those nodes as discussed in Section 6.1.2. The known mass flux along the boundary line is incorporated in the {F} vector. The procedure of the allocation of a constant line source to nodal points is given in Appendix I. Since it is assumed that at the boundaries either the concentration or the flux of the tracer is defined, the matrix [E] in this study is assumed to be zero. In order to construct Equation (4.3.13), three parameters (i.e., dispersion coefficient, velocity vec- tors, and velocity gradients) should be known. Velocity vectors are calculated at the nodes, and the other two ‘ parameters are evaluated at the integration point. All 1 three are then introduced into Equation (4.3.12). The dispersion equation for flow in a confined or an unconfined aquifer is similar. With the fixed node technique described in Section 6.3.3, the piezo— metric heads in a phreatic aquifer can be obtained without repositioning the nodal coordinates of the elements; thus the velocity vectors are known for a given point at any time step. The movement of a tracer is predicted by introducing the calculated velocity vectors in the convective—dispersion equation. The nodal coordinates of the grid system for the tracer will coincide with the nodal points of flow. Thus, either the grid system for flow and dispersion will be identical, or the grid system for predicting the tracer concentration will lie within the grid system of flow prediction. The procedure for calculating the tracer concentration in a confined or unconfined aquifer is given in Figure 7—1. 125 .meesv< pwcflmsooco Ho pmcflmcoo m GH GOHHMHpcwocoo uwomue map mcfluma Isoamo How GOflpMpcmmmHmwm oauwfiwaomll.alh wusmflm :owmummmfip How wEflu Esfiflxmz Bx¢29 Ewumxm m oucfl pavzwouucfl ms quHOm umsu mass Bm¢eme meoeocme cam .mosoams . so u son cons cosuSHOm coflmhwmmflp now down wEHH Boo GesuSHOm 30Hm wow moum wEHB HQ :0flmnwmmflp now wEHu cmmeHm EB 30am now wEflu pwmmmam B unmno 30am :fl muwquMHmm new mcowuflcflmwa mmuv.mwsé A '3... —:0Husnauumac coHumuucmocou mamasoamwg , GOHmhmwnmmHmv w>HUU®>COU HOW mwufluumE Hmnon pmumHou m5 paflzm _Bpm|u as u QB_ Tcmaoawwwoo coflmuwmmflp oflEmcwcoupxc wusmfiou mm% .2 .6? mmpo: mnu um muouuo> xufloon> oumHDUHmo .wpmon uauumEonHm wumHDUHmo o H DB EEG 9Q + B H E _30Hw HOM mwofluumE cwumaou Q: paazm fimcmfiocwcm conummmflplw>fluom>coo can 30Hw Mow sump usmcH 126 7.4 Discretization of Time Derivatives and Solution of System Equatibn Equation (4.3.13) is a set of first order ordinary differential equations. Different methods to discretize the time derivative of the family of equations similar to Equation (4.3.13) were discussed in detail in Chapter V. Using Equation (5.1.5), Equation (4.3.13) takes the following form: [B]{C(t+At)} = {G} (7.4.1) where [B] is a nonsymmetric banded matrix, C is the unknown concentration for time t + At, and {G} is a known force vector. Since [B] is nonsymmetric, the Cholesky method (see Section 6.1.1) is not applicable. The Gauss elimination technique [Carnahan 1969, pp. 270—272] is used to solve the system of simultaneous equations. The subroutine GELB [IBM Application Program, 1968] has been adOpted for solving the system of equations for the dis- persion equation. Initially, a two—point backward finite-difference scheme is used, where the initial con- dition of C at time t is given and C(t + At) is sought. As the computation progresses, solutions of C at a few prior time steps become known and they may be stored in the core. Then a higher order time approximation of Equation (4.3.13) is used as derived in Chapter V, to give a more accurate solution or to permit a larger time 127 step without additional penalty of excessive computer time. As can be seen from Equation (5.1.5) the [A] and [H] matrices are fixed, and only their coefficients will change if the order of approximation or the time step changes. The additional required memory for each increased order of approximation will be 2 x NNDS, where NNDS is the number of nodes. The results of using dif- ferent methods of discretization of the time step are discussed in Chapter IX. 7.5 Stability and Convergence Criteria There are two important concepts closely associ- ated with the convergence of a particular numerical pro- cedure, namely, those of consistency and stability. Carnahan et a1. [1969] define the stability and consist— ency as the following: "In general, a solution is said to be unstable if errors introduced at some stage in the calculation . . . are propagated without bound throughout subsequent calculations." The term consistency means that the numerical procedure may in fact approximate the solution of the partial differential equation under con- sideration and not the solution of some other equations. There is no definite rule for defining stability criteria by the finite element method. For the finite difference explicit scheme, Fried and Combarnous [1971] give a 128 relation for the stability of longitudinal dispersion which can be written as follows: #2 DL At_g Axl s 2 DL/Vl (7.5.2) The stability function for the two—dimensional dispersion equation solved by the explicit finite difference method was obtained by Reddell and Sunada [1970], and is sum— marized as Dll ; D22 > 0 (7.5.3a) 4 D D > (D + D )2 (7 5 3b) 11 22 12 21 ' ' D At D At 11 + 22 5 l (7.5.3c) 2 2 2 (AXl) (AXZ) Although Equation (7.5.3) was derived for finite differ- ences, it reveals the limitation of choosing Ax and At. 1 As discussed in Chapter VIII, for 6 3 1/2 (Chapter V) the result of exceeding the time step restriction is a stable but oscillatory solution, and as 6 approaches 2/3, the oscillation decreases substantially. However, giving any specific 6 for stability criteria will be premature. Using the higher order time derivative approxima— tion of Equation (4.3.13) will improve the convergence of the finite element solution. CHAPTER VIII NUMERICAL RESULTS FOR SIMULATION OF GROUNDWATER FLOW In this chapter the validity of the finite element numerical simulation of the flow in confined and uncon- fined aquifers with a phreatic surface is discussed. The mechanisms of locating the free surface with fixed nodes and with modifiable elements are verified. Also, the results of the two methods for velocity calculation, i.e., direct and simultaneous procedures, are presented. 8.1 Flow in a Confined Aquifer The finite element solution of the flow equation for a confined aquifer has been accomplished successfully many times and will not be the subject of detailed review in this study. Since success of the convective-dispersion model is highly dependent on the hydrologic simulation model, one has to be certain that the hydrodynamic simu- lator provides accurate results. 8.1.1 Pumping in a Single Well Field Consider a single well pumping from a confined aquifer of nearly infinite area (Figure 8—1). The 129 130 ._osma couamz HmnwHMUCDOLflwsofl>HmmEH \N \\\\\\\~J»fl\\\\\\\\\\\\-\\\\\\~\m-«\~\--~\\\--h~ _: H ___ HmMstm 8 I III, 8 _ _.|| 303 __ I'. Hafipmm \ \\\\\\\\~\..~\\\\\-\\~\\\\\\\\~\~\. qqqquqqqququqqqu humpcson msofl>ummEH . . .1. GOHmmmummp I mo wcou NW mommudm pcsouu \ momwusm HHwS soapm>uwmno casmeonflm HHwB COHuOSpOHm HafiuflcH 131 analytical solution for transient radial flow in an isotropic nonleaky artesian aquifer with fully penetrat- ing well and a constant-discharge condition is: P S = W ['Ei ('u)] (8.1.1) = r2_S u 4Tt where s is drawdown, r is the distance from pumped well to the observation point, P is discharge, t is time, T is transmissibility, S is the storage coefficient, and -Ei(—u) is the exponential integral [Davis and DeWiest 1966]. Due to symmetry of the flow field, it is suffi- cient to model only one—fourth of the aquifer, shown in Figure 8-2. No-flow boundaries are taken 4,800 meters from the well located at point W, and initially the drawdown s is zero. The following parameters are used in the modeling: T = 929 mz/day, s = 0.01, and p = 946 m3/day. Figure 8-3 shows dimensionless drawdown versus dimensionless time at locations A and B (Figure 8—2), at distances of 1,697 and 2,163 meters from the well, respectively. Numerical results compare favorably with the analytical solution, and the deviation of numerical results from the theoretical curve is within acceptable range. It has been observed that the size 132 JL 1~——-——— 4800 m Figure 8—2.——Grid System Which Is Used in Simu- lating a Single Well Field. .Hawz one Eosm >w3< mumpmz mwam pom hmma .csopzmuo How chHDSHom Havau>amc< paw HMUHHmEDZ wcu mo cOmHHMQEou|I.muw mucmflm o+m.¢ mNH\qu .wEflH mmmacoflmcwfiflo o+m.H HIM.H m.o 133 _ _ _ 263 cap Eoum hmzm E mmam ‘ HHmB map Eonm hmzm E mama . _ "COHU.5._..Om HMOHHOESZ coH#5HOm amusuwamc< Ill 1 K0 6 l v I ’umopmsxa sseruorsuemrq d/SLMfi 134 and type of elements affect the accuracy of the numerical results. 8.1.2 Effects of Various Time Approximations on the Accuracy of the Results In order to examine the effects of the kinds of time approximation on the results of the numerical simu- lator, a small confined aquifer, 1828.8 X 1219.2 m (6000 X 4000 ft), was modeled (Figure 8-4). Initially, the system is at the steady state with zero piezometric heads. Along boundary lines AC and BD the system is maintained at zero potential, while no-flow boundary con— ditions are assumed along the AB and CD sides. A well is located at the center of the medium at point W. The following parameters were used in simulating the ground- water movement: P = 556.4 m3/day (2.0E + 4 ft3/day), S = 0.01, and T = 929.6 m2/day (10,000 ftz/day). Three different 0's, i.e., 0 = 1/2, 6 = 2/3, and 0 7/12 which is the average of 6 = 1/2 and 0 = 2/3, were used for first order time approximation (see Equa- tion 5.1.5). The related coefficients for 0 = 2/3 and 1/2 are given in Table 5-2, and the coefficients for 9: 7/12 are as follows: a11 = I2' a12 = 1" a13 = 1" a = ' 12' SD [.4 22 = .mucmfimam Hmmmumaflupmso UHHumEMHMQOmH Hmmcflq swag somasqm pmcflmcoo m mo coaumucmmmnmmm pcmfimam muflcflmlu.w1m wusmam Hx Hm mm am as as m H . ® ® ‘ . mm mm mm 5H NH s m 3 mm mm mm as ma m m l 3 am am am as an a a . 6 6 9 ® 0 ® . mm om mm om ma OH m mm “mmpoc Mo .02 $ «N "mucwfiwaw MO .02 136 The other coefficients of Equation (5.1.5) will be zero. The results of the calculated piezometric heads at point W with different 6's for At = 0.05 day are depicted in Figure 8-5. For 8 = 1/2, the piezometric heads oscillate substantially when compared with the results of 6= 2/3. Although 0 = 7/12 is the average of the two previously mentioned 6 values, the observed oscillation is less than for 0 = 1/2. All three 6's provide the same results for steady state as one might have anticipated. To see the effects of 0 more clearly, the time interval is increased 20 times to At = 1.0 day. The piezometric heads for a maximum time of 50 days for point W are shown in Figure 8—6. As can be observed, piezometric heads oscillate at an early time for 0 = 2/3 and 7/12; however, the oscillations diminish very rapidly and converge to steady-state condition. In contrast, for e = l/2, the piezometric heads oscillate about the mean but the oscillation does not vanish. From this example and similar studies that have been carried out, it is believed that for groundwater problems, 6 = 2/3 provides fewer oscillations when compared with 0 = 1/2 and 7/12. The 0 = 2/3 and 7/12 solutions con— verge, but 0 = 1/2 oscillates about the mean for con- tinued time. _ . . 137 ITII] I I r l I II] 22" j 23 - ‘ — \~ 24- - -—"‘—‘ e = 2/3 25— — ——————- 0 = 7/12 26 P ’\ — -—-“——— 6 = l 2 27 — \\ / ‘ \ 28— \/\ — 529- \ “ "830- __ (D a: 031 — \ ‘ H H 3332 P’ - E 8 ...1 \ 94 (A) as I / I 35... ‘ - 36- _ 37- q 38 — _ 39 ' ‘ 11111 1 l 11111| 0.05 0.1 0.5 1.0 Time, days Figure 8-5.--Calculated Piezometric Heads Versus Time at Point W Using Different 0 With At = 0.05 Day. .smo o.H n u< guns o #cwHwMMHQ mGHmD 3 bosom pm mEHB wSme> mpmmm OHMquonflm pmpmHDoHMUII.mIm wusmflm m>d© .wEHB om me ow mm om mm om ma OH m o n: _ _ _ . _ _ _ _ _ (N I. IINN «\H u a I am 3? n A III m. e I 1 m mxm u o ..... . om m i a 1 1 I IS m. > w I > . I I II II II I I I I I I . \ zximam J l o I .IsH m coausaom wumumlhpmwum I .fimfi _ _ _ fl - .. L - 139 8.2 Comparison of Two Different Methods of Solution of the Velocity Vectors As discussed in Section (6.2), the velocity vec- tors can be calculated either simultaneously or by the direct method. To show the discontinuity of velocity vectors at element interfaces or nodes calculated by the direct method, a regional aquifer is taken and shown in Figure 8—4. The area and applied boundary conditions are shown in the figure, with a well located at point W. The system has reached a steady-state drawdown condition due to pumping and the piezometric heads are known. Equation (4.4.3) is used to calculate the velocity vec- tors at the nodes and at the centroid of the elements for the direct method, while Equations (4.4.11) and (4.4.13) provide simultaneous calculation of the velocity vectors at the nodes. Figures 8-7 and 8—8 show the velocity vectors at the nodes, computed by direct and simultaneous methods for the aquifer depicted in Figure 8—4. The numbers in each box represent the velocity vectors calculated simul— taneously, while the ones in the corner of the elements are obtained by the direct method. It is obvious that direct calculation does not provide continuous velocities at the nodes. For example, if one takes element number 11 and evaluates V at node 18, a value of 7.6 will be 1 obtained. But V1 for the same node from element number 15 140 2.8 2.7 2.7 109 .81 0810.0—.81 -.81 109 2.8 2.8 7.6 7.6 -7.6 -7.6 104 2.8 2.87047.6 7. O.O_7.6 -7.6 704 2.7 2.7 .81 .81“ ~08]. -081 1 2‘8 2.7 2.7 1-9 .81 .81V°O—.81 ...81 ‘9 1.9 1.9 .80 .8‘ .80 -.80 2.3 1.5 0.0 1.5 Figure 8-7.--Computed Values of V1 at the Nodes by Direct and Simultaneous Methods. 141 “'01 -05 -6.2 ’05 -.O9 —.18 -.18 -6.9P6.9 -.18 -°09 -018 -018 “6095'609 -018 0.0 0.0 0.0 0.0 .09 .18 .18 6.9609 018 .09 .18 .18 6.9 6.9 .18 .1 .5 6.2 .5 .03 .80 .80 .81 .81 .80 .03 .80 .80 .81 .81 .80 000 095 1‘9 ’95 Figure 8-8.-—Computed Values of V2 at the Nodes by Direct and Simultaneous Methods. 142 will be -7.6. The true velocity at that node is zero, which is obtained by simultaneous calculation. Averag- ing the velocity vectors computed by the direct method might yield the values obtained by the simultaneous method, but the procedure is not valid for every location. The computation time for both methods is almost the same. The simultaneous method requires more core for storing the capacitance matrix [H] and force vectors {FX} and {Fy}. It is believed that the Galerkin formulation of the Darcy law will yield a reasonable estimate of velocity vectors (or flux) at the nodes, provided the calculated piezometric heads are good representations of the actual field conditions. An important point which has been observed in the different numerical results is that the value of the velocity vectors at the boundaries might depart slightly from the true answer. The deviation might be due to the construction of the element matrices, because some nodes along the boundary are weighted less in the formation of the global matrix and force vectors. In practice this would not cause any problem if the hydrodynamic model is used in the prediction of the tracer concentration. The grid system for the dispersion model might be placed within the hydrodynamic model, or taken slightly smaller in order to minimize the effect of this discrepancy. 143 In this study, simultaneous velocity vectors are used in the convective-dispersion program to calcu— late convective terms and dispersion coefficients. Simultaneous and direct calculated velocity vectors are used for locating the phreatic surface. 8.3 Flow in a Phreatic Aquifer The mechanisms of locating the phreatic surface with movable and fixed node techniques are presented in detail in Chapter VI. In this section the existing experimental and field observation data is used to verify the validity of the two methods. In all numerical examples illustrated hereafter, linear isoparametric quadrilateral elements are employed, and DTMAX stands for the maximum specified time step. 8.3.1 Transient Buildup of a Mound Due to Accretion Using the linearized technique, Marino [1967], following Hantush [1963], gives the analytical solution for growth and decay of a mound due to recharge. He verified the analytical solution with a Hele-Shaw model. A strip of finite height and infinite length was chosen for numerical studies as shown in Figure 8-9. Due to symmetry, only half of the system was modeled with linear isoparametric quadrilateral elements. The following data is employed in the computer program: NELS = 54, 144 (a) I Recharge = 5.6 E-2 cm/sec I‘IIIIII /,Initial Watertable L3=j 36 cm IF—g— ——— _-— 3‘84H cm X L]. -' 140 cm H 3 (b) X1 Figure 8-9.—-Simulation of Transient Buildup of a Mound Due to Accretion. (a) Grid system repre- sentation; (b) scheme of the vertical cross section. The phreatic surface is shifted along the vertical lines. 145 NNDS = 70, Kll = K22 = 25.2 cm/min, ne = 1.0, and I = 3.36 cm/min. At the beginning, a small time step equal to 0.01 minute was chosen and increased to DTMAX e 0.2 minute within several time intervals. Both the movable technique (MNT) and the fixed node technique (FNT) were analyzed, and the location of the phreatic surface at different times is shown in Figure 8-10. The results obtained by FNT and MNT compare favorably with existing analytical and experimental results. At early times, both methods provide almost the same results. As time increases, they differ to some extent, but it can be seen that the results are within a reasonable range. In these runs the velocity vectors are computed simul- taneously. The effects due to the use of simultaneous and direct calculated velocity vectors are also examined in this study. Velocity vectors are calculated and averaged at each node using the results obtained from two adjacent elements. For example, consider Figure 8-11 which shows a portion of a grid system which might be used to locate the phreatic surface. For node 6, (V1)6 = 0.5 (Vi + V?) where the superscripts represent the element numbers. Similarly, (V2)6 = 0.5 (V3 + V3). The velocity vectors computed in this manner are introduced to FNT. The cal- culated location of the phreatic surface is also compared with the observed curve as shown in Figure 8-12. The 146 I I I I I 25 ..l E u \\ - ; \.‘ c- '9 “ d 520 \ \ 2 " \ 6\ I=9mIn. '1 s - ~ -* ‘\ II- )- \ _ .\ 'I o _ \\I-5num _ " AN \ filS“ \ \ '- £ . \ .‘ \ _ b ‘~.; t=hnm. ‘ ‘~ - d \ in II— - - ‘ q - 1 O A '0 1 l 1 L J 20 4O 60 80 IOO x‘m an. Experimental . O Movable node technique ---—- Analytical }[M°“n°"967] AFixed node technique Figure 8-lO.-—Prediction of Transient Buildup of a Mound Due to Accretion With Different Techniques. The velocity vectors are calculated by the simultaneous method. G) <9 G <9 Figure 8-ll.—-Portion of a Grid System. 147 .anpmfi uomuflp men an webmasoamo mum mHOpom> hpfloon> one .coflumuooé Oh man peso: 8 mo mapaflom ucmflmcmue mo GOHpOHpmsmII.mHIm whomflm . .ponuwfi uomuflp one he pmumasoamo mum mnouom> muflooam> “Bzm NH Hecapmamcm IIIII P83 .8336 Hopcmeflummxm .Eo .Hx ooa om om ow cm 0 O H Ln r—l 'mo ’queL 1819M go qureH O N 148 results follow the expected curve and are slightly higher than those obtained using the simultaneous velocity solution. Depending on the technique chosen for calculating the velocity vectors and locating the phreatic surface, one might arrive at slightly different results. However, employed techniques provide reasonable solutions which are within an acceptable range with FNT giving the best results for this case of study. It is also shown in the literature that MNT has the capacity to simulate flow field problems which involve a phreatic surface [France et al. 1971, Desai 1972]. It can be observed that FNT is also capable of solving similar problems, if not with greater accuracy, at least with the same degree of accuracy. 8.3.2 Numerical Modeling of a Field Problem Recently, Bianchi and Haskell [1975] presented a series of field observation data on the shape of groundwater mounds produced by artificial recharge water spreading. Detailed descriptions of the experimental recharge ponds and the location of the observation wells are presented by the authors. In one of the experiments a square pond, 90 x 90 m (295 x 295 ft), was chosen and the rise of the mound due to recharge was measured. The initial position of the water table was determined from 149 well observations prior to flooding, and corrected for barometric fluctuation. The recharge rate (I) was taken to be the rate that water enters the surface of the pond, and was assumed constant. The aquifer perme— ability and the "fillable void" were evaluated based on a single pumping test at the pond. In the remainder of this section effective porosity or specific yield is used instead of the term fillable void. The water table was 7.62 m (25 ft), and the impeding layer of lower perme— ability was observed between 5.18 and 5.49 m (17-18 ft) below the ground surface. Although information regarding the numerical simulation of the system is notsufficient,it is modeled to show the capability of the proposed FNT. A vertical cross section of the site is chosen. The constant head boundary is assumed to be 244 m away from the center of the pond as shown in Figure 8-13a. Smaller elements are used in the recharge zone and its vicinity. Initially, the water table is horizontal at 22.9 m from the datum. A small time step equaltx)0.01 day is chosen and increased gradually until it reaches 0.2 day, and is kept constant until 25 days have elapsed. Taking this cross section implies that the pond is rectangular and its longitudinal length is long enough to make the two—dimensional assump- tion valid. However, this assumption does not really represent the actual field situation. 150 . Ce ter line Lr/2p Ground surface Recharge "Illlll Initial water table 57 ml _—— ———— NELS = 76 K11 = K22 = 31.7 m/day l 22.9 m NNDS = 100 ne = 0.2 x3 451n_+4 244 m (a) X1 feet O 100 200 300 U) a I i I If I I _ 8 ii 5 2.0 ~ ..6 'U . g 1.5 r» 0 M44.) 53 83 m 1.0 a, C) M o '\@\ O u 5_ ...2 ‘9 . I l l l l l l I l l g 0 o O 10 20 3O 4O 50 60 7O 80 90 100 Distance from the Center of Pond, meters (b) -——' Computed @ Observed Figure 8-13.--(a) Sketch of the Numerical Model for a Field Problem, and (b) Location of Water Table for Different Times. 151 The numerical results obtained for this model after 5.15 and 25 days are shown in Figure 8-13b. The results compare favorably with observed data in the pond area, but deviate further away from the pond. This devi- ation might have been caused by a number of factors. Some of the apparent reasons for the discrepancy between the numerical results and field observation data are as follows: 1. Aquifer parameters: Only a single measure- ment for hydraulic conductivity and specific yield is available. The specific yield is dependent upon moisture content, degree of saturation, and temperature; thus its value will differ in and outside the pond with respect to position and time. This subject is discussed in detail by Bear [1972] and has also been recognized by Bianchi and Haskell [1975]. The value of hydraulic conductivity has to be known in the recharging zone as well as outside the pond in order to obtain a reasonable estimate for the rise of the water table. The hydraulic conductivity appears in Equations (6.3.8) and (6.3.26), where it par- ticipates in locating the phreatic surface and defining the new Dirichlet boundary condition for piezometric heads, respectively. The value of hydraulic conductivity is highly dependent upon moisture content. It follows that in order to predict the location of the phreatic surface it is logical to replace K.. by K..K , where K 13 1] r r 152 (0'< K 5,1) is the relative hydraulic conductivity — r [Neuman 1973] and its variation with respect to moisture content for a typical situation is given in Figure 8-14. 2. Observed data: At the time that was taken to be t = 0, the water table was not horizontal and a rise of 0.04 to 0.39 m at the wells was recorded. The field data shown in Figure 8-13 are the average values of head rise of four wells which are located an equal distance from the center of the pond. The recorded values for each distance have some fluctuation. 3. Other factors: The boundary conditions, absence of instantaneous uniform recharge, and size of the elements might also be considered as additional r Relative Conductivity, K .1 .2 .3 Sand 0 .2 .4 I .6 Clay Moisture Content Figure 8-14.-—Variation of the Relative Conductivity With Moisture Content for Two Soils [After Neuman 1973]. 153 reasons for the discrepancy of numerical results from observed data. In summary, although some deviation exists between observed field data and simulated results in the region away from the recharge pond, the rise of the water table beneath the pond compares favorably with recorded data. It is believed that the fixed node technique is capable of predicting the rise of the free surface due to accretion in unconfined aquifers, provided sufficient information concerning the field parameters is available. 8.3.3 Steady-State Solution In Section 6.3.4 it was stated that the system will reach the steady-state condition when the velocity normal to the phreatic surface approaches zero. To show this condition, a vertical cross section of an unconfined aquifer is taken (Figure 8-15). The height of the phre- atic surface was kept constant at 126 m (413 ft) away from the center of the recharge zone in order to obtain a rapid steady-state condition. The computer run started with a time step equal to 0.01 day and gradually increased to DTMAX = 0.2 day. When the average shifting distance a (Section 6.3.4) was less than 5.0E-6 m the steady-state condition was assumed to be satisfied. This condition was reached after 19 days. The location of the phreatic surface at this state is shown in Figure 8-l6a. 154 enter line 0.1 m/day [ [ [ [ [ [ [ [ /Phreatic surface '7 22.9 m ll ‘l‘lllll u: <3 ‘*——U“-—.- B (7'17 IriTVVIv" (7' (7’7 I777Trvrlrvrv vvvrr V""' ’7’ "" 126 m i“ X3 NNDS = 55 n8 = 0.2 X1 Figure 8-15.-—Representative Sketch of a Verti- cal Cross Section of a Phreatic Aquifer and Grid System. 155 r. r. 23.5, m 3 I’ u r %’ F i‘ I: “1 22,5L_.___J__.___L____1____..I-____L____ 1-..- 0 20 40 60 80 100 120 Distance, meters (a) U) 23.5— H _ m ,p .. m E F 4323.0 5 m 4 W4 _ g 22.6 I I J I L I I I I l 0 2 4 6 8 10 12 l4 16 18 20 Time, days (b) Figure 8—16.--(a) Location of the Phreatic Surface at the Steady—State Condition, and (b) Rise of Water Table With Respect to Time at the Center Line. 156 The rise of the water table at the center line with respect to time is shown in Figure 8—16b. The free surface at first rises very rapidly, and then gradually slows while approaching the steady state. In Figure 8—17 the equipotentialljxmusand calculated velocity vectors at the nodes are shown. One observes from this figure that beneath the recharge zone and in its immediate vicinity the vertical velocities play an important role; but far from the accretion area the velocity components are almost horizontal, and thus the Dupuit approximations can be used in this region. 8.3.4 Two-Layer Aquifer To investigate the effects of two-layer aquifers on the flow regime and location of the water table, an aquifer with the following characteristics is chosen and its vertical cross section is depicted in Figure 8-18. Initially, the water table is horizontal with a height of 20 m and located 10 m below the ground surface. The size of each element is 20 m along the x1- and 5 m along the x3-direction. At zero time the recharge with 0.1 m/ day intensity is introduced, along the first 60 m of the center line of the system. A constant piezometric head is assumed at the other end. Four different cases are examined as shown in Table 8-1. At the beginning of each run a small time 157 7.10 7.0 .7 6.50 .2 6.00 .7 5.50 .2 l_____.._..l.__ L___ ___ __J__L_J__L___LA_ [‘ 122 m‘ (a) __JL_ \ ,//Location of phreatic surface &____:_____ __:flht::!:::iPL“;::lE;_::1!L::5EL_:=iI=J (b) Figure 8-17.--(a) Equipotential Lines, and (b) Specific Discharge Vectors Beneath and in the Vicinity of the Recharge Zone at the Steady-State Condition. Velocity compo- nents are calculated at the nodes and the length of the arrows represents the approximate magnitude of the vectors. 158 Center line .Accretion HII 1 Zone 1 ‘ .5} s 1 : o g T; m 8 ~ Zone 2 o 1 cu m , l .. . + I "I’I’rl’l'lr’l"fiTr'If'flllfl'IVIVI’llIIIfTWTTIVVITIVVIIIII’fiT’Trrriif’rTT 60 X3 m _‘1 400 m ’1 NELS = 120 NNDS = 147 x1 Figure 8-18.-—Vertical Cross Section of a Two-Layer Phreatic Aquifer. 159 TABLE 8-l.--Different Combinations of Hydraulic Conduc- tivities for Aquifer Shown in Figure 8-18 and Maximum Specified Time Step. Hydraulic Conductivity (m/day) Case No. ?§:§f Zone 1 Zone 2 l 30 30 0.2 2 75 75 0.1 3 75 30 0.0625 4 30 75 0.1 step equal to 0.01 day was chosen. Depending on the kinds of materials, the time step was gradually increased to a maximum specified time step as presented in Table 8-1. In all of these cases the effective porosity is assumed constant and equal to 0.3. The rise of the phre- atic surface corresponding to these conditions after 5.10 days is shown in Figure 8-19. From this figure, it can be observed that with the same effective porosity the rise of the water table at the recharge site and in its immediate vicinity is greater for less permeable aquifers. At early times, in the area away from the recharge zone the water table rises faster in the highly conductive than in the less conductive soil. As time increases this process reverses (not shown in Figure 8-19). At any time the phreatic surface remains higher for smaller K at the recharge zone. The equipotential 160 .msmo H.m umumm oHHMOMm magma young one so Homesofi Homquoze mo muowmwMII.mHIm musmflm mumpmfi .mcfiq Hoodoo EOHM occupmflo oov com com OCH 0 fl _ _ q v mmmu ....... m mmmo .l.l.l Io. HIm manna No A H N mmmo IIII sxeqem ’queL 1919M go esru H mmmo IIIIIII JN.H Akw _y' » 161 lines and discharge vectors for a two-layer aquifer (Case No. 4) are shown in Figure 8-20. As expected, the discharge vectors are greater in the zone possessing high hydraulic conductivity. 8.3.5 Change of the Effec- tive Porosity To show the influence of the effective porosity on the rise of the water table due to the recharge, an aquifer similar to the one described in Section 8.3.4 is chosen. Case No. 4 of Table 8-1 represents the hydraulic conductivities employed in this section. Three different effective porosities, namely, ne = 0.20, 0.25, and 0.3, were examined. The rise of the water table for these effective porosities is shown in Figure 8-21. With simi- lar conditions, the height of the water table is greater for the aquifer with a small effective porosity. The obvious reason for this occurrence can be seen by exam- ining Equation (6.3.5). 8.3.6 Effect of the Depth of the Water Table The depth of the water table also has a major influence on the location of the phreatic surface. In order to observe this effect an aquifer similar to the one described in Section 8.3.4 is selected, except that the height of the initial phreatic surface is lowered to 10 m. The rise of the water table after 10 days for two 162 .mmoo 0H. m Houmo czonm ohm moSHm> one .HomHson mowed I038 o How muouoo> omuonomHa OHHHoomm any new .mocHH HMHpcouomHDom AvaI. omI m ousmHm any ..II III I! IIIIIIIAI III III III III IAHIIIIIIAHIIIHII .‘IIIIAIII {IIJ I A. I A: I II I. I 0.. I I I I I I I ‘II 4. 4. 4.— z _ 0 l. o. 0.. b. 0. Al I I 0.. 8... 0.... ‘II 0.... 3.: IT .‘I 0! .f D m _ Q 0. O. O. O. In 6.. DI 0. II 0. O..- I 3.. 5| 0| «I 4! g L _ \\\ A oommHSm 0Humonnm . _ -___fr_ IIIHIIIAIIIIIIIILIIlL lllll lllIIILIIHIll II_ mNo. mo. H. mlfi. N. m. w. m. m. K.. 163 .moHuHmOHom o>Huoommm ucouommHo Mom whoa H.m Hopmn ooowusm OHHMoHnm on» mo coHuoooqII.HmIm oHomHm muopoE .ocHH Honnou Scum oocoumHQ ooa com com OCH 0 _ _ _ H I o H T. S 8 m. m. M E a. a I L m. 1% TL 3 see}. me m w . 1. 8 mega om u was I is I 0.1“ 164 different initial depths is given in Figure 8-22. As one might have expected, the mound in the recharge zone and its vicinity will be more pronounced for a shallow aquifer than for a deep one. 8.3.7 Decay of a Groundwater Mound The investigation of the decay of a mound is an interesting subject in groundwater hydrology. To show the capability of the proposed technique of handling such cases, an aquifer similar to Figure 8-18 is taken and recharge is applied up to 4.9 days before being halted. The rise and fall of the water table and the location of the phreatic surface at the center line, up to 9.9 days, are shown in Figures 8—23 and 8-24, respec- tively. In this example the hysteresis and variation of the specific yield with time are neglected. It can be observed that more time is required for the water table to decay to the initial steady-state condition than the time necessary to build the mound. 8.3.8 Maximum Applicable Time Interval Specifying the time interval is an important task in solving transient phreatic problems. There is no par- ticular rule for selecting the time step. A general procedure which is useful for defining At is given in Section 6.3.5. In all the different examples illustrated 165 .msowcmmoEon mfl Bantu: .msmo o.oa “mama mommusm mmum map mo msumma amHuHcH ucmnmmmao you magma “mums mo mmflm--.mmnw masons mumums .mcflq Hmucmo Eoum mocmumflo oov oom .oom OCH 0 _ _ _ a was m.o u x425o «mm Ham %MU\E cm H sxeqem 'atqem 1819M go 6513 mumumE om I. la! mumume 0H nzpmmc magma umumz Hafiuflsfl Sufi: 166 .msomammoEon ma Esflpmz .HmMst< oapmwmsm m ca canoe Houmz mo Hash pcm mmflmul.mmlm musmflm mumpmfi .oCHA Hmucmu Eouw monopmflo oov com com OOH mmc\s om n mmm mmmv m.m.:ll.2ll.fll mhmw m.mnl mmmw m.q mNMU o.m.2..1..2.3. saw o.H.|u ..... mac o.o ......... Hmumm momuusm.owummnzm on“ mo coflumooq v.0N w.om m.om o.HN Slalom ’SIQQL 1319M go qureH 167 .mcoN wmumnomm m mo mafia Hmucmu mnu pm wanna Hmumz mo Hawm cam mmflmll.wmlm mnsmflm mhmw .mEHB OH m m n m m w m m H o _ _ _ _ _ _ _ _ ~ _ smw\s 0.0 smw\s H.o u H n H sxeqem ’queL 1849M 30 iqbreH — o — o — p — o _ u — O — c _ u — o Fri, . . 168 in Section 8.4 a small time step was chosen, and using Equation (6.3.28) it was increased to a maximum specified time step and maintained constant thereafter. Defining the maximum time step requires some experience which might be obtained after a few runs for a specific problem. In general, for highly conductive porous media or nonhomof geneous aquifers, a smaller maximum time interval should be employed. Usually, the rise of the water table fol- lows a trend similar to that depicted in Figure 8-l6b. When it starts to oscillate with respect to time, the maximum time step should be decreased. For example, consider a phreatic aquifer repre- sented by Figure 8-15, in which Kll = K22 = 31.7 m/day. In Section 8.3.3 it is stated that 0.2 day was used to obtain steady—state conditions for this specific problem. To show the sensitivity of the numerical solution to the time step, the hydraulic conductivities were increased 2.5 times. For the first time, DTMAX equal to 0.2 day was used to find the location of the phreatic surface at the steady-state condition. As expected, this value was too large and the results were not correct. Then DTMAX was reduced to 0.1 day. Although the numerical results were realistic for early times, after 4 days the values of the piezometric heads started to fluctuate. The oscillation increased with time as shown in Figure 8—25. By further reduction of DTMAX to 0.05 day, the oscillation 169 .mwum mEHB wmumq m on wso mHQmE kum3 m0 coflpmaaflomoll.mmlm madman .mnsmflw mflau an c30£m mfl mafia Hmpcmo map um manna HmuMB may mo mmflm mhmc .mEHB c I saw mo.o u u<_apflz 4. COHuSHOm mumuwlmpmmum:nl was mo.o saw H.o saw ~.o u< #4 HQ 4. O :--l O N 'mo 'etqe; 1919M go 35:3 0 m ow om 170 was eliminated. It can be concluded that 0.05 day is an appropriate DTMAX for this specific example. 8.4 Summary In this chapter it was shown that the finite ele- ment technique is capable of solving flow problems, both in confined and unconfined aquifers. Under transient conditions the fixed node technique (FNT) was tested and it was verified that this method can predict the location of the phreatic surface and can yield piezometric heads and velocity vectors with a reasonable degree of accuracy. It was shown that the simultaneous solution of velocity functions at the nodes is continuous and acceptable. Finally, the effects of the order of time approximation on the accuracy of the predicted results were examined. CHAPTER IX NUMERICAL RESULTS FOR PREDICTION OF CONCENTRATION OF A TRACER The computer model developed in this study solves a set of partial differential equations. One equation is the combined equation of motion and continuity of flow which describes the piezometric head distribution of the aquifer; in turn, the velocity components and hence the dispersion coefficients are computed. Then the sec- ond partial differential equation, the mass-transport equation (convective—dispersion equation), is solved to yield the concentration distribution in the flow domain. The numerical results for the flow equation were discussed in Chapter VIII. In this chapter the finite element model for simulating mass transport is verified by comparing the numerical results with several existing analytical solutions. If the results of a known analyti- cal solution can be approximated, a great deal of confi- dence in the numerical simulation can be gained. The feasibility of the model to predict the concentration of a tracer in field problems is shown. In this chapter concentration is used as a synonym for dimensionless l7l 172 concentration, i.e., C/Co, CO being the initial tracer concentration. The DL and DT are longitudinal and trans- versal dispersion coefficients, respectively, and are assumed constant for the porous medium independent of velocity. 9.1 Longitudinal Dispersion in Steady Uniform One-Dimensional Flow A semi—infinite cross section of a homogeneous and isotropic porous medium with a plane source maintained at x1 = O is shown in Figure 9-l. C(oo,t) = 0 ..".“-‘V = 11...n —‘>Xl Figure 9-l.-—Cross Section of a Homogeneous and Isotropic Porous Medium With a Plane Source Maintained at x1 = O. The flow is maintained at a constant flux ql in the x1- direction. For an isotropic medium, the axis of the dispersivity tensor is assumed to coincide with the velocity vector. Equation (3.5.1) reduces to D ——-—-V __:___ (9.1.1) 173 where DL is the longitudinal dispersion coefficient, and V is the velocity vector along the xl-direction = ql/ne. l The solution of Equation (9.1.1) with the following ini- tial and boundary conditions, C(xl, 0) - 0 x1 3 O C(°°, t) = 0 tZO is given [Bear 1972] as: x-Vt Vx x+Vt %— = % erfc —£——l— + exp erfc —l;—l—- (9.1.2) o 2/D t D 2/D t L L L where erfc(u) = l - erf( ). °° 2 erf(u) = —2— e a d a (9.1.3) /? O The prOgress of a concentration front in an infinite column of a porous medium is modeled numerically on the computer, as shown in Figure 9-2. A constant source is maintained at x = O and the following parameters are employed in the finite element model: Ax = 0.4 cm, 1 0.01 Vl = 0.1 cm/sec, L = 10 cm, At = 2 sec, and D cmz/sec. L 174 —*1 Axi r'_‘ 3C _ .' ..I WI ' ° nut L fi’l Figure 9-2.--Finite Element Model to Simulate One- Dimensional Longitudinal Dispersion. The results of runs using the above data are shown in Figure 9-3. The numerical results compare favorably with the analytical solution; however, in some instances they are slightly higher. Reddell and Sunada [1971] have solved one-dimensional longitudinal dispersion with the method of characteristics using 2 and 4 moving points per grid. Since the finite element technique used in this study does not employ moving points per grid and since the Ax in this study is greater than the one used by Reddell and Sunada, it is difficult to compare the degree of accuracy of the two methods. The order of time approximation has a great effect on the accuracy of the numerical solution. A sec- ond order approximation, introduced in Chapter V, was used to obtain Figure 9—3. The first order time approxi- mation provides poor results as shown in Figure 9-4. In Figure 9-4 for 8 = 1 (implicit method), the numerical 175 .muHSmmm HmoHumHmcd mnu nuH3 GOmHHmmEoo cam coHumEonummd mEHB Hoppo ocoomm mchD .owm cm H u can .oom cm H u um conHmmmHo HmchsuHmcoq may mo COHunHom HmoHHwEDZII.mnm wusmHm .Eo .mocmumHo 03/3 ’uorgEJguaouoa OH m m n m m v m m H o 0 [Ho .IlNo 1m. 1?. m\m .ND u o ‘ 1.... H u o .0 u.xonmmm mEHu :10. Hopno ccoomm IN.- GOHusHom 333383 I. 1 m . |.m. o.H 176 .muHsmmm HMOHuhqum on“ aqu somwummfioo cam coHmeonnmm¢ mEHB Hmpuo uwHflm mchD .omm om n u can .omm cm H u um GOHmHmmmHQ HmcHUsuHmcoq map mo coHudHom HMUHHmEUZII.¢Im ousmHm .Eo .mOGMDmHQ . 0H m m h m m w m m «\Hno I m\~uo4 ans 0 u.xoummm mEHu .umcuo umHHm cOHusHOm HMUHumqumall 03/3 ’uorgexguaouoa 177 results fall below the curve close to the point source and then move above the curve as the mass travels further from the origin; 6 = 2/3 gives good results at high con- centrations, but diverges later; 8 = 1/2 yields a higher concentration than the analytical solution at every point. The second order approximation provides a better estimate for tracer concentration compared to first order. Pinder [1974] examines the stability of the first order time approximation by Fourier analysis of the mass-transport equation. He concludes that 6 = 0 provides unstable results for any spatial distance, but 8 = 0.5 and 8 = 1.0 give stable numerical results for any reasonable space mesh. Cheng [1973] mentions using second and third order implicit approximations of the time derivative in the solution of the convective—dispersion equation. Although third order approximations were derived in Section 5.5, they were not examined in this study. 9.2 Sensitivity Analysis for Time Step and Grid Size To investigate the effects of the types of e on the convergence of the numerical solution using larger time steps, several At were chosen and the above data was employed to solve the concentration distribution for a porous medium 14 cm in length. Let the residual be the difference between the analytical and numerical solutions for a given x Then define ZT as the averaged sum of l' the squares of residuals, i.e., 178 NT 2 (c* - c) ZT = “=1 (9.2.1) NT where NT is the number of available data from time zero to 90 seconds, and C* and C are the analytical and numerical values of concentration at x1 = 2 cm, respec- tively. The variations of ZT with respect to At for different second order 6 values, namely, 1/2, 2/3, and l, are given in Figure 9-5. For small time steps up to At = 4 sec, 6 = 1 gives better results compared with the other two 6's. For greater At, 6 = 2/3 is superior and the best results for this 8 are obtained at At = 5 sec. The value of ZT for 8 = 1/2 is always the highest, and for At greater than 4 sec it increases very rapidly. This means that for 8 = 1/2 at larger At the results diverge and are not correct. To see this effect very clearly the variation of concentration using At = 7 sec at x1 = 2 cm with differ- ent 8 is given in Figure 9—6. It can be observed from this figure that poor results are obtained using 8 = 1/2. The‘probable conclusion that one can make from Figures 9-5 and 9-6 is that for larger time steps 8 = 1/2 diverges and the results oscillate, and as At increases the error increases gradually for both 8 = 2/3 and 8 = l, with 8 = 2/3 showing the least error. 10 I I I I I I I I I I I I I III L Jil ll lllll] \ llllll I 1 x. I I llllll llllll NT n: TV I \ I I .OOl At, seconds Figure 9-5.--Variation of ZT With Respect to At for Dif- ferent 8. Second order time approximation is used. 180 01.4.. L Concentration, C/C l l I l l l L L 0 20 40 60 80 100 Time, sec. 0 6 = 1/2 .A 9 = 2/3 E] 9=l Figure 9-6.--Effects of Different 8 on the Accuracy of the Concentration Distribution at x = 2 Cm. With At = 7 Sec. Second order time approximation is used. 181 To show the effect of the grid size on the accuracy of the numerical solution for the convective- dispersion equation obtained by the finite element method, one can define ZX as the average sum of the squares of residuals, i.e., NX * 2 (c - C) zx = “=1 (9.2.2) NX 2 * In Equation (9.1.5) C and C are the analytical and numerical values of the concentration, respectively, at a specific time, and NX is the number of available data between x1 = O and x1 = L. Different grid sizes from Axl = 0.2 to Axl = 2 cm are employed, and second order time approximation with 8 = 2/3 is used. The variation of ZX with respect to Axl for t = 80 sec is given in Figure 9-7. It can be observed that the error is reduced as the size of the elements becomes smaller. 9.3 Two—Dimensional Dispersion With Uniform Flow For steady and uniform flow parallel to the x-axis in an isotropic and homOgeneous porous medium, the convective-dispersion equation can be written as 2 2 39 = D é—E + D 3—9 — v 39— (9.3.1) T 2 1 a 2 X1 182 .omm om u u an x< Op pommmmm QDHB xN mo coHpmHum>II.nIm musmHm o.N m.H m.H v.H .EO N.H ~Ufl< o OH m. w. v. N. I IIIII I IIII H _ L _ 1 _ _ H - 1 11111.1 1 1111] Hooo. Z X __ MT. Hoe "u N .F71x MW _ D x. Z Ho. 183 where DL and DT are the longitudinal and transversal dis- persivities, respectively. Consider a porous medium as shown in Figure 9-8. The seepage velocity V is uniform l and constant throughout the medium. x2 , +’ Porous medium = :2 "_..'. .0 ."_~.'.'....-‘.‘._... . ... . n" I '- '1... ~ “I m. o —~>- .Vl- - ‘ - >- X1 C/CO l _gj-jfjtffj:]::. Impervious boundary Figure 9-8.—-Two—Dimensional Dispersion With One- Dimensional Flow. At x1 = O, the concentration is held constant and equal to CO for x2 < O and equal to zero for x2 > 0. The medium is confined by two impervious boundaries far apart from each other, such that the concentration remains CO and zero at the bottom and top of the boundaries, respectively, for all XI. The following boundary conditions can be stated for steady-state distribution of the tracer con- centration: fi"_—‘ ~~- '- '--~:::.;:z 184 C(O, x2) = O 0 < x2 < w EE— = 0 x = i w for all x 3x2 2 1 (9.3.2) For this Special case Equation (9.3.1) reduces to 2 v 33—: D 33—9 (9.3.3) 1 3x T 8x2 1 2 The solution of Equation (9.3.3) with the indicated boundary conditions is given by Harleman and Rumer [1963]. The solution is X g_.:%erfc 2 (9.3.4) 0 2I/DTxl/Vl In simplifying Equation (9.3.1) to Equation (9.3.3), Harleman and Rumer assumed and experimentally justified 2 2 2 2 that 3 C/Bxl << 8 C2/3x2. The movement of a tracer in two dimensions with one-dimensional uniform flow is modeled numerically, as shown in Figure 9-9a. The following parameters are 7 cmz/sec, D = employed: V = 0.1 cm/sec, D = l x 10- T l L 0.01 cm2/sec, Axl = 1.0 cm, and Ax3 = 0.5 cm. The steady—state solution of the problem by the finite element technique is obtained using the coefficients given in Table 5—3. 185 C/Co=0 " I cxco =I l.0 0.8- Q 0 c .6- 0 g C C G C C t .4)- S X2=2.0cm 7 O 5 ) o .2- 0 O I L I I I I I , 0 2 4 8 IO [2 l4 l6 l8 20 X I Distance Along X| ,cm Analytical solution C) Numerical solution (b) Figure 9-9.-—Steady State Solution of Two—Dimensional Dispersion With One—Dimensional Flow. (a) Representative sketch, (b) concentration at different x2. Cubic quadrilateral ele— ments with NELS = 14 and NNDS = 98 are used. 186 The numerical results using cubic elements are compared with the analytical solution in Figure 9-9b. The results compare quite well with the analytical solution except in the vicinity of X1 = 0, where the results are slightly lower than those obtained by the approximate analytical solution. This slight deviation from the analytical solution for small x might be due to neglecting 32C/8xi in Equation (9.3.3). The type and size of the elements may be considered other reasons for this small discrepancy. 9.4 Point Source With Uniform Flow In this example a source is maintained at a constant concentration for t > O at point A in a porous medium, as shown in Figure 9-lOa. The parameters used in this example are: DL = 0.01 cmz/sec, DT = 0.001 cmz/sec, Vl = 0.1 cm/sec, V2 = 0, At = 5 sec, Axl = 1 cm, and Ax2 = 0.5 cm. The concentration distribution of the tracer at x2 = 0.0, x2 = 0.5, and x2 = 1.0 cm after 180 seconds is depicted in Figure 9-10b. Although V2 in this example is zero and DT is smaller than D tracer move- LI ment along the x2-direction is evident. This implies that velocity components facilitate the spreading of chemical substances, but in the absence of a velocity component the tracer will move due to the dispersion phenomenon. 187 . 2.. ...,..'.......:.( Moianource . . '° . .' i . . . ' . a . . XI o ' ' . . ’ - '0, ...'..°. . ....-....3C g, «.BC__0.——I-V. “.oe-r....__=o I%' . . - _I I“ L,= 20 cm ii (a) o in O) In Concentration C/Co iv l l C) C) 4 6 Distance Along X, , cm (b) Figure 9-lO.-—Two-Dimensional Dispersion With a Point Source. (a) Sketch of the system, (b) con- centration distribution at time 180 sec. along x at different x . Linear quadri— lateral elements with NELS = 120 and NNDS = 147 are used. 188 9.5 Two-Dimensional Dispersion With Transient Flow To show the capability of the numerical model in solving the convective—dispersion problems with transient flow conditions, two examples are given in this section. In the first example a porous medium is taken with boundary conditions and employed parameters shown in Figure 9-11. Water is withdrawn at a constant rate at point W to produce the piezometric gradient. The vari— ation of the piezometric head, magnitude of the velocity vector, and concentration for point A with time is given in Figure 9—12, and the concentration distribution for the system after 100 seconds is depicted in Figure 9-13. An important observation which can be made (see Figure 9-12) is that the simultaneous solution of the velocity vectors is sensitive to the fluctuation of the piezometric heads. Even so, the concentration dis- tribution is a smooth curve. In this study the first order time approximation was used to solve the flow equa- tion, while the second order time approximation was used for dispersion. It is believed that the employment of a higher order time approximation helps in the smoothing process and it might have been used in the hydrodynamic problem. However, this idea needs to be investigated further before any concrete conclusion can be drawn. 189 .COHmHQOHQ pom >DH00H¢> oHQmHHm> oEHB QDHS Esprz msouom m mo cHoEom m>HDmpcomoumomII.HHIm onomHm .oom m bVH omd ED Hoo.o HmHoDmHHupmsv HmocHH "mpcoEoHo mo om%9 Hm u< Eo Ho.o mmzz oom\mEo H.o meZ Hoo.o n HHm omm\NEo m.o ll E—IUJO-I ii EU OH II E0 m.h fifilfio NII. .‘____mo €-——i- l a ; . s 1 . \ f .. . . x 1. .. . 2 .. I 1. . . I Q Q me 1 6.." If." 3 1. I U.® J\I o u s l \ . . 1 . . 1 . I . 1‘ . 0 o 1 \ s \ 1r. . I s s n ‘ C a o . . . n I I I I I I I . I u a u I u - o . I I u . . v . . ~ . . . . . . . . n I \ I . . a a I . u . n x . I . \ . l c I Q Q s \ I u c u I . a I I . n I I . a o o o \ . o . \ . I . o h o I . I n r I u I I n \ I . I . . u . u o o . a u a u I . u o I I . u u I I I \ n p I o n . o In\hLNLNANLPF-§IEIsLINLNIKL \NshnbihblLKLhKKLPNLNbethbkthhhkt . I - ‘ I ...I ll . . . I u I.. .0 . .I u no \ I I I. 1-100 I .p I cl. 0 In. 0 I I .w u I .o I ~. I I I II pa I u ...I .I II... n|u a ...I. - _— OI. I. ‘I... ‘ ‘ “I .Ib -\ \ I. . a I \ .I I '- ... . . . \ . .. . . u s.' .0 AU. . I u \ u I ‘ I I . . . . . .. U U. m I. Is. I uI... .\ .I I.. ‘\~u... . c C.— I . .. . I \\ I I. o I ‘I I § I I n 1 III - o I I v - - ... \ ~\ . . ... . .. .. I \ a . a. . . I . . Q . . .. . On I ‘I ‘ \ p I‘ .II -I .~ I. ‘ I ' I I. Q. . .c I I . .. ~.. . .n. . . I a \ 0 I III 0 n . . n I .. . s s . . o . . ..~ I ‘ .0. ~ - Is. '0 Q .0 ‘ \.. \ n a. I I. c - - . . . . . . “1II“‘V‘~‘INV“QWV‘N‘V‘fi‘“‘QI“I““‘I‘CNN‘N“V“I“‘VQV““Q‘VVVI‘I L Nx 190 —.l _- 0 WI . I: 5 g -012"- . . o q; Piezometric head N to .1 £3;%-.14 b m -.16 P 8.0361. m ‘H \ ° '5‘ '8 .034 - Velocity .‘3 >1 4: u '3 '3 032 tn 0' - ‘2“ '3 ;> .030 b Concentration 1.0 )- .8 - c o -H JJ .6 I- o H 2 (D .4 — o c 0 U .2 I- 0 _ I I J I I I I I I 10 20 3O 40 50 60 70 80 90 100 Time, sec. Figure 9-12.--Variation of the Piezometric Head, Velocity, and Concentration With Time at Point A. 191 .mpcooom ooH Hmpw¢ pHmHm son ucmHmcmne m SDHB EsHpmz a How mmCHH COHumuuchGOUIOmHII.MHIm ousmHm I <2 °mo g ——————j_i 192 In the second example a confined aquifer, 80 x 80 m, is chosen. The concentration and the piezometric head at point C are held constant, and a well is located at point W, of Figure 9-14. A high pumpage value is chosen to produce a steep gradient in order to obtain high velocity vectors. This in turn will cause a tracer buildup at a lesser time. Initially, the system is at zero potential, and after three days of pumping a constant source of tracer is introduced at point C. The variations of concentration at points A and B, 10.00 and 14.14 meters away from the point source, respectively, are shown in Figure 9-15. The concentration distribution at the aqui- fer after 10 and 20 days of initial pumping is shown in Figure 9-16. Since dispersivity along the xl—direction is greater than along the xz—direction, the tracer will advance more in the x1- than in the xz—direction. The curves are non—symmetric as one might have anticipated. As time increases the asymmetry becomes more pronounced. As far as it is known by the author, no published results are available to compare with the obtained numerical results. The only conclusion that can be made at present is that the simultaneous calculation of the velocity vectors and the use of higher order time approximations can be considered to be a step in predict- ing tracer concentration in transient groundwater flow problems. 193 = K22 = 30 m/day T = 900 mZ/day s = 0.005 aII = 5 m e E 3 o P = 100 m /day m NELS = 64 NNDS = 81 _LE [3 B SI A 10 m —+8O m Figure 9-14.--Domain of Two—Dimensional Dispersion With a Transient Flow Field. 194 .m can ¢ mucHom um mEHB zqu COHDMHucmocou Houmue mo GOHDMHHM>II.mHIO whomam mmmp .oEHB ON ON mm ON OH OH OH NH OH O O O N O _ H _ _ _ _ H _ _ _ _ _ _ O O u o .m N 4 1 IL I w... v I. 0 1m. u m uQHom pm mm 1%. D o ation um Io; 195 if 50 I o o c o q o o _.__—— 10 days 40_ 20 days 30 m H o ..IJ o E 20 10 -.-n—a—IO - —-—O ’“x 0 10 20 30 40 50 meters Figure 9-16.--Iso-Concentration Lines After 10 and 20 Days. 196 9.6 Dispersion in a Phreatic Aquifer With Accretion To investigate the capability of FNT (fixed node technique) in solving the convective-dispersion equation, a portion of a phreatic aquifer was taken. A sketch of a vertical section of the chosen model is depicted in Figure 9—17. Initially, the water table is maintained horizontal at 20 m, and both ends of the system are kept at a constant piezometric head. It is assumed that the recharge site is a continuous source of a tracer with constant concentration, and the tracer maintains the same concentration until it reaches the water table. The rea- son for choosing this model is to develop the conditions such that piezometric heads and velocity components vary considerably within the period of interest and trends of velocity components will differ in the system. The procedure outlined in Section 6.3 was used to locate the phreatic surface and define the velocity vectors at the nodes above the free surface. It was assumed that recharge directly reaches the phreatic sur- face. In calculating the location of the phreatic sur- face the primary interest was to find the velocity vectors on the surface, and the above assumption was perfectly applicable. In the dispersion process the primary source of the pollutants is located on the ground surface and the material is washed and carried 197 I = 0.1 m/day iiiii Tb S7 lyInitial water table 2 I x O E A (O O i ‘i 40 m ‘1? 400 m OH.” X NELS = 120 = = 3 K11 K22 30 m/day NNDS = 147 ne = 0.25 aI = 20 m aII = 5 m X1 Figure 9—17.——Representative Model Used in Simulating Tracer Movement in a Phreatic Aquifer. 198 downward by infiltration. This event suggests that in the dispersion model at the accretion zone, the value of V3 has to be equal to the recharge rate. This modifica- tion of the assumption is realistic for the dispersion phenomena and improves the numerical model. The following parameters are used in the system: NELS = 120, NNDS = 147, K11 = K22 0.25, I = 0.1 m/day, aI = 20 m, aII 0.2 days. As outlined in detail in Chapter VII, velocity = 30.0 m/day, ne = = 5 m, and DTMAX = vectors and hence dispersion coefficients are calculated and then used in the convective-dispersion equation. A constant concentration of a tracer is introduced at the end of the second day. The reason for choosing the second day is that usually at most sites prior to recharge a potential gradient exists in the porous media, and secondly, it takes some time for a tracer to travel from the ground surface to the water table beneath the recharge site. However, the starting time for computing the dispersion of a tracer can be any time, including t = 0. The discharge vectors at t = 2 and t = 8 days and the iso-concentration lines at t = 4, 8, and 12 days are depicted in Figures 9-18 and 9-19, respectively. From the consideration of physical aspects, a solute cannot move in a completely dry soil. This condition in turn implies that the concentration of a tracer at the nodes 199 . Center line I | ...-_-_-..L—uh—an-nan—uarH—Dh—fli-h-F'4*.’-*¥' { :::: |‘__nL—pn—Hr=—4-—u-——rL—4>-—au i \N\ -———.L—an’—.P——f'——i>——’-——'M——.~ ‘ ~H§.___.r—4R——.~—4P——¥D——.h—4>-—iF——.- O (b) Figure 9-18.—-Discharge Vectors in a Phreatic Aquifer With Recharge After (a) 2 and (b) 8 Days. [aria—7x Dry soil wWater table \.~. 200 Recharge = 0.1 m/day HHHHHHH/Gmundwrface 1754: 57/1,” ----0—5-—--—" 0.1 L.___ 40m ___... A‘HHHHHHH —' .M 0—0 0-.-. .- ..- -.—.-. ‘0', ____i__g ——-——__L_ __ ° \ \\\ \ _______________ .// o\ ./ ‘ ........... 0. 1 ,,,,,,,, ’ (b) | Ab Inc /% no <:;::(::i:;::<:;//' I“; no 26 i I _ ‘ \ / / 1 l I k \\. .// j l ' x / | I \I\. ‘——___‘_O__5_--__,_..—’ ’. 1 .‘0‘°~.‘._'_.__._._.-9:l.-.=‘ __ ’’’’’’’’ I l (C) Figure 9-19.—-Iso—Concentration Lines in a Phreatic Aquifer After (a) 4 Days, (b) 8 Days, and (C) 12 Days, 201 above the phreatic surface should remain constant. To investigate the capability of the proposed numerical technique for handling such conditions, two different computer runs are examined. In both tests all condi- tions are kept the same, except that in one run the con- centration is set equal to zero at the nodes above the phreatic surface, while in the second run the computer is allowed to calculate the concentration at these nodes. The numerical results for concentration distribution at the interior nodes are the same for both runs (see Figure 6-10 for definitions). But it is observed that a slight concentration of the tracer appears at the nodes above the phreatic surface for the second run. The numerical technique might account for this slight devi- ation. At each element adjacent to the recharge zone the concentration is kept equal to l at two nodes, so the numerical solution will compute a small concentra- tion at the other two nodes. Since the computing of the concentration distribution beneath the phreatic surface is the major concern, the deviation discussed above will not contribute any error. As far as it is known by the investigator, no other results exist to compare with the computed velocity vectors and tracer distribution. How- ever, the results are promising and seem realistic. 202 9.7 Reasonable Time Step for Calculating Tracer Concentration in a Phreatic Aquifer As discussed in Section 9.2, the numerical solu- tion of the mass transfer equation is sensitive to the time step. Usually, an allowable time step for the solu- tion of the convective-dispersion equation in confined aquifers is smaller than that for flow. In the process of calculating the location of the phreatic surface, choosing the prOper time step is very important and instability and oscillation occur for larger time inter- vals (see Section 8.3.8). Experience shows that in phreatic aquifer problems, the same time step can be used for the solution of both flow and mass transfer phenomena. For example, it was found that DTMAX = 0.2 days is adequate for obtaining the location of the phre- atic surface in the study of the previous section, and the same time step can be used for the solution of the mass transfer equation. Reduction of the time step to 0.1 day does not have any effect on the results, as shown in Figure 9-20. In this figure the concentration distri- bution at point A (of Figure 9-17) versus time for two different time steps is given and the results are identi- cal. It is obvious that the time step for the solution of the convective-dispersion equation should not exceed the DTMAX . 203 .mmoum wEHB uanDHMHQ mchD oEHB QDHB 4 ucHom um aoHDmuucwocoo Mo COHDMHHM>II.ONIO oHDOHm mhmp .wEHH 3 3 m m n m m s m m H o . _ _ _ . _ q n . o l l N0 I JV. .II 1m. I was N6 n 2 o I m. I was To I 2 II I 04 "SDHB muHsmwu HMUHMoESZ _ _ . _ _ _ _ _ p _ 03/3 'uorgexgueouoa 204 9.8 Summary In this chapter numerical examples are presented to show the applicability of the employed techniques in predicting tracer movement, both in confined and phreatic aquifers with a transient flow domain. The sensitivity of the numerical results to time steps and orders of time approximation was examined. Primary results reveal that the fixed node technique is capable of solving the convective—dispersion equation in an unconfined phreatic aquifer. CHAPTER X SUMMARY AND CONCLUSIONS In this study the movement of a tracer in an aquifer is investigated. The tracer may be introducedeus 1 a constituent of artificial recharge, for example, a chloride ion present in treated sewage water. It is assumed that the tracer remains unaltered in the aquifer. Both two-dimensional regional (horizontal) and two- dimensional local (vertical) flows are considered. In order to accomplish the movement of the free surface within the grid system without repositioning the nodal coordinates of the elements, a procedure for locat- ingeaphreatic boundary of an unconfined aquifer is adapted. In order to obtain continuous flow across ele- ments and at the nodes, the Galerkin formulation of the Darcy law is constructed and velocity vectors are calcu- lated simultaneously at the nodes. These transient velocities are subsequently used in shifting the phreatic surface as well as computing dispersion coefficients and convective terms of the mass-transport equation. Finite element formulation of the flow and convective-dispersion equations leads to a set of first order partial differ- ential relations. Using the finite element concept, 205 206 higher order time approximations for the system of equations are derived. The validity of the proposed techniques is established by first comparing the numeri- cal flow results with existing analytical, experimental, and field data. Upon verification of the solution of the flow equation, the prediction of the movement of a tracer in an unconfined aquifer with a transient phre- atic boundary and in a confined aquifer with a transient flow condition is conducted. Numerical examples have been presented to demonstrate the capability of the proposed techniques. It is shown that the Galerkin finite element method can be used to solve the flow and convective- dispersion equations, both in confined and phreatic aquifers under time—variable flow conditions. On the basis of the present study the following conclusions can be made: 1. The fixed node technique is capable of locating the transient phreatic aquifer due to accretion. The numerical results compare favorably with experimental and analytical solutions. The major advantage of this method is that the solution of the convective—dispersion equation for an unconfined aquifer with a movable free surface is possible. The primary results of predicting the tracer movement seem realistic. Much work yet 207 remains to be done in order to complete this investiga- tion and cover all related aspects of the problems. 2. The Galerkin formulation of the Darcy law provides continuous velocity vectors across element boundaries. Calculation of the transient velocity vec- tors based on known piezometric heads becomes straight— forward. 3. The time approximation of equations describ- ing the transient behavior of the field problems is an important factor on the stability and convergence of the numerical results. A second order time approximation gives more accurate results for convective-dispersion problems than a first order one. 4. The finite element numerical technique pro- vides facility for solving field problems related to flow and mass-transport situations. It yields accurate and realistic results provided that the physical behavior of the phenomena under investigation is well understood, and related parameters and initial and boundary condi- tions are properly specified. CHAPTER XI RECOMMENDATIONS FOR FUTURE STUDIES 1. This study is the first step in predicting the tracer movement in an unconfined aquifer with a transient phreatic surface. The fixed node technique is used to find the location of the phreatic surface. The flow solution of this technique agrees favorably with experimental and analytical results and field data, and is believed to provide reasonable results particularly when the rise of the water table is caused by accretion. The calculated tracer distributions appear to be realistic. Still there is more to be done, especially in the areas outlined below: a. Equation (6.3.26) was used to obtain the value of the piezometric head at the phreatic nodes (see Figure 6-10) with constant effective porosity and hydraulic conductivity. In Section 8.5 it was stated that one oftflmapossible ways to reduce the error is to modify the Darcy law for unsaturated flow. This means that in order to improve the technique, especially for the area away from the recharge zone, it is required to 208 . .- 209 extend the work for saturated and unsaturated porous media. b. Assigning the value of the piezometric heads at nodes above the phreatic nodes needs to be further pursued until more realistic condi- tions can be obtained. As was already discussed, this task becomes more important when the velocity vectors are calculated simultaneously and used in shifting the phreatic nodes, calculating the dis- persion coefficients, and evaluating the convec- tive term of the dispersion equation. c. Extending this work to three dimensions will be a significant contribution and will reduce the errors that are associated with two— dimensional assumptions. d. Obtaining field and laboratory data regarding the movement of dissolved chemical sub- stances in unconfined aquifers will help to verify the numerical models, so that increased confidence will be gained in more sophisticated problems dealing with convective-dispersion phenomena. 2. It was shown that the second order time approximation provides more accurate results than the first order approximation. Investigating the effects of 210 ' second and third order time approximations on flow and third order time approximation on dispersion is highly recommended. 3. It was observed that the accuracy of the numerical results depends on the type and size of ele— ments. Further investigations regarding both size and type of elements will be very useful. LI ST OF REFERENCE S 211 LIST OF REFERENCES Amar, A. C., Ground-water recharge simulation, J. of the Hydraulics Division ASCE, 101(HY9), 1235- 1247, 1975. Bahr, T. G., R. C. Ball, and H. A. Tanner, The Michigan State University Water Quality Management Program, presented at the conference on The Use of Waste- water in the Production of Food and Fibers, Oklahoma City, 1974. Bear, J., Scales of viscous analogy models for ground— water studies, J. of the Hydraulics Division ASCE, 86(HY2), 11423, 1960. Bear, J., On the tensor form of dispersion in porous media, J. of Geophys. Res., 66(4), 1185-1198, 1961a. Bear, J., Some experiments in dispersion, J. of GeOphys. Res., 66(8), 2455-2467, 1961b. Bear, J., Dynamics of Fluids in Porous Media, American Elsevier Publishing Company, Inc., New York, 1972. Bianchi, W. C., and E. E. Haskell, Jr., Field observa— tions of transient groundwater mounds produced by artificial recharge into an unconfined aquifer, Agricultural Research Service W—27, U.S. Department of Agriculture, 27 pp., 1975. Bredehoeft, J. D., and G. F. Pinder, Mass transport in flowing groundwater, Water Resources Res., 9(1), 194-210, 1973. Bruch, J. C., and R. L. Street, Two dimensional dis- persion, J. of the Sanitary Eng. Division ASCE, 93 (8A6), 17—39, 1967. Carnahan, B., H. A. Luther, and J. O. Wilkes, Applied Numerical Methods, John Wiley and Sons, Inc., New York, 1969. 212 213 Cheng, R. T., On the study of convective dispersion equation, State University of New York at Buffalo, 1973. Cheng, R. T., and C. Y. Li, On the solution of transient free surface flow problems in porous media by finite element method, J. of Hydrology, 20, 49-63, 1973. Cooley, R. L., A finite difference method for unsteady flow in variably saturated porous media: appli- cation to a single pumping well, Water Resources Res., 7(6), 1607-1625, 1971. Davis, S. N., and R. J. M. DeWiest, HydrOgeology, John Wiley and Sons, Inc., New York, 1966. de Josselin de Jong, G., and M. J. Bossen, Discussion of paper by J. Bear, on the tensor form of diSper- sion, J. of GeOphys. Res., 66(10), 3623-3624, 1961. Desai, C. S., Seepage analysis of earth bank under draw- down, J. of the Soil Mech. and Found. Division ASCE, 98(SM11), 1143-1162, 1972. Desai, C. S., and J. F. Abel, Introduction to the Finite Element Method, Van Nostrand Reinhold Co., New York, 1972. Doctors, L. J., An application of the finite element technique to boundary value problems of potential flow, Int. J. for Num. Meth. in Eng., 2, 243—252, 1970. Donea, J., On the accuracy of finite element solutions to the transient heat—conduction equation, Int. J. for Num. Meth. in Eng., 8, 103—110, 1974. Douglas, J., D. W. Peaceman, and H. H. Rachford, A method for calculating multidimensional immiscible dis- placement, Trans. AIME, 216, 297-308, 1959. Eagleson, P. 8., Dynamic Hydrology, McGraw-Hill Co., Inc., New York, 1970. Ellis, B. G., The soil as a chemical filter,in Recycling Treated Municipal Wastewater and Sludge Through Forest and CrOpland, ed. by W. E. Sopper and L. T. Kardos, The Pennsylvania State University Press, 1973. 214 Ergatoudis, I., B. M. Irons, and O. C. Zienkiewicz, Curved isoparametric quadrilateral elements for finite element analysis, Int. J. Solids Struc- tures, 4, 31-42, 1968. Fattah, Q. N., Investigation and verification of a model for the dispersion coefficient tensor in flow through anisotropic, homoqeneous, porous media with application to flow from a recharge well through a confined aquifer, Ph.D. thesis, Univ. of Wisconsin, Madison, 1974. France, P. W., Finite element analysis of unconfined seepage problems, Ph.D. thesis, Univ. of Wales, Swansea Glamorgan, England, 1971. France, P. W., C. J. Parekh, J. C. Peters, and C. Taylor, Numerical analysis of free surface seepage prob- lems, J. of the Irrigation and Drainage Division ASCE, 97(IR1), 165-179, 1971. France, P. W., Finite element analysis of three-dimensional groundwater flow problems, J. of Hydrology, 21, 381-398, 1974. Fried, J. J., and M. A. Ckmbarnous, Dispersion in porous media, in Advance in Hydroscience, ed. by Ven Te Chow, 7, 170—280, Academic Press, New York, 1971. Gallagher, R. H., Finite Element Analysis Fundamentals, Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1975. Gardner, A. 0., Jr., D. W. Peaceman, and A. L. Pozzi, Jr., Numerical calculation of multidimensional miscible displacement by the method of charac— teristics, Soc. Petrol. Eng, J., 4(1), 26-36, 1964. Gray, W. G., and G. F. Pinder, Galerkin approximation of the time derivative in the finite element analysis of groundwater flow, Water Resources Res., 10(4), 821-828, 1974. Guymon, G. L., A finite element solution of the one— dimensional diffusion—convection equation, Water Resources Res., 6(1), 204-210, 1970. Guymon, G. L., V. H. Scott, and L. R. Herrmann, A general numerical solution of two-dimensional diffusion— convection equation by the finite element method, Water Resources Res., 6(6), 1611-1617, 1970. ”CD 215 Hantush, M. 8., Growth of a groundwater ridge in response to deep percolation, symposium on Transient Groundwater Hydraulics, 223 pp., Colorado State University, Fort Collins, 1963. Hantush, M. 8., Growth and decay of groundwater-mounds in response to uniform percolation, Water Resources Res., 3(1), 227-234, 1967. Harleman, D. R. F., and R. R. Rumer, Jr., Longitudinal and lateral dispersion in an isotropic porous medium, J. of Fluid Mech., 16(3), 385-394, 1963. Herbert, R., Time variant groundwater flow by resistance network analogues, J. of Hydrology, 6, 237~264, 1968. Hoopes, J. A., and D. R. F. Harleman, Wastewater recharge and dispersion in porous media, J. of the Hydrau- lics Division ASCE, 93(HY5), 51-71, 1967. IBM Application Program, System 360 Scientific Package (360A-CM03X), Version III, Program Manual, H20—0205-3, 137—141, 1968. Isaacson, E., and H. B. Keller, Analysis of Numerical Methods, John Wiley and Sons, Inc., New York, 1966. Javandel, I., and P. A. Witherspoon, Application of the finite element method to transient flow in porous media, Soc. of Petrol. Eng. J., 8(3), 241—252, 1968. Javandel, I., and P. A. Witherspoon, A method of analyz- ing transient fluid flow in multilayered aquifer, Water Resources Res., 5(4), 856—869, 1969. Konikow, L. F., and J. D. Bredehoeft, Modeling flow and chemical quality changes in an irrigated stream- aquifer system, Water Resources Res., 10(3), 546—562, 1974. Marino, M. A., Hele-Shaw model study of the growth and decay of groundwater ridges, J. of Geophys. Res., 72(4), 1195—1205, 1967. Marino, M. A., Water-table fluctuation in semipervious stream—unconfined aquifer system, J. of Hydrology, 19, 43-52, 1973. 216 Marino, M., Longitudinal dispersion in saturated porous media, J. of the Hydraulic Division ASCE, 100(HY1), 151-157, 1974a. Marino, M. A., Rise and decline of the water table induced by vertical recharge, J. of Hydrology, 23, 289- 298, 1974b. Marino, M. A., Water-table fluctuation in response to recharge, J. of the Irrigation and Drainage Division ASCE, 100(IR2), 117-125, 1974c. Marino, M. A., Distribution of contaminants in porous media flow, Water Resources Res., 10(5), 1013— 1018, 1974d. Meissner, U., A mixed finite element model for use in potential flow problems, Int. J. for Num. Meth. in Eng., 6, 467—473, 1973. Moench, A. F., Analytical solutions to the one-dimensional nonlinear equation for flow through porous media, Water Resources Res., 9(5), 1378-1384, 1973. Nalluswami, M., R. A. Longenbaugh, and D. K. Sunada, Finite element method for hydrodynamic dispersion equation with mixed partial derivatives, Water Resources Res., 8(5), 1247-1250, 1972. Neuman, S. P., and P. A. Witherspoon, Finite element method of analyzing steady seepage with a free surface, Water Resources Res., 6(3), 889-897, 1970. Neuman, S. P., and P. A. Witherspoon, Analysis of non— steady flow with a free surface using the finite element method, Water Resources Res., 7(5), 611—623, 1971. Neuman, S. P., Saturated—unsaturated seepage by finite elements, J. of the Hydraulics Division ASCE, 99(HY12), 2233-2250, 1973. Norrie, D. H., and G. deVries, The Finite Element Method, Academic Press, New York, 1973. Ogata, A., and R. B. Banks, A solution of the differen- tial equatiOn of longitudinal dispersion in porous media, Professional paper 4ll—A, U.S. Geological Survey, U.S. Government Printing Office, Washington, D.C., 7 pp., 1961. 217 Peaceman, D. W., and H. H. Rachford, Numerical calcula- tion of multidimensional miscible displacement, Soc. of Petrol. Eng. J., 2(4), 327-339, 1962. Pennington, R. H., Introductory Computer Methods and Numerical Analysis, The Macmillan Co., New York, 1970. Pinder, G. F., and H. H. Cooper, Jr., A numerical tech- nique for calculating the transient position of the saltwater front, Water Resources Res., 6(3), 875-882, 1970. Pinder, G. F., and E. O. Frind, Application of Galer- kin's procedure to aquifer analysis, Water Resources Res., 8(1), 108-120, 1972. Pinder, G. F., A Galerkin-finite element simulation of groundwater contamination on Long Island, New York, Water Resources Res., 9(6), 1657-1669, 1973. Pinder, G. F., E. O. Frind, and S. S. Papadopulos, Functional coefficients in the analysis of ground— water flow, Water Resources Res., 9(1), 222- 226, 1973. Pinder, G. F., and W. G. Gray, Notes on the finite ele- ment method in surface and subsurface hydrology, School of Engineering and Applied Science, Princeton Univ., 1974. Polubarinova-Kochina, P., Theory of Groundwater Movement, Princeton Univ. Press, Princeton, N.J., 1962. Poreh, Michael, The dispersivity tensor in isotropic and axisymmetric mediums, J. of Geophys. Res., 70(16), 3909-3913, 1965. Price, H. S., J. C. Cavendish, and R. S. Varga, Numerical methods of higher order accuracy for diffusion- convection equation, Soc. Petrol. Eng. J., 8(3), 293-303, 1968. Reddell, D. L., and D. K. Sunada, Numerical simulation of dispersion in groundwater aquifers, Hydrol. Pap. 41, 79 pp., Colo. State Univ., Fort Collins, 1970. 218 Robertson, J. B., Digital modeling of radioactive and chemical waste transport in the Snake River Plain Aquifer at the national reactor testing station, Idaho, USGS Open-File Report IDO-22054, U.S. Department of Commerce, Sprifigfield, Virginia, 41 pp., 1974. Rumer, R. R., Longitudinal dispersion in steady and unsteady flow, J. of the Hydraulics Division, 88(HY4), 147-172, 1962. Sandhu, R. S., I. S. Rai, and C. S. Desai, Variable time step analysis of unconfined seepage, in Finite Element Methods in Flow Problems, ed. by J. T. Oden et a1., 573-579, UAH Press, 1974. Scheidegger, A. E., General theory of dispersion in porous media, J. of Geophys. Res., 66(10), 3273-3278, 1961. Segerlind, L. J., Applied Finite Element Analysis, John Wiley and Sons, Inc., New York (in press). Segol, G., G. F. Pinder, and W. G. Gray, A Galerkin- finite element technique for calculating the tran— sient position of the saltwater front, Water Resources Res., 11(2), 343-347, 1975. Shamir, U. Y., and R. F. Harleman, Dispersion in layered porous media, J. of the Hydraulics Division ASCE, 93(HY5), 237—260, 1967. Smith, I. M., R. V. Farraday, and B. A. O'Connor, Rayleigh—Ritz and Galerkin's finite element for diffusion-convection problems, Water Resources Res., 9(3), 598-606, 1973. Taylor, G. S., and N. J. Luthin, Computer methods for transient analysis of water—table aquifers, Water Resources Res., 5(1), 144—152, 1969. Taylor, R. L., and C. B. Brown, Darcy flow solutions with a free surface, J. of the Hydraulics Division ASCE, 93(HY2), 25-33, 1967. Tinsley, P. S., and R. M. Ragan, Investigations of the response of an unconfined aquifer to localized recharge, Technical Report No. 20, Univ. of Maryland, College Park, 80 pp., 1968. 219 Todd, D. K. Groundwater Hydrology: John Wiley and Sons, Inc., New York, 1959. Todsen, M., On the solution of transient free-surface flow problems in porous media by finite- difference methods, J. of Hydrology, 12, 177- 210, 1971. Tseng, M. T., and R. M. Ragan, Behavior of groundwater flow subject to time-varying recharge, Water Resources Res., 9(3), 734-742, 1973. Walton, W. C., Groundwater Resource Evaluation, McGraw- Hill Co., Inc., New York, 1970. Wang, M., and R. T. Cheng, A study of convective- dispersion equation by isoparametric finite ele- ment, J. of Hydrology, 24, 45-56, 1975. Weaver, W., Jr., Computer Programs for Structural Analysis, D. Van Nostrand Co., Inc., Princeton, N.J., 1967. Welty, J. R., C. E. Wicks, and R. E. Wilson, Fundamentals of Momentum Heat and Mass Transfer, John Wiley and Sons, Inc., New York, 1969. Wilson, E. L., and R. E. Nuckell, Application<1ffinite element method to heat conduction analysis, Nucl. Eng. Des., 4(3), 267-286, 1966. Yen, B. C., and C. H. Hsie, Viscous flow model for ground- water movement, Water Resources Res., 8(5), 1299—1306, 1972. Zienkiewicz, O. C., and Y. K. Cheung, Finite elements in the solution of field problems, The Engineer, 220, 507—510, 1965. Zienkiewicz, O. C., P. Mayer,enuiY. K. Cheung, Solution of anisotropic seepage by finite elements, J. Eng. Mech. Division ASCE, 92(EM1), 111—120, 1966. Zienkiewicz, O. C., and C. J. Parekh, Transient field problems, two—dimensional and three—dimensional analysis by iSOparametric finite elements, Int. J. for Num. Meth. in Eng., 2, 61-71, 1970. Zienkiewicz,(D.C., The Finite Element Method in Engineer- ing Science, McGraw-Hill Publishing Co., London, 1971. APPENDICES 220 '—.u.‘l:_..d-—_..'.'.a.-P—_‘vb ...,-___. 2.1:" APPENDIX I FINITE ELEMENT DEVELOPMENT 221 APPENDIX I FINITE ELEMENT DEVELOPMENT 1.1 Introduction In the finite element technique, a continuum is divided into a finite number of subdomains which are called "elements." Each element is designated by "nodes." It is possible to define a functional such that it will describe uniquely the state of a parameter within an element based on its values at the nodes. Polynomials are most commonly used in deriving such functionals, which are termed "shape functions"[Seger1ind, in press]. A detailed formulation of the finite element method is given in the literature, e.g., Zienkiewicz [1971], Norrie and de Vries [1973]. Let the dependent variable C in the domain De be approximated by e (I.1.l) where Nn are the appropriate shape functions defined piecewise, element by element; Cn are the nodal values of C in the discretized domain; and M is the number of 222 223 the degree of freedom (number of nodes in each element). For example, for an element with three nodes _ e 2 C3 where [ ]e and { }e denote a row matrix and a column vector respectively, which contain prOperties of the three nodes associated with one element. In this section the shape functions for the different types of elements used in this study are described, and the numerical integration of element matrices for isoparametric elements is discussed. The integrated element matrices for one-dimensional quadratic and two-dimensional triangular elements are given in Appendix II. Finally the procedure for allocation of a constant line source to the boundary nodes is given. I.2 Types of Finite Elements and Their Shape Functions One cannot subdivide a continuum into elements without first knowing what general shapes are permissible. In this part some of the more common finite elements which are used in the analysis of flow and dispersion phenomena are given. For the derivation of the shape functions, the reader is referred to Zienkiewicz [1971]. 224 £.2.1 One-Dimensional Element a. Simple e1ement.--The simplest one-dimensional element has two nodes, one at each end (Figure I-l). Figure I-l.--Simp1e One-Dimensional Finite Element. The shape functions for this element are x E] (1....) b. Quadratic e1ement.--This element has three nodes, two nodes at each end and one in the center of the element (Figure I-2). L————L—-l a; 2 I '——-’-—X Figure I—2.--One-Dimensiona1 Quadratic Element. H-“_—‘-.L:l_l"'-E '- 225 The shape functions for a quadratic element are N=[———) [Hg-1 [1 - 35] (1.2.2) c. Quadratic element for axisymmetric case.-- This element also has three nodes, but the shape func- tions are written with respect to the origin of the global coordinate (Figure I-3). («_———-4 11.4 rzj r3 F. Figure I-3.--Quadratic Element for One-Dimensional Axisymmetric. The shape functions can be written N1 = 2(r - r2)(r - r3)/L2 N — 4 - ) - 2 2 — (r rl (r r3)/L N3 = 2(r - rl)(r - r2)/L2 (1.2.3) 226 1.2.2 Two-Dimensional Simplex Element The two-dimensional simplex element is the tri— angle shown in Figure I-4. The triangular element has rightly become more popular due to the ease with which the subdivision can be graded and the boundary shapes approximated. The evaluation of the element matrices is simple. #3, (X31173) Centroid (X21172) ’X Figure I-4.--Simp1e Triangular Element. The shape functions for the triangular element are given below: _1 Nl_7\— N 227 N = —1— (a + b x + c y) (1.2.4) 3 2 3 3 3 a1 = X231’3‘X3Y2 ’ b1:3’2‘3’3 ’ C1="3‘X2 a ==xy -y)< ;b =y -y ;c =x -x 2 3 1 31 2 3 1 2 1 3 a3 = xly2-X2yl ; b3==yl--y2 ; C3==x2--Xl (I.2.5) _ l A — 3(le1 + b2x2 + b3x3) = area of triangle (1.2.6) 1.2.3 Two-Dimensional Iso- parametric Elements The use of a curvilinear coordinate system has definite advantages when considering two- and three- dimensional elements, because it allows the boundaries of these elements to be distorted. a. Linear quadrilateral isoparametric e1ement.-- Consider Figure I—5, and define E,Iqsuch that —1 fi E i 1 and -1 i n i 1. The shape functions for the linear quadrilateral element are (1-6) (l-n) and N : 228 ,5 X Global 7] (-l,+l)¢4 j!) (+1,+1) ""'E 6i ii (-l,-l) (+l,-l) Local Figure I—5.-—Linear Quadrilateral Isoparametric Finite Element on Global and Local Planes. — “- rd- whigflfi' 229 1 N2 = I (1 + a) (1 - n) N—-1—(1+t:)(1+) 3‘4 ” N=l-(l-E)(l+) (127) 4 4 n .. b. Quadratic quadrilateral element.--The shape functions for the quadratic quadrilateral element (Figure I-6) are given [Zienkiewicz 1971, p. 109]: _ 1 . _ l _ 2 _ N1 - " Z (17g) (l-n) (§+T)+l) I N2 - 2 (l E ) (1 T1) N3 = %-<1+a> (l-n) (E-n-l) ; N4 = %-(l-n2) (1+a) N5 = §-(1+§> <1+n> (€+n-l) ; N6 = %-(1-g2) (1+n) N = l-(l—E) (1+n) (-€+n-l) - N = l-<1-n2) <1—a> (1: 2 8) 7 4 ' 8 2 ° ° c. Cubic quadrilateral element.--The shape functions for the cubic quadrilateral element (Figure I-7) are given [Zienkiewicz 1971, p. 109]: 1 g, (1 — a) (1 - n) {-10 + 9 (£2 + n2)] _ 9 2 N2 — 32— (l - 0) (l ‘ E ) (l - 3E) N =—9—(1-n)(1-g2)(1+3g) LA) w 2 230 Figure I-6.--Quadratic Quadrilateral Element. 7 8 1 9 + n 6 11 _____ ”"— a 12 5 \ 1 ' A 2 3 4 Figure I-7.--Cubic Quadrilateral Element. 231 __ 1 _ _ 2 2 N4 - 35 (l + E) (1 n) [ 10 + 9 (E + n )1 _ 9 _ 2 _ N5 — 32— (1 + E) (1 n ) (1 3n) _ 9 2 N6 - g‘z' (1 + 5) (l ' T] ) (l + 3H) _ 1 _ 2 2 N7 _ 33 (1 + E) (l + n) I 10 + 9 (E + n )] _ 9 2 N8—3—2—(1+n)(l-€)(l+3€) _ 9 _ 2 _ N9 - 37 (1 + 71) (l E) (l 32:) _ l 2 2 N10 — 3—2 (1 "' E) (l + Tl) ['10 + 9 (E + n )] N =3—(1—a) (1-n2) (1+3n) 11 32 N = i <1 - a) (1 - n2) (1 — 3n) (I-2—9) 12 32 1.3 Numerical Integration Usually the element matrices are in the integral form which has to be evaluated. For example, the element matrices for fow in a two-dimensional horizontal plane are: e aNk BNn [B] = ge rij 5;; 5;? dD (4.2.7a) [H]e '{F}e = Integration of 232 Ie s Nk Nn dD D ée Nk Q2 d5 - D fe P Nk dD (4.2.7b) (4.27c) Equation (4.2.7) or similar equations for simple elements is straightforward, and some of the inte- grated forms are given in Appendix II. For isoparametric elements, shape functions are described in the local (E, n) coordinates, but Equation (4.2.7) is written in global (x, y) coordinates. To perform the transformation of the shape function derivatives BNk/BX and aNk/By, the following relationship is used. 0} Nk Q.) “I W [w ‘< L In which [J] is 8x 36 RH = 8x am where x1, x2, . [J]— 0) ml W I W i [m J the Jacobian matrix: _— fix a 22 an a—4 nodal coordinates. m1 - 2 H O) \f‘ ‘1 80 .— Another transformation BN 2 3E 8N2 n 8 ., XM and yl, y2, . M 8n_ _ 8N BS EN ° I (1.3.1) _ y1 Y 2 (1.3.2) Yb: yM are the required is 233 the replacement of the element of area, dxdy(dD), by the expression: dD(e) dxdy = det [J] dEdn (1.3.3) The limits of integration in the local coordinate system become -1 and +1. For example, Equation (4.2.7a) will change to +1 +1 , Be=f I T 33324.11 3%. kn xx'ax 8x yy3y 3y -1 -1 det EH dndE (1.3.4) Similar expressions are developed for the remaining terms of Equation (4.2.7), Equation (4.3.13), and Equation (4.4.7), etc. Equation (1.3.4) and similar equations are inte- grated by the Gaussian quadrature integration technique [e.g., see Zienkiewicz 1971, pp. 144—149]. For a poly- nomial of degree 2n = 1, the number of sampling points will be n. 1.4 Parameter Definition The parameters such as hydraulic conductivity, dispersion coefficients, storage coefficient, resistivity of aquitard, etc., can be specified either for each node or each element. Because the numerical solution requires slightly less computation time when parameters are assumed constant over an element, usually the parameters 234 such as storage coefficient which do not change con- siderably are specified for each element. Since velocity components vary within the element, dispersion coefficients are specified at each node. 1.5 Allocation of a Constant Line Source to Boundary Nodes Dirichlet (specified head or concentration) and Neuman (specified flux) conditions are the two boundary conditions generally encountered in field problems. Introducing the Dirichlet boundary condition was dis- cussed in Section 6.1.2. When the normal flux qn is assumed constant along an element face of length L, the integration of £8 qusds (1.5.1) will result in the constant flux to allocated at the boundary nodes. The results for elements with two, three, and four nodes along the boundary line are given in Figure (I—8). One—dimensional linear, quadratic, and cub shape functions are used in this integration. 235 .mucmEon mo mommsm pcmumMMHo cfl mmpoz mumocsom ou condom mcflq pcmumcoo 6 mo coflumooaamln.le mnsmflm ncmfimam vanso unmaoam OHUMHUMSO ucmEmHm Hammad m m ...m... III. III. a m Alull qcvm m m t II II .Hch 4 All scam m m ‘ N ‘II I and CT ace 0 AI gag APPENDIX II INTEGRATED ELEMENT MATRICES FOR ONE- DIMENSIONAL QUADRATIC AND TWO- DIMENS IONAL TRIANGULAR ELEMENTS 236 —I--TI.:T_’ I W. APPENDIX II INTEGRATED ELEMENT MATRICES FOR ONE- DIMENSIONAL QUADRATIC AND TWO- DIMENSIONAL TRIANGULAR ELEMENTS In this section some of the integrated terms that are used in the construction of the element matrices for flow and dispersion phenomena are provided. 11.1 One-Dimensional Quadratic Element The one-dimensional quadratic element is depicted in Figure 1—2. Its shape functions are given by Equation (I.2.2), and can be written as 2 Nl‘l‘EEJ'giz" L L 4x 4X2 N2=—"7 L L 2 N32357—5 (11.1.1) L L 237 332:1- 8x L 3N3 _ 4x 8x L2 238 (11.1.2) In Equations (11.1.1L (11.1.2), and in the following equations, L represents the length of an element. 1. I. 8Nn EE&.dX 8x 8x I NnNk dx RIP I. 8N2 x 8 8Nl x 8 8Nl x 8 8N2 x 8 8N3 x 8 (11.1.3) (II.l.4) 3N Nan 53?— dx I 3. 0\|!—‘ S'an8f8x x I Bhd 239 1 2 Nf%.5§‘ Nflkzgi‘ 8N SN 1 2 N2N1 8_x. N2% 3‘? 8N 8N 1 2 fiffié?‘ N955} F10 8 -fl -6 o 6 _ 1 -8 19‘ (; 8N1 N 8N2 1 5x 1 8x N 8N1 N 8N2 2 5x 2 8x N 8N1 N 8N2 _3 8x 3 8x F—3 ~17 -4 4 _1 -4 1. [; 8N1 8N1 N 8N1 1 8x 8x 2 8x 1 8x 8x 2 8x __1 3x 8x 2 8x aug' N1MB x 8N3 N21"3a'x 51" M 3313. 3 3 8x_)( (11.1.5) N BN3 1 8x 8N3 N2 3;— dx N BN3 3 8Xd (11.1.6) ‘7 8N2 N 8Nl 3N3— 8X 3 8x x 8x 3 8x 8x 8N2 N 8N—3 8N3 8x 3 8x 8x_J 37 __1__ ‘ 30L 44 _ 7 _ N1N1N1 '2 6. f Nk(Nn) dx - £ NZNlNl _P3N1N1 —39 _ L “Z§’ 20 _:3 N1 7 f Nk dx — f N2 X N3 240 -32 16 192 16 7 -44 3zj 1N2N2 2N2N2 3N2N2 -3 20 39 1 L g 4 1 (11.1.7) N1N3N3 N2N3N3 dx 1x13113113-J (11.1.8) (11.1.9) 11.2 Two-Dimensional Simplex Element The two-dimensional simplex element shown in Figure I-4 and the shape functions and related terms are defined in Section I.2.2. The first derivatives of the shape function with respect to x and y are ..b Dh‘ Q) “I ..b J}. Dh‘ 241 1, 2, 3 (11.2.1) blb2 blbgl bzbz b2b3 (11.2.2) b2b3 b3b%_ Clc2 ClC3 c2c2 c2c3 (11.2.3) C203 °3Cq_ blc2 blcil b2c2 b2c3 (11.2.4) b302 b3cg_ b2C1 b301— b2C2 b3c2 (11.2.5) b2c3 b3c3‘ 5. 242 BN . n . Evaluation of i Nk 5§—~dA. —' "1 N 3N1 N 3N2 N BN3 1 5x 1 5x 1 5x 8N 3N 3N 3N3 8N EN EN N ——$ N 2 N 3 L_3 3x 3 3X 3 5x_‘ N1b1 N1102 N1133 1 = 7K f szl sz2 sz3 dA (II.2.6) A N3bl N3b2 N3b%J Each term in Equation (11.2.6) can be integrated separately, e.g., l _ l 21 f Nlbl(fl\_ __§.f (al-+blx-+cly) bl dx dy (11.2.7) A 4AA f dx dy = A = area of triangle A f X dX dy = f y dx dy = 0 A A Then equation (11.2.7) is reduced to ' ‘Cv'JVVn.- ...,“ L1: 1‘ ; . . 243 _Lfab 4A2 l 1 dx dy = (11.2.8) If it is assumed that the origin of the coordinates is taken at the centroid of the element, i.e., 3 = 0 y + y + y l 2 3 = 0 (11.2.9) 3 Th - 39 - — [2' k ' 1971 E 4 8] en a1 — 3 — 32 -‘ a3 len erCZ , q. o . Substituting the value of a into Equation (11.2.8) 1 yields albl _ b_1 4A ‘ 6 1 2 3 aNn 1 i Nk 51— dA = g- bl b2 b3 (11.2.10) b b b J_1 2 {J Note that the assumption leading to Equation (11.2.9) is automatically satisfied in the computer program, regardless of whether local or global coordinates are used. 244 Nlcl NlCZ Nlc3 3Nn 1 6. I}; Nk Ey— dA = fl; NZCl N2C2 N2C3 .P3C1 N3c2 N3c%_ r -( Cl C2 C3 = i c c c (11 2 ll) 6 1 2 3 ' ° f1 C2 C3‘ )2 1 1‘ ._ A 7. £ Nk Nn dA — 12 1 2 1 (11.2.12) 8. Evaluation of g N qn ds: 1f qn is the constant flux k along the line 2—3, as shown in the following figure, then where L is the length of line 2-3. APPENDIX III DEVELOPED COMPUTER PROGRAMS 245 APPENDIX III DEVELOPED COMPUTER PROGRAMS For the numerical investigation of this thesis, several computer programs have been developed for use on a CDC 6500 computer with FORTRAN IV Extended language. The finite element technique is used in formulating all of the programs. Some of these programs will be docu- mented and available. The more pertinent prOgrams are as follows: 1. One-dimensional plane or axisymmetric medium; a. flow b. tracer movement with uniform or transient flow Two-dimensional horizontal plane medium; a. flow with mixed elements (triangu- lar or quadrilateral elements) b. dispersion with simplex triangular elements c. dispersion with quadrilateral ele- ments Two-dimensional vertical plane medium with transient phreatic surface; a. flow MNT quadrilateral alements b. flow FNT linear quadrilateral elements c. dispersion linear quadrilateral elements 246 247 In all programs, options are provided to obtain the velocity vector either directly or simultaneously. Different G-values (see Chapter V) can be used in approximation of time-dependent functions, with a first or second order time approximation available for the convective-dispersion solution. Uniform, steady-state, or transient flow can be used in the calculation of mass-transfer phenomena. A M'Tllifillflfil1fifllfllflth[[fiflflTLfitfiilMWflTflfififilfl