I'neef »a The Library dliit* Ontversity 8t Urlnra' Ifnoair.n Distance East (feet) ILLINOIS STATE WATER SURVEY CHAMPAIGN each LosrBooKi^^S" P^e ,o, its’^reKo S^?brfr^fror®l^^' BULLETIN 69 GWPATH: Interactive Ground-Water Flow Path Analysis by JOHN M. SHAFER Title: GWPATH: Interactive Ground-Water Flow Path Analysis. Abstract: GWPATH is an interactive software package for estimating horizontal fluid pathlines and travel times in fully saturated ground-water domains. GWPATH was developed for the IBM PC-AT, and compatibles, microcomputing environment and takes full advantage of interactive graphical display. The code features interactive data entry, forward and reverse pathline tracking, time-related capture zone analysis, multiple pathline capture detection mechanisms, multiple configurations for pathline starting locations, and variable time stepping. GWPATH is applicable to inhomogeneous, anisotropic complex flow domains. Data requirements include definition of the ground-water flow domain configuration, a hydraulic head distribution, hydraulic conductivity, and effective porosity. Both hardcopy printable output and six-pen plotter output are generated by GWPATH. Reference: Shafer, John M. GWPATH: Interactive Ground-Water Flow Path Analysis. Illinois State Water Survey, Champaign, Bulletin 69, 1 987. Indexing Terms: Computer codes, computer programs, contaminant source identification, ground water, ground-water models, microcomputers, nonuniform flow, numerical integration, pathline analysis, reverse tracking, time-related capture zones. STATE OF ILLINOIS HON. JAMES R. THOMPSON. Governor DEPARTMENT OF ENERGY AND NATURAL RESOURCES DON ETCHISON, Ph.D., Director BOARD OF NATURAL RESOURCES AND CONSERVATION Don Etchison, Ph.D., Chairman Walter E. Hanson, M.S., Engineering Seymour O. Schlanger, Ph.D., Geology H. S. Gutowsky, Ph.D., Chemistry Robert L. Metcalf, Ph.D., Biology Judith Liebman, Ph.D. University of Illinois John C. Guyon, Ph.D. Southern Illinois University Funds derived from University of Illinois administered grants and contracts were used to produce this report. STATE WATER SURVEY DIVISION RICHARD G. SEMONIN, B.S., Chief CHAMPAIGN 1987 1‘rinted by authority of the State of Illinois (7 87-500) TABLE OF CONTENTS Page DISCLAIMER.iv ABSTRACT.v TABLE OF NOTATION.vi PROGRAM SYNOPSIS.1 METHODOLOGY.2 Overview.2 Theoretical Background.3 Numerical Implementation.4 Assumptions/Limitations.13 References.14 INPUT REQUIREMENTS.15 Ground-Water Flow Domain Definition.16 Ground-Water Flow Parameters.17 Output Specifications.20 GWPATH Execution Control.21 GWPATH Dimensionality.23 OUTPUT.24 GWPATH VALIDATION.25 GWPATH DEMONSTRATION.31 Example No. 1: Forward Tracking.36 Example No. 2: Reverse Tracking.37 Example No. 3: Capture Zone Analysis.37 APPENDIX: Software Installation Instructions 41 DISCLAIMER Neither the developer nor any person or organization acting on behalf of him: a) makes any warranty, express or implied, with respect to this software; or b) assumes any liabilities with respect to the use, or misuse, of this software, or the interpretation, or misinterpretation, of any results obtained from this software, or for damages resulting from the use of this software. ABSTRACT GWPATH is an interactive software package for estimating horizontal fluid pathlines and travel times in fully saturated ground- water flow domains. GWPATH was developed for the IBM PG-AT, and compatibles, microcomputing environment and takes full advantage of interactive graphical display. The code features interactive data entry, forward and reverse pathline tracking, time-related capture zone analysis, multiple pathline capture detection mechanisms, multiple configurations for pathline starting locations, and variable time stepping. GWPATH is applicable to inhomogeneous, anisotropic complex flow domains. Data requirements include definition of the ground-water flow domain configuration, a hydraulic head distribution, hydraulic conductivity, and effective porosity. Both hardcopy printable output and six-pen plotter output are generated by GWPATH. TABLE OF NOTATION 1. Mathematical Symbols h q K f(x) g(x,y) da/dx da/dx grad h A nonboldface lower case letter: a scalar. A boldface lower case letter: a vector in two dimensions with components and qy. A boldface upper case letter: a tensor. A function of one variable. A function of two variables. The first order derivative of a = f(x) with respect to x. The first order partial derivative of a = g(x,y) with respect to X . (5h/3x, dh/dy): the gradient of scalar h. 2. Variables h K q s(t) V x,y Total hydraulic head (L). Hydraulic conductivity tensor oriented along principal horizontal axes of anisotropy (L/T). Effective porosity (dimensionless). Darcy velocity vector with two-dimensional horizontal components qjr and qy (L/T) . Ground-water pathline trajectory: s(t) = [x(t),y(t)] (L). Average pore ground-water velocity with two-dimensional horizontal components v^ and Vy (L/T). Horizontal position coordinates (L). PROGRAM SYNOPSIS Program Name: GWPATH, Version 3.0 Purpose: Estimate two-dimensional horizontal ground-water flow paths in complex flow domains. Track flow paths down-gradient from source to sink; or up-gradient from sink to source. Compute ground-water travel time (days) for each pathline. Estimate time-related capture zones for hydraulic sinks. Approach: Lagrangian Language: IBM Professional FORTRAN, Version 1.14 Mode: Interactive Recommended Hardware : IBM AT w/ hard disk / math co-processor EGA / high resolution color monitor HP-7475A or IBM-7372 six-pen plotter Software Requirements: DOS 2.x or greater Input Requirements: Flow domain configuration (i.e., grid structure) Hydraulic head distribution Hydraulic conductivity Effective porosity Source/sink locations (optional) Output: Screen graphics display Disk file of pathline(s) trajectory and travel time Optional plot (HP-7475A/IBM-7372 six-pen plotter) Assumptions: Regular, nonuniform, mesh-centered rectangular grid Steady state hydraulic head distribution Two-dimensional regional flow Features: Interactive data entry Inhomogeneous, anisotropic flow domain Continuous velocity distribution Fourth-order Runge-Kutta numerical integration Variable time stepping Multiple configurations for pathline starting locations Forward tracking; reverse tracking Multiple pathline capture detection mechanisms Time-related capture zone analysis Note: GWPATH supports only IBM standard devices. 1 METHODOLOGY The following sections describe the methodology implemented in GWPATH for estimation of two-dimensional, horizontal ground-water flow paths. For more in-depth discussions of ground-water hydraulics and/or ground-water flow modeling, several informative general references are available including: Bear, J., 1979. Hydraulics of Groundwater, McGraw-Hill Book Co., New York, New York, 569 pp. de Marsily, G., 1986. Quantitative Hydrogeology, Academic Press, Inc., Orlando, Florida, 440 pp. Mercer, J.W. and C.R. Faust, 1981. Ground-Water Modeling, National Water Well Association, Worthington, Ohio, 60 pp. Wang, H.F. and M.P. Anderson, 1982. Introduction to Ground- water Modeling, W.H. Freeman and Co., San Francisco, California, 237 pp. Overview There are four basic types of ground-water models, each addressing different physical phenomena; i.e., ground-water flow, solute transport, heat transport, and structural deformation (Mercer and Faust, 1981). Ground-water pathline analysis is an extension of ground-water flow modeling. It seeks a balance between predictive sophistication (i.e., solute transport modeling) and the amount and quality of available data. Solute concentrations are not directly modeled; however, center-of-mass flow path trajectories and ground-water travel times are modeled. There are a variety of regional and sub-regional ground-water contamination problems where flow path determination (especially for contaminant source identification) is a logical first step of analysis (Nelson and Schur, 1980). Further, pathline/travel time analyses, used as precursors to more elaborate evaluations, can aid in optimizing data collection efforts. The reverse pathline calculation capability of GWPATH provides the capability for estimating time-related capture zones of hydraulic sinks throughout the flow domain. A large number (i.e., usually between 100 and 300) of reverse calculated pathlines, all with the same travel time, can be used to approximate the envelope of influence of hydraulic sinks over the specified time period. For a detailed discussion of reverse pathline calculation of time-related capture zones, refer to Shafer (1987) . 2 Theoretical Background The ground-water flow system kinematics provide the geometry of saturated ground-water movement. A pathline is the projected two-dimensional representation of the mean trajectory and travel time an individual fluid particle* follows through the flow domain. The pathline describes the expected planar location and time of travel of the fluid from any initial location, along its trajectory, to its terminus (e.g., well or flow domain boundary). Trajectories, or pathlines, can be tracked forward (i.e., down-gradient) or in reverse (i.e., up-gradient) through the flow system. Multiple pathlines can be used to determine the envelope of influence of sinks within the flow domain; i.e., the time-related capture zone. The basic solution approach used to calculate pathlines is La- grangian. An infinitesimal mass of fluid (the particle) is followed as it would move through the horizontal, two-dimensional representation of the flow domain. The equations for fluid pathlines are derived from Darcy's law for incompressible fluids (see de Marsily, 1986), which is q = -K grad h (1) where q = Darcy velocity vector with two-dimensional horizontal components q^ and qy (L/T), K = hydraulic conductivity tensor with two-dimensional components and Ky, assuming the coordinate system is aligned with the principal horizontal axes of anisotropy (L/T), h = total hydraulic head, h = h(x,y) (L), and x,y = two-dimensional horizontal Cartesian coordinates (L). The average pore ground-water velocity is where V = 3_ n e V = average pore velocity vector with two-dimensional horizontal components v^, and Vy (L/T) , and ng = effective porosity (dimensionless). ( 2 ) Combining equations 1 and 2 results in v grad h (3) This "particle" has no physical meaning and mathematical representation for an infinitesimal mass of is used fluid. as a 3 with individual horizontal velocity components f(x,y) = = Vx (4.1) and , g(x,y) = ^ = ^y (^-2) A pathline, s(t) = [x(t),y(t)], in the flow domain is one for which the velocity, v, associated with the path satisfies equation 3 everywhere along the path. The determination of pathlines is achieved by simultaneously solving the initial value problems expressed by the set of ordinary differential equations (equations 4.1 and 4.2), which represent an equivalent form of equation 3. Initial conditions are x(0) = Xq y(0) = yo- Numerical Implementation The solution for two-dimensional ground-water fluid pathlines is accomplished by numerically implementing the approach outlined in the previous section. There are three basic components to the numerical determination of pathlines: 1) calculation of Lnternodal average pore ground-water velocities, 2) interpolation of intracell pore ground- water velocities, and 3) numerical integration of equations 4.1 and 4.2. The methods used for each component are discussed in the following paragraphs. GWPATH is capable of solving for pathlines in inhomogeneous, anisotropic flow domains. Consequently, two sets of internodal average pore ground-water velocities are calculated. The first set is velocities between nodes that are horizontally adjacent (i.e., west to east in the x-direction) throughout the flow domain. The second set of velocities is between nodes that are vertically adjacent (i.e., north to south in the y-direction). If the forward (i.e., down-gradient) mode is selected, all internodal velocities are multiplied by -1.0. Conversely, if reverse (i.e., up-gradient) tracking is desired, all internodal velocities are multiplied by +1.0. Figure 1 shows the mesh-centered flow domain grid configuration with orientation of flow parameters. The hydraulic conductivities input to CWPATH are assumed to be pointwise directional, measured at the nodes. An arithmetic mean is computed to represent the internodal hydraulic conductivity. Therefore, two sets of hydraulic conductivities are used; i.e., one set in the x-direction and one set in the y-direction. Effective porosities are assvuned to be representative of the porosity at each node and are scalar quantities. Internodal Velocity Calculation. The method used to calculate the internodal pore ground-water velocities in the x-direction is 4 y I i-i i+i ◄-DELX(I)-^ nj) +' > UJ rv i. I ■hdj) Kx(lJ) Ky(lJ) ne(IJ) Figure 1. Flow domain grid configuration and orientation of ground- water flow parameters where VELX(i,j) Kx(i-.j) 5h(i,j)' .ne(i.j). . 5x(i) ^ (5) Kx(iJ) = [Kx(i+l,j) + Kxd.j)] / 2 ne(i,j) “ (ne(i+l,j) + ne(i,j)] / 2 3h(i,j) » [h(i+l,j) - h(i,j)] ax(i) = DELX(i) 5 Similarly, in the y-direction internodal velocities are calculated as where VELY(i,j) Ky(i,j)^ '3h(i,j)' ne(i.j). . 3y(j) ^ ( 6 ) Ky(i,j) = [Ky(i,j+1) + Ky(i,j)] / 2 ne(i-.j) = [ne(i,j+l) + ne(i,j)] / 2 ah(i,j) = [h(i,j+l) - h(i,j)] ay(i) = DELY(j) The internodal velocities are calculated once and saved. They are the basis for interpolation of intracell velocities throughout the entire flow domain. Intracell Velocity Interpolation. The method of interpolating intracell velocity components at any (x,y) location within the flow domain is based on an approach developed by Prickett, Naymik, and Lonnquist (1981). It is an elaborate scheme created to maintain continuity of velocity across cell boundaries. A continuous velocity field is important in assuring high accuracy, especially for divergent flow which occurs, for example, away from the center of an injection well. Prickett, Naymik, and Lonnquist (1981) developed a three-step procedure to interpolate velocities using linear basis functions (i.e.. Chapeau functions) borrowed from finite element theory. The first interpolation uses the four immediate internodal velocities surrounding the current (x,y) location. These are the outlined arrows in figure 2. The relative position within a grid cell is used to calculate weighted averages of x-direction and y-direction internodal velocities for the initial interpolated velocities, VX and VY, at (x,y). The second interpolation uses the four next closest internodal velocities (figure 3). For the demonstration presented in figure 3, these velocities are the outlined arrows. In exactly the same manner as the first interpolation, weighted averages, based on the (x,y) position, are used to calculate the two intermediate velocity components, vk and i/y. Finally, the third step combines the first and second interpolations (figure 4) . The a and b coefficients from the linear basis functions are used to combine VX and i/x (VY and uy) to estimate the intracell velocity components at position (x,y). With the x- direction for explanation, XZ represents a weighting, or proportioning, between VX and i/x. As the position x nears I, XZ increases to a maximum of 0.5 as b goes to zero at I. Hence, the interpolated velocity is continuous across I. The same procedure is repeated for the y-direction velocity component. This process results in a continuous velocity field across cell boundaries. Figure 5 shows the internodal velocity components used in the interpolation scheme based on location within the flow domain. 6 VELOCITY VECTORS BETWEEN DELX(I) COLUMNS I AND I+l >- CHAPEAU FUNCTIONS FOR VELOCITY VECTORS BETWEEN ROWS J AND J + 1 VX = aVELX(I,J + l) + (l-a)VELX(IJ) VY = bVELY(I+l,J) + (l-b)VELY(IJ) Figure 2. Initial interpolation of intracell pore ground-water velocities (after: Prickett, Naymik, and Lonnquist, 1981) 7 VELOCITY VECTORS BETWEEN COLUMNS I-l AND I DELX(I) aVELX(I-l,J + l) + (l-a)VELX(I-l,J) bVELY(I+lJ-l) + (l-b)VELY(IJ-l) I-l I I+l CHAPEAU FUNCTIONS FOR VELOCITY VECTORS BETWEEN ROWS J-1 AND J + 1 C’ 0 ux = uy = Figure 3. Intermediate interpolation of intracell pore ground-water velocities (after: Prickett, Naymik, and Lonnquist, 1981) 8 CHAPEAU FUNCTIONS FOR y DIRECTION VELOCITY VECTORS XZ 1-XZ CHAPEAU FUNCTIONS FOR X DIRECTION VELOCITY VECTORS XZ = |b-0.5| YZ = |a-0.5| FINAL VELOCITY COMPONENTS VX = (l-XZ)VX + XZ-ux VY = (l-YZ)VY + YZ'uy Figure 4. Final interpolation of intracell pore ground-water velocities (after; Prickett, Naymik, and Lonnquist, 1981) 9 >• X INTERIOR OF MODEL M I I + l J=1 SIDES OF MODEL 1=1 I=NC CORNERS OF MODEL Figure 5. Internodal velocity components used in GWPATH interpolation scheme (after: Prickett, Naymik, and Lonnquist, 1981) 10 Numerical Integration. The determination of fluid pathlines is achieved by numerically solving the linked system of equations (i.e., equations 4.1 and 4.2) using a fourth order Runge-Kutta numerical integration scheme. There are numerous references to this approach in the literature. A discussion specific to pathline analysis can be found in Daly and Morel-Seytoux (1980). A variably calculated time step, At (days), is used to calculate successive locations [x(n),y(n)], [x(n+l),y(n+ 1 )] according to x(n+l) = x(n) + l/ 6 [ap + 2 a 2 + 2 a 3 + ai^] (7.1) y(n+l) = y(n) + l/ 6 [b 3 + 2 b 2 + 2 b 3 + b 4 ] (7.2) where a^ = At X f [x(n),y(n) ] bx = At X g[x(n),y(n)] a 2 = At X f [x(n)+ax/2,y(n)+bx/2] b 2 = At X g[x(n)+ax/2,y(n)+bx/2] a3 = At X f [x(n)+a2/2 ,y (n)+b2/2 ] b 3 = At X g[x(n)+a2/2,y(n)+b2/2] a 4 = At X f[x(n)+a 3 ,y(n)+b 3 ] b 4 = At X g[x(n)+a 3 ,y(n)+b 3 ] The Runge-Kutta algorithm requires the starting location (xQ,yQ). Time Stepping. The time step. At, is allowed to vary according to user-defined criteria. At the beginning of each Runge-Kutta iteration the time step is calculated to be approximately equal to the time required for the fluid to linearly traverse the longest dimension of the current cell divided by the user-supplied number of moves per cell. The larger the number of moves, the more accurate the solution (i.e., the shorter the time step). However, the shorter the time step, the greater the computational time necessary to reach a solution. The user is also required to enter the maximum desired time step which is invoked if the variable time step exceeds the maximum allowable At. All travel times are in days. Pathline Capture. Pathlines are often captured by sinks (forward tracking) or sources (reverse tracking) within the flow domain. If numerical capture mechanisms were not included in the model formulation, pathlines would oscillate indefinitely around sinks (or sources). A dual approach is used in GWPATH to terminate pathlines at capture points with travel times accurate to At. First, point locations (i.e., x,y coordinates) of all identifiable sinks and sources should be determined for the flow domain. These locations are part of the optional input to GWPATH. At the beginning of each time step, GWPATH calculates the distances from the current pathline position to the location of each sink (or source). GWPATH then projects the travel distance of the pathline using the interpolated Vj^ and Vy velocity components associated with the current position. If the projected travel distance is less than the distance to all sinks (or sources), the trajectory is extended At with no capture. However, if the projected pathline travel distance 11 I I+l di > di ;Continue d2 < di ;Backup time step, reduce A t. Figure 6. Schematic representation of pathline capture detection procedure 12 is greater than one or more of the distances to each sink (or source) , the pathline is assumed to be captured by the sink (or source) closest to the current pathline location. GWPATH has an override feature to the pathline capture detection method described above, which is diagrammatically presented in figure 6. If during any particular time step the distance from the previous position [x(n-l),y(n-l)] to the next position [x(n+l),y(n+1)] is less than the distance from the previous position to the current position [x(n),y(n)] a flag is raised. This situation indicates that there has probably been a 90 degree or more reversal in the direction of flow. The dashed line in figure 6 is orthogonal to the current trajectory. Such an occurrence is indicative of trajectory oscillation around a sink (or source). In this event, the entire numerical integration procedure is backed up one time step and a new mode is entered. The time step is reduced by 10 percent and the trajectory is recalculated beginning at location [x(n-1),y(n-1)]. This process is repeated, each time reducing the time step by 10 percent, until either the flag is no longer raised or the time step has been reduced to less then one minute. If At has been reduced to less than one minute, the terminus of the pathline is assumed to be [x(n-l),y(n-l)]. Assumptions/Limitations There are several tacit assumptions associated with two-dimensional numerical pathline analysis. The ground-water fluid flow paths produced by GWPATH are projected average pathlines. In reality ground water flows in three dimensions. With two-dimensional horizontal pathline analysis the vertical flow component is assumed zero. This implies that along any vertical section of the flow domain all particles are moving horizontally in the same direction and at the same speed. As configured, GWPATH assumes that the hydraulic head distribution is steady throughout the time plane of pathline calculation. However, successive time planes can be averaged to approximate the impact of a time-varying hydraulic head distribution on flow path trajectories. In this sense, the resulting pathlines become an averaged estimate of trajectories and travel times, not a true approximation of a transient traj ectory. As previously discussed, GWPATH requires a regular, although not necessarily uniform, mesh-centered rectangular grid to operate. Large irregularly shaped flow domains can be subset into regular components with GWPATH implemented on the various domain subsets. An important consideration is that the accuracy of the calculated pathlines is a function of the quality and quantity of available data as well as the grid density and length of time step. 13 References Daly, C.J. and H.J. Morel-Seytoux, 1980. An Analytical/Numerical Procedure for Predicting Mass Transport in Groundwater, CEP80-81CJD-HJM7, Department of Civil Engineering, Colorado State University, Ft. Collins, Colorado. de Marsily, G., 1986. Quantitative Hydrogeology, Academic Press, Inc., Orlando, Florida. Mercer, J.W. and C.R. Faust, 1981. Ground-Water Modeling, National Water Well Association, Worthington, Ohio. Nelson, R.W. and J.A. Schur, 1980. PATHS Groundwater Hydrologic Model, PNL-3162, Battelle, Pacific Northwest Laboratory, Richland, Washington. Prickett, T.A., T.G. Naymik, and C.G. Lonnquist, 1981. A Random-Walk Solute Transport Model for Selected Groundwater Quality Evaluations, ISWS/BUL-65/81, Illinois State Water Survey, Champaign, Illinois. Shafer, J.M., 1987. "Reverse Pathline Calculation of Time-Related Capture Zones in Nonuniform Flow," Ground Water, 25(3), pp. 283-289. 14 INPUT REQUIREMENTS GWPATH is designed for interactive data input via query response. The code is organized so that default uniform conditions (i.e., homogeneous hydraulic conductivity, effective porosity, and flow domain discretization) can be input directly or variable conditions can be read from disk files. Hydraulic head distributions are input from ASCII files. Figure 7 shows the series of data input screens. Units of measure are user-defined as either U.S. Customary or SI. Time units are always days. The following consistencies must be observed. U.S. Customary units Time: days Flow domain definition: Hydraulic conductivity: Velocity: feet/day Pathline: feet SI units Time: days Flow domain definition: Hydraulic conductivity: Velocity: meters/day Pathline: meters feet feet/day meters meters/day 2-D HORIZONTAL GROUND-WATER PATHLINE ANALYSIS Version 3.0, 1987 DISCLAIMER Neither the developer nor any person acting on behalf of him: (a) makes any warranty, express or implied, with respect to this software; or (b) assumes any liabilities with respect to the use, or misuse, of this software, or the interpretation, or misinterpretation, of any results obtained from this software, or for damages resulting from the use of this software. Execution suspended : Press ENTER to Resume Execution. Figure 7a. Screen 1: Program title and disclaimer 15 There are five categories of required data input for GWPATH: 1. Title and units 2. Pathline domain definition 3. Flow parameters 4. Output file specifications 5. Program execution control information Ground-Water Flow Domain Definition The flow domain is always assumed oriented so that the positive x-direction is from west to east, and the positive y-direction is from north to south on a flat page. Therefore, the location of the minimum x [XMIN] , minimum y [YMIN] dimension is in the upper left corner. The user is asked to input XMIN and YMIN, the number of column nodes [NC], and the number of row nodes [NR] of the flow domain. The user is also prompted for x-direction spacing [DELX] and y-direction spacing [DELY]. The fields provided for response must be completely filled. If data are to be entered from files, pressing (instead of typing a numeric value) will reveal the filename input query. In this manner a nonuniform grid can be used to describe the ground-water flow GWPATH Version 3.0, 1987 Title For This Session [This is an Input Demonstration] Units of Length (Feet or Meters) [Feet ] PATHLINE DOMAIN DEFINITION West (XMIN) Boundary [0.000000] No. of Column Nodes (NC) [51] North (YMIN) Boundary [0.000000] No. of Row Nodes (NR) [51] GRID SPACING Enter Numeric Real Value If Uniform, Otherwise Press ENTER For Variable Spacing Input File Delta-X, East-West Spacing [160.00] Delta-Y, North-South Spacing [160.00] Figure 7b. Screen 2: Flow field definition input 16 domain. Variable x-direction and/or y-direction flow domain grid spacing may be input from ASCII files constructed for list-directed input according to the following example. Record No. 1 2 3 I NC-1 EOF Value DELX(l) DELK(2) DELX(3) DELX(I) DELX (NC-1) There must be NC-1 entries in the DELX file and NR-1 entries in the DELY file. Ground-Water Flow Parameters Hydraulic conductivity data and effective porosity data are entered in a manner identical to pathline domain grid definition input. GWPATH Version 3.0, 1987 FLOW PARAMETERS Enter Numeric Real Value If Uniform, Otherwise Press ENTER For Variable Parameter Input File Hydraulic Conductivity (feet/day) X-Direction [50.000] Y-Direction [ ] File Name [HYDY1.DAT ] Effective Porosity (%) [15] File Containing Hydraulic Head Distribution [GWHEAD1.DAT ] File Containing Source/Sink Information [SSI.DAT ] Press ENTER If None Figure 7c. Screen 3: Flow parameter input 17 If hydraulic conductivity is uniform in either the x- or y-direction the actual value is input (the field must be completely filled). For inhomogeneous conditions, pressing will reveal the filename query. The name of a disk file that contains the variable hydraulic conductivity data is entered. The same procedure is used for effective porosity. ASCII grid data matrices must be formatted using a ROW, BEGINNING COLUMN, ENDING GOLUMN arrangement with up to six grid-node values per 80 column record. In this manner, a typical record would have the matrix row number for the first entry followed by a beginning column number and then an ending column number. The remainder of the record would contain the grid matrix numbers corresponding to elements defined from (row, beginning column) to (row, ending column). An example of this file structure is provided on the following page. Before beginning the grid-data input, the file must contain two header records. These records can be blank but typically contain GWPATH Version 3.0, 1987 OUTPUT CONTROL INFORMATION Printable File of Results [EXAMPLE.DAT ] Press ENTER If Option Not Desired Plot Results With HP7475A/IBM7372 Color Plotter, (Y)es or (N)o? Figure 7d. Screen 4: Output control menu 18 Illinois State Water Survey 2204 Griffith Drive Champaign, Illinois 61820-7495 ORDER FORM To receive a S?", 1.2 Mb diskette containing GWPATH3.0, IBM device drivers, and example data sets, please send check (or money order) for $125.00 payable to University of Illinois, c/o Illinois State Water Survey to John M. Shafer at the above address. Name _ Organization _ Address _ City _ State_Zip Date Phone (Please allow 2 to 3 weeks for delivery) information regarding the flow domain. The following FORTRAN source code is representative of how the files are read by GWPATH. OPEN(IUNIT,FILE='filename.ext') READ(lUNIT,100,END=300) NC, NR, NULL, XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX 100 FORMAT(3I5,6F10.3) READ(lUNIT,100,END=300) NULL, NULL, NULL, ZMIN, ZMAX 200 READ(lUNIT,100,END=300) J, II, 12, (Z(J,I), I = 11,12) GO TO 200 300 CLOSE (lUNIT) An example ASCII data file with the proper structure for GWPATH input ; is shown using the hydraulic conductivity data listed below, 5 1 5 1 0 0.000 8000.000 0.000 8000.000 10.000 20.000 0 0 0 10.000 20.000 1 1 6 11.258 11.686 1 1 .887 12.025 12.221 12.513 1 7 1 2 12.792 13.010 13.168 13.260 13.280 13.227 1 1 3 1 8 13.104 12.923 12.700 12.450 12.184 11.913 1 1 9 24 11.645 11.386 11.142 10.917 10.713 10.534 1 25 30 10.379 10.251 10.149 10.074 10.024 10.000 1 31 36 10.000 10.024 10.070 10.136 10.221 10.323 1 37 42 10.440 10.569 10.708 10.853 11.000 11.148 1 43 48 11.291 11.425 11.548 11.653 11.737 11.796 1 49 5 1 11 .826 11.842 11.858 2 1 6 11.148 11.580 11.954 12.176 12.405 12.666 2 7 1 2 12.951 13.189 13.370 13.494 13.553 13.544 2 1 3 1 8 13.507 13.340 13.087 12.868 12.628 12.376 2 1 9 24 12.117 11.866 11.626 11.403 11.200 11.019 2 25 30 10.861 10.727 10.619 10.535 10.477 10.442 2 3 1 36 10.430 10.441 10.472 10.522 10.591 10.674 2 37 42 10.771 10.880 10.997 11.119 11.244 11.368 2 43 48 11.488 11.601 11.702 11.788 11.856 11.902 2 49 5 1 11.925 11.940 11.954 5 1 1 6 16.254 16.224 16.200 16.183 16.170 16.141 5 1 7 1 2 16.105 16.081 16.079 16.102 16.139 16.180 5 1 1 3 1 8 16.224 16.270 16.318 16.367 16.416 16.465 5 1 1 9 24 16.513 16.561 1 6.608 16.654 16.698 16.741 5 1 25 30 16.783 16.824 16.863 16.901 16.938 16.973 5 1 3 1 36 17.007 17.038 17.068 17.096 17.120 17.142 5 1 37 42 17.159 17.172 17.179 17.180 17.172 17.155 5 1 43 48 17.128 17.087 17.033 16.963 16.879 16.781 5 1 49 5 1 16.669 16.543 16.404 EOF 19 GWPATH also allows an optional ASCII data file containing the ground-water flow domain source/sink location information. This file must also be prepared for list-directed input. The first record of the source/sink file must contain the number of sources followed by the number of sinks. The real x-y coordinate locations of all sources are listed first followed by the real x-y coordinate locations of all known sinks. The following example shows how the ASCII source/sink data file must be organized. Record No. 1 2 3 NSOUR+1 NSOUR-l-2 NSOUR-l-3 NSOUR+NSINK-t-1 EOF Output Specifications The user is also the results of the pathline analysis. If hardcopy output is not required for any particular analysis session, pressing before any other response will eliminate hardcopy output. If a filename is provided, GWPATH writes the results (i.e., the pathline(s) trajectory and travel time) to this file which can be printed after the program ends. Value NSOUR NSINK SOURCE(l,x) SOURCE(l,y) SOURCE(2,x) SOURCE(2,y) SOURCE(NSOUR,x) SOURCE(NSOUR,y) SINK(l,x) SINK(l,y) SINK(2,x) SINK(2,y) SINK(NSINK,x) SINK(NSINK,y) prompted to name an ASCII disk file for storing The final GWPATH control information input requirement is a "Yes" or "No" response to whether or not a hardcopy plot (i.e., HP-7475A or IBM-7372 six-pen plotter output) is desired. If the response to this question is "Yes", the plot will mirror the screen graphics until program execution is halted. The following pen stall numbers correspond to the indicated components of a hardcopy plot. Stall No. 1 2 3 4 5 6 Component of Plot Flow domain border, legend box Text, marker symbol outline Flow domain grid Sources, forward tracked paths Sinks, reverse tracked paths Capture zone Recommended Pen 0.7mm Black 0.3mm Black 0.3mm Yellow 0.3mm Blue 0.3mm Red 0.7mm Blue 20 GWPATH Execution Control GWPATH execution control enables the user to interactively specify several options for running the program, especially the manner in which pathline starting locations are determined. At the beginning of the pathline calculation procedure, a graphical representation of the flow domain is displayed on the screen followed by a command line displayed at the bottom. Figure 8 shows the display screen arrangement of the flow domain, legend, time - stepping control information, and the command line. The user must respond via keyboard input to the various options available for program execution. The first option to choose is whether pathline analysis or capture zone analysis is to be performed. Pathline analysis allows for forward or reverse tracking, multiple starting locations, and display of pathlines. Capture zone analysis automatically invokes reverse tracking and suspends display of individual pathlines. However, capture zone analysis enables the feature that connects the end points of all pathlines and will cross-hatch the interior of the resulting polygon. If pathline analysis is the option selected, the user is asked if forward pathline tracking or reverse pathline tracking is to be performed. The next three questions relate to the time specifications for the analysis. The first question asks for the maximum allowable travel time, in days. After the maximum travel time is entered, the user is asked for the minimum number of moves per grid cell and the absolute maximum time step, in days. The previous discussion on time¬ stepping explains how these parameters affect the calculation of pathlines. Once the time - stepping specifications have been entered, the user is asked to choose the mode of pathline starting location data entry. The five different ways to enter pathline starting locations are: 1. Manual cursor entry (pathline analysis only) 2. Manual keyboard entry (pathline analysis only) 3. Circular starting configuration from cursor entry 4. Circular starting configuration from keyboard entry 5. ASCII file of starting locations (pathline analysis only) If cursor entry is selected, a cross hair appears on the screen in the middle of the flow domain. The arrow keys are used to move the cursor to the desired pathline starting location. The key toggles between large cursor increments and small cursor increments. At any time the real coordinate location can be determined by pressing . To select any particular location to start a pathline press . Multiple pathlines can be calculated by repeating this process a maximum of 299 times. To end selecting pathline starting locations and begin their calculation the user must press . To select precise starting locations, the user should use the keyboard input option. With this option the user types the exact 21 0 8000 LEGEND • Hydraulic Source —Foruard Tracked Path ■ Hydraulic Sink — Reverse Tracked Path Foruard Tracking Travel Tine = 10000 Moves per Cell =2 Tine Step = 100 No. of Path lines = 10 Units = Feet (C)apture Zone ftnalysisi (F)1 om Path ftnalysisi or (S)top? Figure 8. Example display screen for GWPATH (integer) x-y coordinates of pathline starting locations. In either case, once all starting locations have been entered the user is asked if the locations should be saved in an ASCII file. If the answer is "Yes" the user must enter a filename for the locations. Instead of entering each pathline starting location individually, pathline starting locations can be automatically calculated around a circle defined by the user. This option requires the user to enter the x-y coordinates of the center of a circle and then provide a radius from the center. GWPATH will then automatically calculate pathline starting locations at equal arc lengths around the circle. As is the case for manual starting location entry, the center of the circle can be input by using the cursor or by direct keyboard entry. 22 In addition to entering pathline starting locations either manually or with a circular configuration, GWPATH allows for starting location input from an ASCII file. If the ASCII disk file option is selected for input of pathline starting locations, the file must be formatted for list-directed input according to the following example. Record No. 1 2 Value BECIN(l,x) BECIN(l,y) BECIN(2,x) BECIN(2,y) NPATH BEGIN(NPATH,x) BEGIN(NPATH,y) EOF GWPATH Dimensionality GWPATH dimension limitations are as follows: • Maximum number • Maximum number • Maximum number • Maximum number • Maximum number • Maximum number of columns = 99 of rows = 99 of sources = 20 of sinks = 20 of pathlines per execution = 300 of discrete points per pathline = 1000 23 OUTPUT By entering a filename at the "Printable File of Results [ ]" prompt, an ASCII file of the pathline trajectories and travel times will be created. The file will automatically be saved in the current directory when execution of GWPATH is terminated. The file can then be printed in the same manner as any ASCII file. If is pressed at the prompt, the results of the current session will not be saved. following Point listing is an example of GWPATH file output DEMONSTRATION OF OUTPUT TO FILE Pathline No. = 1 Travel Time (days) = 500.00 No. of Points = 51 X ( feet ) Y ( feet ) Time (davs) 1 1943.00 1814.00 0.00 2 1946.24 1823.08 10.00 3 1949.45 1832.66 20.00 4 1952.53 1842.73 30.00 5 1955.35 1853.23 40.00 11 1961.26 1928.26 100.00 12 1960.82 1943.37 110.00 13 1960.66 1959.35 120.00 14 1960.77 1976.26 130.00 15 1961.17 1994.15 140.00 25 1979.11 2161.50 240.00 26 1981.47 2176.13 250.00 27 1983.87 2190.54 260.00 28 1986.29 2204.74 270.00 29 1988.75 2218.73 280.00 30 1991.23 2232.52 290.00 45 2030.26 2419.51 440.00 46 2032.92 2431.10 450.00 47 2035.60 2442.60 460.00 48 2038.30 2454.01 470.00 49 2041.01 2465.33 480.00 50 2043.74 2476.57 490.00 51 2046.47 2487.74 500.00 Ik GWPATH VALIDATION GWPATH was validated by forward tracking a pathline and comparing the estimated trajectory with a known calculated trajectory. The reverse tracking option was subsequently validated by reverse tracking the known pathline. Finally, the effective porosity was divided by 2.0, doubling the velocity. Therefore, it should take one-half the original travel time for the same trajectory, all other factors remaining the same. To accomplish the validation a regular 8000 ft by 8000 ft homogeneous, isotropic flow domain was constructed with a hydraulic head distribution as shown in figure 9. The validation flow domain is symmetric about the north-south centerline; i.e., x = 4000 ft. Consequently, pathlines along the centerline are linear. This characteristic is used to validate GWPATH. The following parameters describe the validation flow domain and the imposed aquifer properties. •No. of column nodes = 51 •No. of row nodes = 51 • Delta-X = 160 ft • Delta-Y = 160 ft • Kjj = 50 ft/day • Ky = 50 ft/day • ne = 0.5 The following time - stepping constraints were used for the code validation. • Maximum allowable travel time = 10,000 days • Minimum number of moves per grid cell = 1 • Maximum allowable time step = 25 days He.d - 150 ft 120 - 110 - 1 f 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 I |4U I : 1 - - 120 - 110 - 8000 ‘ ‘ ‘ ‘ ‘ ‘ ‘‘ ‘ ‘ ‘ ‘ ‘‘ ‘ ‘ ‘'''' ‘‘‘ 0 " 8000 Figure 9. Hydraulic head (feet above MSL) for GWPATH validation 25 The forward tracking validation pathline computed by GWPATH begins at location (4000 ft, 1120 ft). The terminating coordinate (i.e., s[10,000 days]) equals (4000 ft, 7352.87 ft). The pathline trajectory is presented in the following listing and is shown in figure 10. FORWARD TRACKING VALIDATION Pathline No. = 1 Travel Time (days) = 10000.00 No. of Points = 401 lint X ( feet ) Y ( feet ) Time (dav: 1 4000.00 1120.00 0.00 2 4000.00 1136.09 25.00 3 4000.00 1152.17 50.00 4 4000.00 1168.25 75.00 5 4000.00 1184.33 100.00 6 4000.00 1200.41 125.00 7 4000.00 1216.49 150.00 8 4000.00 1232.56 175.00 9 4000.00 1248.64 200.00 10 4000.00 1264.71 225.00 11 4000.00 1280.78 250.00 12 4000.00 1296.85 275.00 13 4000.00 1312.92 300.00 14 4000.00 1328.98 325.00 15 4000.00 1345.05 350.00 16 4000.00 1361.11 375.00 17 4000.00 1377.17 400.00 18 4000.00 1393.23 425.00 19 4000.00 1409.29 450.00 385 4000.00 7.110.79 9600.00 386 4000.00 7125.93 9625.00 387 4000.00 7141.07 9650.00 388 4000.00 7156.21 9675.00 389 4000.00 7171.35 9700.00 390 4000.00 7186.48 9725.00 391 4000.00 7201.61 9750.00 392 4000.00 7216.75 9775.00 393 4000.00 7231.88 9800.00 394 4000.00 7247.01 9825.00 395 4000.00 7262.13 9850.00 396 4000.00 7277.26 9875.00 397 4000.00 7292.38 9900.00 398 4000.00 7307.51 9925.00 399 4000.00 7322.63 9950.00 400 4000.00 7337.75 9975.00 401 4000.00 7352.87 10000.00 26 0 8000 o (fc.OOOft, 1120ft) LEGEND ^Hydraulic Source — Forward Tracked Path ■ Hydraulic Sink ---Reverse Tracked Path t- O o o CD (<.000ft, 7352.87ft) Units • Feet Figure 10. Forward tracked GWPATH validation pathline To validate the code, this numerically determined pathline was analytically derived. The average hydraulic gradient over the length of the pathline was calculated by taking the mean of the hydraulic gradient in the north-south direction between each pair of nodes in the flow domain grid between the starting point and the ending point of the pathline. The resulting average gradient is ^ = 0.0062339 With equation 3, the average pore ground-water velocity trajectory of the pathline can be calculated as X 0.0062339 = 0.62339 ft/day along the Multiplying v^ by 10,000 days results in a linear travel distance of 6,233.9 ft. Adding the travel distance to the starting location gives the pathline ending location (4000 ft, 7353.9 ft). Comparison of this value with the ending location of the validation pathline calculated by GWPATH shows that within roundoff error the code correctly estimates ground-water pathlines in the forward tracking mode. The error between the GWPATH numerically calculated pathline and the same pathline analytically determined is approximately 0.01%. After the forward pathline tracking methodology implemented in GWPATH had been validated, the reverse tracking was easily validated by simply reverse-tracking the forward tracked pathline. The same time stepping constraints used for the forward pathline validation were used for the reverse pathline validation. The pathline starting location for reverse pathline validation was the ending location of the forward tracked pathline; i.e., (4000 ft, 7352.87 ft). The total allowable travel time was 10,000 days. The following pathline trajectory is the result of the reverse tracking validation. The pathline for reverse calculation validation is presented in figure 11. 27 Point 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 REVERSE TRACKING VALIDATION Pathline No. = 1 Travel Time (days) = 10000.00 No. of Points = 401 X ( feet ) 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 Y ( feet ) 7352.87 7337.75 7322.63 7307.51 7292.39 7277.26 7262.14 7247.01 7231.88 7216.75 7201.62 7186.48 7171.35 7156.21 7141.08 7125.94 7110.80 7095.65 7080.51 Time (days) 0.00 25.00 50.00 75.00 100.00 125.00 150.00 175.00 200.00 225.00 250.00 275.00 300.00 325.00 350.00 375.00 400.00 425.00 450.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 4000.00 1377.18 1361.11 1345.05 1328.99 1312.92 1296.85 1280.78 1264.71 1248.64 1232.57 1216.49 1200.41 1184.33 1168.25 1152.17 1136.09 1120.00 9600.00 9625.00 9650.00 9675.00 9700.00 9725.00 9750.00 9775.00 9800.00 9825.00 9850.00 9875.00 9900.00 9925.00 9950.00 9975.00 10000.00 28 0 8000 Figure 11. Reverse tracked GWPATH validation pathline The results of the reverse tracking exercise clearly validate the reverse pathline analysis capability of GWPATH. The pathline ending location (in a backward sense) after 10,000 days is exactly the starting location for the same pathline forward tracked; i.e., (4000 ft, 1120 ft). The final code validation exercise was to decrease the effective porosity by a factor of 2, thereby doubling the average pore ground- water velocity. Instead of 0.50, ng equalled 0.25 for this case. The time stepping constraints remained the same as for the previous validation exercises. However, the total allowable travel time was also reduced by a factor of 2 to reconstruct the original pathline. The same forward tracking pathline starting location was used; i.e., (4000 ft, 1120 ft). The anticipated results were that the length of the pathline calculated by GWPATH should be identical to the previous pathlines because the mean velocity was doubled while the travel time was halved. The resulting pathline trajectory is as follows. 29 REDUCED EFFECTIVE POROSITY Pathline No. = 1 Travel Time (days) = 5000.00 No. of Points = 201 oint X ( feet ) Y ( feet ) Time (davs 1 4000.00 1120.00 0.00 2 4000.00 1152.17 25.00 3 4000.00 1184.33 50.00 4 4000.00 1216.49 75.00 5 4000.00 1248.64 100.00 6 4000.00 1280.78 125.00 7 4000.00 1312.92 150.00 8 4000.00 1345.05 175.00 9 4000.00 1377.17 200.00 10 4000.00 1409.29 225.00 191 4000.00 7050.21 4750.00 192 4000.00 7080.50 4775.00 193 4000.00 7110.79 4800.00 194 4000.00 7141.07 4825.00 195 4000.00 7171.35 4850.00 196 4000.00 7201.61 4875.00 197 4000.00 7231.88 4900.00 198 4000.00 7262.13 4925.00 199 4000.00 7292.38 4950.00 200 4000.00 7322.63 4975.00 201 4000.00 7352.87 5000.00 As expected, the pathline is identical to the previous pathline, only the travel time is one-half the original 10,000 days. 30 GWPATH DEMONSTRATION Several example applications of GWPATH are presented in order to demonstrate its use. For demonstration purposes, a complex ground- water flow domain was developed with the following specifications. •No. of column nodes = 51 •No. of row nodes = 51 • Delta-X = 160 ft • Delta-Y = 160 ft • Inhomogeneous, anisotropic aquifer conditions • Constant head boundaries • Steady state hydraulic head • Five hydraulic sources; five hydraulic sinks Figure 12 is a contour plot of the steady state hydraulic head distribution of the demonstration ground-water flow domain. Figure 13 shows the spatial variability of x-direction hydraulic conductivity Figure 12. Hydraulic head (feet above MSL) for GWPATH flow path analysis demonstration 31 which ranges from 10 ft/day to 20 ft/day. In similar fashion, figure 14 shows the y-direction spatially varying hydraulic conductivity which ranges from 10 ft/day to 200 ft/day. Finally, the effective porosity is assumed to range from a low of 0.15 to a high of 0.25 over the flow domain (figure 15). Figure 13. X-direction hydraulic conductivity (ft/day) for GWPATH flow path analysis demonstration 32 y > X Figure 14. Y-direction hydraulic conductivity (ft/day) for GWPATH flow path analysis demonstration 33 Figure 15. Effective porosity (dimensionless) for GWPATH flow path analysis demonstration 34 0000 are; The five hydraulic sources and their steady state source strengths X-Coordinate 1760 ft 3680 ft 1120 ft 1440 ft 2880 ft Y-Coordinate 1120 ft 1760 ft 2240 ft 3200 ft 3200 ft Injection Rate 200 gpm 400 gpm 600 gpm 500 gpm 500 gpm The five hydraulic sinks and their steady state source strengths are: X-Coordinate Y-Coordinate Pumping Rate 3040 ft 4640 ft 400 gpm 4320 ft 4320 ft 400 gpm 5440 ft 3680 ft 400 gpm 4640 ft 6560 ft 500 gpm 5760 ft 5600 ft 600 gpm Figure 16 shows the orientation of sources and sinks throughout the flow domain. LEGEND A Hydraulic Source — Forward Tracked Path ■ Hydraulic Sink -Reverse Tracked Path -► X y Units = Feet Figure 16. Source/sink locations for GWPATH flow path analysis demonstration 35 Example No. 1: Forward Tracking This example demonstrates the forward pathline tracking feature of GWPATH. Both the manual selection and circular configuration options for pathline starting locations were used to perform the pathline analysis shown in figure 17. The pathlines beginning in the upper righthand (i.e., northeast) boundary of the flow domain were individually selected using the cursor input option. The starting locations of the forward tracked pathlines emanating from the source located at (2880 ft, 3200 ft) were calculated using the circular starting location configuration option. Equal arc distances around an imaginary circle centered at (2880 ft, 3200 ft) with a 100 ft radius were used to automatically compute 20 pathline starting locations. The following time - stepping constraints were used for this example. • Maximum travel time = 5000 days • Minimum moves per cell = 2 • Maximum allowable time step = 25 days Comparison of the pathlines calculated by GWPATH (figure 17) with the hydraulic head distribution (figure 12) confirms the sensitivity of the methodology to complex flow domains. 0 0000 LEGEND A Hydraulic Source — Forward Tracked Path ■ Hydraulic Sink — Reverse Tracked Path Units = Feet Figure 17. GWPATH forward pathline tracking demonstration 36 Example No. 2: Reverse Tracking To demonstrate the reverse pathline tracking option of GWPATH, a circular pathline starting configuration around a hydraulic sink was employed. In a backward sense, however, the input pathline starting locations are, in reality, the end points of the pathlines in the direction of flow. GWPATH simply tracks these pathlines up-gradient (i.e., backward) as opposed to forward, or down-gradient, through the flow domain. Gonsequently, instead of being convergent toward a hydraulic sink, the reverse pathline tracking flow paths are divergent away from a sink. The same flow domain used for Example No. 1 (figures 12 - 16) was used for the reverse tracking example application. Twenty equally spaced pathline starting locations were automatically calculated by GWPATH along an imaginary circle of 160 ft radius around the sink located at (4640 ft, 6560 ft). The following time - stepping constraints were imposed on the analysis. • Maximum travel time = 10,000 days • Minimum moves per cell = 2 • Maximum allowable time step = 25 days The results of the reverse tracking demonstration are presented in figure 18. The ability of GWPATH to adequately evaluate complex flow domains is well represented in this figure. The pathlines calculated by GWPATH and displayed in figure 18 are not necessarily obvious upon visual inspection of the hydraulic head contour map presented in figure 12. However, the reverse pathline analysis clearly shows the hydraulic relationship of the sink at (4640 ft, 6560 ft) to the hydraulic sources in the flow domain and to the other sinks. Further, given that by definition no flow can be orthogonal to a pathline, the contribution from various sources of flow can be estimated for any given sink from the distribution of reverse pathlines emanating from the sink. If, at a hydraulic sink, reverse tracked pathlines are started along even arc distances around the sink (as was the case for this example) , there is equal flow contribution in each flow tube represented by any adjacent pair of pathlines. Example No. 3: Capture Zone Analysis In addition to forward and reverse pathline tracking, an extension of reverse tracking is included in GWPATH which allows the calculation of time-related capture zones in ground-water flow domains. For a complete discussion of reverse pathline calculation of time-related capture zones using GWPATH, refer to Shafer (1987). A t-year capture zone in a three-dimensional flow regime is the bounded volume described by the set of all flow paths that intercept a hydraulic sink in t years. For the case of two-dimensional horizontal flow the vertical flow component is ignored and the capture zone becomes the planar area described by the set of all horizontal pathlines with t- year travel times to the sink. 37 0 8000 LEGEND A Hydraulic Source — Forward Tracked Path ■ Hydraulic Sink — Reverse Tracked Path Units = Feet Figure 18. GWPATH reverse pathline tracking demonstration GWPATH estimates the time-related capture zone by reverse tracking a large number of pathlines (typically between 100 and 300) all with the same travel time. GWPATH connects the resulting end points to form a polygon around the sink which represents the t-year capture zone. GWPATH then fills the polygon with a hatch pattern to display the enclosed area. The flow domain represented by the steady state hydraulic head distribution shown in figure 19 was used by GWPATH to calculate a time- related capture zone for demonstration. The flow domain was assumed to be inhomogeneous and isotropic with block hydraulic conductivities shown in figure 20. Effective porosity was assumed to be a constant 0.25. Four steady strength hydraulic sinks were located in the flow domain at the following coordinate locations. X-Coordinate Y-Coordinate Pumpine Rate 2400 ft 1600 6400 ft 2880 3840 ft 4640 2880 ft 6720 ft 200 gpm ft 400 gpm ft 500 gpm ft 300 gpm The 10-year (i.e., 3652 day) capture zone was calculated for the hydraulic sink located at (3840 ft, 4640 ft). The following time stepping constraints were used for the analysis. 38 Figure 19. Hydraulic head (feet above MSL) for GWPATH capture zone analysis demonstration • Maximum travel time = 3652 days • Minimum moves per cell = 2 • Maximum allowable time step = 20 days For this demonstration, 300 reverse tracked pathlines were used to delimit the capture zone. Figure 21 shows the resulting 10-year capture zone calculated by GWPATH around the sink at (3840 ft, 4640 ft). Although not performed for this example, once the capture zone has been calculated, the flow path analysis mode can be re-entered and a limited number of reverse pathlines (e.g., 20) overplotted on the capture zone for verification purposes or to estimate flow contributions from specific points in the flow domain. 39 8000 0 8000 Figure 20. Block hydraulic conductivities (ft/day) for GWPATH capture zone analysis demonstration 0 or 8000 rifT'irrrrTii i iiji i itiiiii ii iiirfrrr t t LEGEND A Hydraulic Source — Forward Tracked Path ■ Hydraulic Sink — Reverse Tracked Patii Units = Feet Figure 21. 10-year capture zone for hydraulic sink located at (3840 ft, 4640 ft) for GWPATH demonstration 40 APPENDIX Software Installation Instructions 41 Installing GWPATH 1. Establish directory on hard disk for GWPATH. mkdir c:\gwpath 2. Change directory. chdir c:\gwpath 3. Place diskette in drive A. 4. Copy diskette to hard disk. copy a:\*.'^ /v Running GWPATH To run GWPATH type GWPATH NOTE: All data files must reside in the current directory from which the call to run GWPATH is made. 42