A computer program for simulating geohydrologic systems in three dimensions

This document is directed toward individuals who wish to use a computer program to simulate ground-water flow in three dimensions. The strongly implicit procedure (SIP) numerical method used to solve the set of simultaneous equations is based on the computer programs published by Trescott (1975) and Trescott and Larson (1976). New data processing techniques and program input and output options are emphasized. Theoretical development is minimized, referring to other sources for this information. The aquifer system to be modeled may be heterogeneous and anisotropic, and may include both artesian and water-table conditions. Systems which consist of well defined alternating layers of highly permeable and poorly permeable material may be represented by a sequence of equations for two dimensional flow in each of the highly permeable units. Boundaries where head or flux is user-specified may be irregularly shaped. The program also allows the user to represent streams as limited-source boundaries when the stream flow is small in relation to the hydraulic stress on the system. The data-processing techniques relating to "cube" input and output, to swapping of layers, to restarting of simulation, to free-format NAMELIST input, to the details of each subroutine's logic, and to the overlay program structure are discussed. The program is capable of processing large models that might overflow computer memories with conventional programs. Suggestions are given for the modeler wishing to convert this code to computers not manufactured by Control Data Corporation. The program is written completely in the FLECS structured programming language, and produces standard FORTRAN source statements. Detailed instructions for selecting program options, for initializing the data arrays, for defining "cube" output lists and maps, and for plotting hydrographs of calculated and observed heads and/or drawdowns are provided. Output may be restricted to those nodes of particular interest, thereby reducing the volumes of printout for modelers, which may be critical when working at remote terminals. "Cube" input commands allow the modeler to set aquifer parameters and initialize the model with very few input records. Appendixes provide instructions to compile the program, definitions and cross-references for program variables, summary of the FLECS structured FORTRAN programming language, listings of the FLECS and FORTRAN source code, and samples of input and output for example simulations. INTRODUCTION The computer program documented in this report simulates the flow of ground water in a three-dimensional aquifer system. The numerical method used to solve the set of simultaneous equations is the strongly implicit procedure (SIP) described by Stone (1968) for two dimensional problems, extended to three dimensions by Weinstein, Stone, and Kwan (1969) and published as computer code documentation by Trescott (1975) and Trescott and Larson (1976). The modeler must evaluate the applicability of the program to the specific geohydrologic system that he intends to simulate. The aquifer system may be hetergeneous and anisotropic and may include both artesian and water-table conditions. The physics of the flow field may be approximated by either of two basic sets of simultaneous equations. For an aquifer system in which storage in confining beds is not significant, the flow field may be described by the equation for the flow of ground water in three dimensions. For an aquifer system in which storage in confining beds is significant, the flow field may be described by a sequence of equations for two dimensional horizontal flow in each of the more permeable beds coupled by the equations for one-dimensional vertical flow through the poorly permeable beds. The computer program allows both specified-head and specified-flux boundaries. At a specified-head boundary, the head is not allowed to change during any given interval of time. At a specified-flux boundary, the flow into or out of the aquifer system is not allowed to change during any given interval of time. The computer program also allows head-dependent source/sink boundaries at which the flow rate is calculated to represent recharge from or discharge to a stream. The program documented in this report differs from earlier computer codes in several major respects. The data array structure and data processing techniques used in this code allow the economical processing of extremely large grids, essentially unconstrained by the physical limits of a computers' memory. A layer swapping algorithm is utilized, based on the techniques developed for the "HULL" programs by Matuska, Durrett, etal (1973). The program is written completely in structured FORTRAN and uses the FLECS structured precompiler written by Beyer (1975). The output from FLECS includes standard FORTRAN IV for the CDC Cyber-176 computers. The FLECS code is listed in Appendix VII and the FORTRAN code is listed in Appendix VIII. The input/output logic within the program is totally different than the schemes used in Trescott and Larson (1976). The user selects run-time options using the free-format NAMELIST I/O. Data grid initialization uses "cube" input techniques, which drastically reduces the number and complexity of input records. Output "cubes" may be defined to reduce the volume of printed output and focus the output on areas of greatest interest to the modeler. Operation of the program from a remote terminal is thereby made more efficient. Extensive backup and restart logic has been designed into the code to allow the modeler to break up long or complex simulations into stages. Thus, the user may review results and modify the model parameters as desired. The I/O structure and program design ease the logistics of changing grid schemes and reinitializing the estimates of geohydrologic parameters compared to earlier simulation programs. The data processing techniques designed into the program are outlined in an effort to assist the user and to indicate potential problem areas if attempting to run this code at computer facilities that do not have CDC Cyber computers available. Detailed instructions on the input data structure, together with examples of input and output, should5 facilitate the actual use of this program. References to other publications are provided which discuss the physics of ground-water flow, the finite-difference solution of partial differential equations, structured programming languages and the specifics of Control Data Corporawion Cyber computer ccftvnrc. This report is the result of cooperation between the U.S. Geological Survey, the U.S. Bureau of Indian Affairs, and the Office of the New Mexico State Engineer. CAUTION TO USERS AND DISCLAIMER The user of this program is cautioned to verify that the program is in fact functioning as intended for the specific data being used. Not all options have been exhaustively tested; therefore, the program is subject to revision as any subsequent errors are encountered. The user is advised to contact the authors and obtain information as to revisions in the program which may have been made since this document was released. In spite of these cautions, the program is believed to be operational as described. The data processing techniques used in the program were developed on computers manufactured by Control Data Corporation. The use of the brand name in this report is for identification purposes only and does not imply endorsement by the U.S. Geological Survey. REPRESENTATION OF GEOHYDROLOGIC SYSTEMS If the geohydrologic system to be simulated is a simple aquifer, a two-dimensional model may be adequate. The three-dimensional model is capable of simulating the response of complex aquifer systems. An aquifer system is a heterogeneous body of intercalated permeable and less permeable material which acts as a hydraulic unit of regional extent. The computer program documented herein evolved from that of Trescott (1975) as modified by Trescott and Larson (1976). Representation of the aquifer system The digital simulation of geohydrologic systems requires that the continuous functions of space and time be made discrete. A rectangular block-centered grid (in which variable grid spacing is permitted) is used to form the finite difference approximations for the derivatives in the flow equations. Geohydrologic systems in which storage in confining beds is not significant If confining beds are thin relative to the vertical dimension of the cells of the model, each cell will represent several beds of both permeable and less permeable material. In the macroscopic scale of the model the cell is homogeneous although possibly anisotropic. The flow field may be described by the equation for ground-water flow in three dimensions. (A)— !-(K Hi) + J_(K J^)+-!(K i5. ) = S |£ + W (x,y,z,t) ax x ax' 3y y ay az z ^z' s 3t in which K , K , K are the hydraulic conductivities in the x, y, and x' y z ' J z directions (LT~~ ); h is the hydraulic head (L); S is the specific storage (L ); 6 -i W is the volumetric flux per unit volume (T ). This is equivalent to Trescott's equation 3 (Trescott, 1975, p. 3), and may be solved by the computer program. Alternately, equation A can be multiplied by the thickness of the hydraulic unit and expressed as: (B) JL(T ——.) + —— (T l2l) + b -MK Hi) = slA+ bW (x,y,z,t) v ' 9 x v x a x' 3 y y 3 y 9 z v z 3 z' 31 in which 2 -1 are the transmissivities in the x and y direction (L T ); _ i K is the hydraulic conductivity in the z direction (LT ); Z h is the hydraulic head (L); S is the storage coefficient (dimensionless); b is the thickness of the hydraulic unit (L); W is the volumetric flux per unit volume (T ). This is equivalent to Trescott's equation 4 (Trescott and Larson, 1976, p. XV) and may be solved by the computer program. The finite-difference approximation and the solution algorithm used by the computer program are described by Trescott (1975). Geohydrologic systems in which