vvEPA United States Environmental Protection Agency National Risk Management Research Laboratory Ada, OK 74820 Research and Development EPA/600/S-00/001 May 2000 ENVIRONMENTAL RESEARCH BRIEF Analytic Element Modeling of Ground-Water Flow and High Performance Computing H.M. Haitjema1, V.A. Kelson2, and K.H. Luther3 Abstract Several advances in the analytic element method have been made to enhance its performance and facilitate three- dimensional ground-waterflow modeling in a regional aquifer setting. First, a new public domain modular code (ModAEM) has been developed for modeling ground-water flow on a regional scale by use of the AEM method. The code can be compiled as Fortran 90 for use on a single processor machine, such as a PC, or it can be compiled with a high performance Fortran (HPF) compiler for use on a massively parallel processor machine (MPP). The performance increase from the MPP machine enhances the applicability of the code to large regional aquifers including local three-dimensional flow effects. Second, new algorithms were developed, coded, and tested to expand an earlier solution to a partially penetrating well in aconfined aquifer to cases of unconfined flow. The new three-dimensional solution to a partially penetrating well includes the computation of the phreatic surface and, when appropriate, a seepage face. The solution also includes local aquifer stratification. The newSD solution is suitable for embedding in a Dupuit-Forchheimer model based on the analytic element method. However, the 3D solution is mathematically involved and computationally demanding; it has not been included in ModAEM. Introduction The analytic element method (AEM) has been shown to be uniquely suited for modeling large regional aquifers. The method employs superposition of analytic functions to create an ' School of Public and Environmental Affairs, Indiana University, Bloomington, IN 47405. 2 Now at Wittman Hydro Planning Associates, Inc., Ellettsville, IN 47429. 3 Now at Department of Mathematics and Computer Science, Valparaiso University, Valparaiso, IN 46383. approximate, but analytic solution to the regional ground-water flowproblem. The analyticfunctions are called "analyticelements," each of which represents some flow feature in the aquifer. For instance, a line sink is used to model a section of a stream. Since only the hydrologic features (streams, lakes, etc.)arediscretized in the analytic element method and not the flow domain in between these features (no finite difference grid or finite element mesh), the AEM is relatively insensitive to scale. Regionally, a rather coarse representation of streams and lakes may be employed, while locallythe hydrologicfeatures can be represented in more detail in order to obtain the required model precision. In fact, three-dimensional solutions can be embedded in the regional Dupuit-Forchheimer model (horizontal flow model) to further enhance the realism of the model in the area of interest. This flexibility in scale and model resolution makes the analytic element method efficient for modeling large regional aquifers. Priorto this research, local three-dimensional flow in the analytic element method was mostly limited to flow near a partially penetrating well in a confined aquifer. As part of the current research project, the solution to a partially penetrating well has been expanded to cases of unconfined flow. Thethree-dimensional solution to flow near a well in an unconfined aquiferis significantly more complex than the earlier confined flow solution due to the presence of a phreatic surface and, in many cases, a seepage face along the upper portions of the well bore. In addition, embedding the 3D unconfined solution is not a simple matter of superposition, as it was for the case of confined flow. Large regional models with much local detail require many analytic elements which all contribute to the flow. This leads to many calculations to arrive at the head or flow vector at a point in the aquifer. The computational effort is particularly large when some of the analytic elements represent three-dimensional flow, like the one for a partially penetrating well. To facilitate large and complex models that remain computationally fast, we tested the suitability of the AEM for parallel processing. In principle, the AEM ------- is highly suitable for parallel processing as the evaluation of the head orthe flow vector at a point involves very many independent calculations that can be executed in parallel. As part of this research, we developed a new public domain code which can be compiled for a single processor machine (e.g., PC) or an MPP machine. The new code, ModAEM, is modular in its design to facilitate future expansions by others. The current research is reported in detail in two different Ph.D. theses. The development and testing of the three-dimensional solution to a partially penetrating well in an unconfined aquifer is described in "Analytic Solutions to Three-DimensionalUnconfined Groundwater Flow near Wells" (Luther, 1998). The development and testing of the code ModAEM and the testing of the AEM on an MPP machine is found in "Practical Advances in Groundwater Modeling with Analytic Elements" (Kelson, 1998). The latter thesis also reports on two alternative procedures to locally improve the realism of the regional analytic element method. The first procedure is the extraction of a local U.S. Geological Survey MODFLOW model out of a regional analytic element model. The second procedure is to model locally three-dimensional flow in an approximate fashion using only two-dimensional functions. The idea is to subdivide the aquifer into layers and model Dupuit- Forchheimer flow in each layer, while modeling vertical flow based on leakage between the layers. Both theses are published by the Graduate School of Indiana University. An electronic copy can be downloaded from the resources library of the OpenAEM webpage(http://www.indiana.edu/~aem/). An additional reference for the 3D well research is Luther and Haitjema [1999]. Three-Dimensional Flow near a Partially Penetrating Well Most aquifers of interest for water resource purposes are relatively thin compared to their horizontal (lateral) extent. Under these circumstances, resistance to vertical flow may be ignored. This is known as the Dupuit-Forchheimer assumption. Locally, however, for instance near a partially penetrating well, resistance to vertical flow is important and a fully three-dimensional solution to the flow problem may be needed. An analytic solution to flow near a partially penetrating well in a confined aquifer has been implemented before in the code GFLOW [Haitjema, 1995]. This solution uses a vertical line sink at the center of the well bore to model the flow into the well. The line sink has a variable sink density distribution to account for the inflow variation along the well bore. In fact, the sink density near the ends of the line sink (and thus well bore) is infinite. The no-flow conditions along the (horizontal) aquifer top and bottom are realized by use of a truncated series of image line sinks, in combination with two semi-infinte line sinks [Haitjema and Kraemer, 1988]. The complete three-dimensional solution is simply superimposed onto the many two-dimensional solutions that generate the regional Dupuit-Forchheimer flow field. The resulting solution is trulythree-dimensional nearthe well and is a Dupuit-Forchheimer approximation (2D) everywhere else. As part of the current research project, a three-dimensional solution for flow near a partially penetrating well in an unconfined aquifer was developed and tested. This solution is substantially more involved than the earlier one for confined flow. In the first place, the unconfined solution involves a phreatic surface near the well whose shape and location is not a priori known. In many cases there is also a seepage face along the upper part of the well bore, the extent (indeed existence) of which is not known in advance. Finally, the three-dimensional solution cannot be simply superimposed onto the two-dimensional solutions, as was done forthe case of confined flow. To further enhance the realism of the solution, we also implemented aquifer stratification in the vicinity of the well. Representation of the Phreatic Surface We set out to create an analytic solution for flow near a partially penetrating well that is everywhere regular and differentiate, including at the boundaries of the flow domain. To accomplish this, the generating functions for flow, which often contain singularities, were placed outside the actual flow domain. Thus, instead of representing the phreatic surface by a facet surface of single or double layers (sink or doublet panels), as is done in the classical boundary element method, we introduced a continuous scalar function %(x,y,z), which is zero along the phreatic surface. This auxiliary function x may be viewed as the three-dimensional version of the Zhukovski function used in complex potential theory to represent a phreatic surface in the vertical plane. In words, the function % is the hydraulic conductivity times the pressure head, and can be expressed in terms of the specific discharge potential (hydraulic conductivity times the head) minus the hydraulic conductivity times the elevation of the point where the potential is evaluated. The phreatic surface is then defined by the surface where the function % is equal to zero, since that is where the pore pressure (pressure head) is zero. We placed point sinks and sink rings of different radii above the phreatic surface to enforce the condition that the flow normal to the phreatic surface is zero. The sink densities of the line sink along the well axis and the auxiliary sink distributions above the phreatic surface are determined in such a manner as to satisfy the following conditions. Along the well bore, equipotential conditions are imposed, and along the %=0 surface, no flow conditions or a specified recharge rate is enforced. These conditions are met nearly exactly at a number of collocation points on these surfaces, while they are approximate in between. Due to the non-linear nature of the equations that must be solved to arrive at the sink density parameters, an iterative solution procedure is employed by which the location of the phreatic surface is adjusted during iterations. In Figure 1, streamlines toward three wells are depicted in a 3D diagram. All streamlines start from a vertical array of points, located in the vertical sheet to the right in the figure. The horizontal grid is the aquifer bottom, while the thin shaded zone above the wells is a side view of the (slightly) curved phreatic surface. The vertical arrays of points above the wells are point sinks outside the flow domain to enforce no flow conditions along the phreatic surface. Some streamlines are shown above the phreatic surface, which enter these point sinks while emanating from points just above the phreatic surface. The latter streamlines are, of course, outside the flow domain and shown only to illustrate the no-flow conditions along the phreatic surface. Representation of the Seepage Face If the head in the well occurs below the top of the well bore, a seepage face may develop along the upper portion of the well bore. The pressure head along the seepage face is zero, but the flow normal to the seepage face is not zero. Consequently, only the condition %= 0 is enforced along the seepage face. The line sink along the well bore opposite the seepage face is used to enforce this condition in a number of collocation points along the seepage face. The solution procedure for the sink density parameters is iterative since neither the location of the phreatic surface northe location of the seepage face is known in advance. In Figure2, we compared our solution to othersolutions presented in the literature. The problem is a sand box model of flow toward a well in an unconfined aquifer as published by Hall in 1955. Our analytic element solution (denoted by AEM in Figure2) compares ------- , , 1 ....._; i.._.j.—• :__^. I t—L ., , i 4-MW4X ^£=3~f-::±::-i- T-^---V-^r-^^--\----v-^-V^---\^Av1V^: Figure 1. Streamlines toward three partially penetrating wells in an unconfined aquifer. Dots above the phreatic surface are point sinks to enforce no-flow conditions along the phreatic surface. AEM, BIEM, FEM solutions Dupuit-Forchheimer surface Figure 2. Comparison with the sand box experiment of Hall (1955). Current analytic solution is denoted by the AEM solution. ------- well with a published boundary integral equation model (BIEM) [Liggett and Liu, 1983] and a published finite element model (FEM) [Neuman and Witherspoon, 1970]. All three numerical solutions, however, exhibit a phreatic surface that is lower than observed in the experiment. It is noted that published numerical solutions that include unsaturated flow do match more closely with Hall's experiment [Hall,1955]. The horizontal lines inside the well bore, opposite the seepage face, and the horizontal lines above the phreatic surface, represent sink rings used for the AEM solution. Sink rings were used inside the well rather than line sinks because of the relatively large well diameter. Representation of Aquifer Stratification Particularly useful are solutions to flow near wells in a stratified aquifer. Such solutions exhibit the preferential paths of ground- water flow in the more permeable strata. Conceptually, the aquifer is subdivided into a number of horizontal layers of varying thickness and hydraulic conductivity: the aquifer strata. Consequently, the discharge potential (hydraulic conductivity times the head) jumps across the interfaces between the layers. The jump in the potential across an interface is modeled with double layers, which are dipole distributions over the interface with their axis normal to the interface (also referred to as doublet distributions). We used doublet discs with a radially varying doublet density to account for the radial variation in the potential jump. In Figure 3, some streamlines are depicted in a cross-section over a well in a stratified aquifer. There are three strata with a conductivity of 10, 50, and 10 feet per day, respectively. The solutions for stratified aquifers are not as robust as those for homogeneous aquifers. Since we use radially varying doublet densities we can only model radially symmetric flow, hence flow toward a single well. Also, contrasts in the hydraulic conductivity in excess of a factor of 10 yield inaccurate solutions for the velocities near the interfaces. Finally, when the phreatic surface intersects the interface, some local inaccuracies in the flow field near the intersection may occur. Embedding 3D Solutions in a Dupuit-Forchheimer Model Under confined flow conditions, the three-dimensional solution to a well reduces to atwo-dimensional (Dupuit-Forchheimer) solution at some distance from the well (3 to 4 times the aquifer thickness). Near the well, where the flow is three-dimensional in nature, the effect of remote flow features (e.g., streams) is two-dimensional due to the horizontal aquifer boundaries. Consequently, the two- and three-dimensional solutions in confined aquifers can simply be superimposed. In so doing, three-dimensional flow is taken into account where it occurs, near partially penetrating wells, while all other flow features are represented by two-dimensional functions. This makes for very efficient modeling of three- = ft/d k = 50 ft/d Figure 3. Streamlines toward a partially penetrating well in a stratified unconfined aquifer. The two interfaces between different hydraulic conductivities are modeled with doublet discs with a radially varying doublet density. ------- dimensional flow in regional aquifers. Under unconfined flow conditions, however, the presence of a three-dimensional phreatic surface near the well prevents the simple superposition of two- and three-dimensional solutions. Under conditions of unconfined flow, three-dimensional solutions may be embedded in a Dupuit- Forchheimer model by defining a 3D zone near one or more partially penetrating wells. Heads obtained from the Dupuit- Forchheimer solution are defined on the boundary of that zone, which contains the partially penetrating well(s), as fully penetrating with the same pumping rate(s). The solution procedure is iterative in case the pumping rate of a well is not known in advance; for instance, when the head at a well is specified instead of the pumping rate. For a well with a known (specified) pumping rate, no iterations are needed to couple the 2D and 3D solutions. Efficiency of the 3D Solution fora Partially Penetrating Well The analytic element solution for a partially penetrating well, developed as part of this research project, issuperiorto previous solutions in that the discharge potential (thus the piezometric head) is everywhere differentiable, including at the well bore and at the phreatic surface. However, the solution requires many analytic functions (line sinks, point sinks, sink rings, doublet discs, etc.) to enforce the boundary conditions at collocation points along the well bore, phreatic surface, and interfaces between aquifer strata. In fact, to obtain accurate solutions to non-axisymmetric flow, a least square method was used to approximate the conditions at many more control points than the number of degrees of freedom (e.g., sink density parameters) available at the auxiliary functions. To maximize performance, some experimentation with the number, nature, and location of the auxiliary functions was required. In summary, the solutions are involved, computationally intensive, and not always straightforward in their implementation. Analytic Elements on a Massively Parallel Processor Machine Analytic element models perform three computationally intensive tasks: (1) generating a set of equations to solve for the strength parameters of the analytic elements, (2) solving these equations, and (3) calculating heads and orfluxes at points in the aquifer. We tested the performance increase foreachofthese tasks separately on Indiana University's Intel-Paragon MPP machine with 64 processors. Tasks (1) and (3) consist of a large number of independent calculations: the computation of the influence of each analytic element at a number of points in the aquifer. All of these computations can be carried out in parallel with little or no communication between processors (no overhead). Consequently, we found that the computational speed scaled up linearly with the number of processors employed. Task (2), the solution of a set of equations with a full coefficient matrix, can also be speeded up using parallel processors, but not without adding overhead. The computational speed up of the matrix solution process, therefore, appeared to be sub-linear. Yet the overall speed up of analytic element models on parallel processors is nearly linear, since the matrix solution procedure accounts for only a small portion of the overall computational effort during analytic element modeling. A Modular Parallel Analytic Element Code We developed a new public domain analytic element code, ModAEM, with the following characteristics: (1) a modular design similar to MODFLOW, (2) containing a limited set of basic but versatile analytic elements, (3) designed to run on both a single processor machine (e.g., PC orworkstation) and an MPP machine. The code requires a GUI for efficient operation. An example GUI is the public domain code WftAEM (U.S. EPA, 2000), which has been developed for capture zone delineation purposes under EPA's wellhead protection program. W/?AEM is a Windows 957 98/NT program, which uses ModAEM as its computational engine. ModAEM has been written in Fortran 90 (F90) and has been expanded with directives for High Performance Fortran (HPF). When compiled using a standard F90 compiler it will create a single processor code, while when compiled with HPF it will run on an MPP machine. ModAEM is intended as a publicly available basis forthe implementation of advances in the analytic element method. ModAEM does not currently include any three- dimensional analytic elements, such as the partially penetrating well discussed above. Summary and Conclusions We developed and tested a new three-dimensional solution for flow toward a partially penetrating well in an unconfined aquifer. The solution includes a phreatic surface and seepage face near thewell. Forthe axisymmetric case of flowtoward asinglewellwe also included aquifer stratification. The analytic solution employs singularity distributions outside the flow domain. Consequently, the resulting discharge potential is everywhere differentiable, including at the domain boundaries (e.g., well bore, seepage face, and phreatic surface). The solution, however, is complex and computationally intensive. To increase the computational efficiency of the AEM, we tested the AEM performance on an MPP machine. The majority of the computations scaled up linearly with the number of processors used due to the independent nature of most AEM computations. The nearly linearoverall performance increase is unusual. Analytic element models appear particularly suitable for MPP architectures. We developed a public domain modular analytic element code, ModAEM, which will run on both a serial machine (e.g., PC) and an MPP machine. The code uses a limited number of versatile analytic elements to perform most groundwater modeling tasks. ModAEM is written in F90 and contains directives for HPF compilers to adapt it to MPP machines. We found that the three-dimensional solution to the partially penetrating well in an unconfined aquifer is involved and computationally intensive. Development of additional three- dimensional solutions, for instance, for a shallow creek or lake, may be expected to be similarly complicated. Kelson [1998] proposes a technique forapproximatethree-dimensional solutions in analytic element models, which is conceptually similar to the procedure used in MODFLOW. The aquifer is subdivided into a number of horizontal layers with resistance to vertical flow represented as resistance to flow between the layers ("quasi three-dimensional flow"). Each layer supports Dupuit-Forchheimer flow, hence only two-dimensional solutions to flow are being used. The approach has been demonstrated for relatively high vertical resistances, such as occur when the vertical hydraulic conductivity is less than one-tenth of the horizontal conductivity. However, it is uncertain if accurate solutions can be obtained under more homogeneous conditions. Acknowledgments Dr. Stephen Kraemer served as the U.S. EPA Project Officer. We thank the following reviewers for their valuable comments: Dr. David Steward of Kansas State University; Mr. Brent Swartz of Lockheed Martin Technical Services, Inc.; Dr. Charlie Fitts of the University of Southern Maine; and Dr. Mark Bakker of the University of Nebraska-Omaha. Quality Assurance Statement This projectwas conducted underan approved quality assurance program plan and the procedures specified therein were used. ------- Information on the plan and documentation are available fromthe Principal Investigator. Notice The U.S. Environmental Protection Agency through its Office of Research and Development partially funded and collaborated in the research described here under Cooperative Agreement CR- 823637 with Indiana University. The work has been subjected to the Agency's peer and administrative review and it has been approved for publication as an EPA document. Mention of trade names or commercial products does not constitute endorsement or recommendation for use. References Haitjema, H.M., 1995. Analytic Element Modeling of Groundwater Flow. Academic Press. Haitjema, H.M., andS.R. Kraemer, 1988. A new analytic element function for modeling partially penetrating wells, Water Resources Research, 24(5):683-690. Hall, H.P., 1955. An investigation of steady flow toward a gravity well, La Houille Blanche, 10, 8-35. Kelson, V.A., 1998. Practical Advances in GroundwaterModeling with Analytic Elements, Ph.D. dissertation, Indiana University, Bloomington (H. Haitjema, advisor). Liggett, J.A., and P.L.-F. Liu, 1983. The Boundary Integral Equation Method for Porous Media Flow. George Allen and Unwin Ltd., London, U.K. Luther, K.H., 1998. Analytic Solutions to Three-Dimensional Unconfined GroundwaterFlownearWells, Ph.D. dissertation, Indiana University, Bloomington (H. Haitjema, advisor). Luther, K., and H.M. Haitjema, 1999. An analytic element solution to unconfined flow near partially penetrating wells, Journal of Hydrology, 226:197-203. Neuman, S.P., and P.A. Witherspoon, 1970. Finite element method of analyzing steady seepage with a free surface, Water Resources Research, 6(3):889-897. U.S. EPA, 2000. Working with W/7AEM2000: source water assessmentforaglacial outwash aquifer, Vincennes, Indiana, EPA/600/R-00/022. ------- United States Environmental Protection Agency Center for Environmental Research Information Cincinnati, OH 45268 Official Business Penalty for Private Use $300 EPA/600/S-00/001 BULK RATE POSTAGE & FEES PAID EPA PERMIT No. G-35 ------- |