Development of a Quasi Three-Dimensional Groundwater Model for a Portion of the Nevada Test Site

Page 1

DEVELOPMENT OF A QUASI THREE-DIMENSIONAL GROUNDWATER MODEL FOR A PORTION OF THE NEVADA TEST SITE Water Resources Center Desert Research Institute University of Nevada System Prepared for the U.S. Department of Energy, Nevada Operations Office, under contract EY-76-C-08-1253.

Page 2

This report was prepared as an account of work sponsored by the United States Government. Neither the United States nor the United States Energy Research and Development Administration, nor any of their employees, nor any of their contractors, subcontractors, or their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness or usefulness of any information, apparatus, product of process disclosed, or represents that its use would not infringe privately-owned rights. AVAILABLE FROM THE NATIONAL TECHNICAL INFORMATION SERVICE U.S. DEPARTMENT OF COMMERCE

Page 3

DEVELOPMENT OF A QUASI THREE-DIMENSIONAL GROUNDWATER MODEL FOR A PORTION OF THE NEVADA TEST SITE Water Resources Center Desert Research Institute University of Nevada System

Page 4

Invaluable assistance from several sources is gratefully acknowledged. For pro viding the original suggestion as to the nature of this study, the author wishes to thank Dr. Paul R. Fenske. Thanks are also due to Dr. Michael E. Campana for assisting with the computer work. The author wishes to express deep regret that his original thesis committee chairman, Dr. George B. Maxey, passed away before completion of this project.

Page 5

FINITE DIFFERENCE EQUATIONS.......................................................... 1 DEVELOPMENT OF A QUASI THREE-DIMENSIONAL GEOLOGY AND HYDROGEOLOGY FROM PAHUTE MESA TO YUCCA FLAT, NEVADA TESTSITE.............................................. 15 Regional Geologic Setting............................................................... 19 Rock Units with Associated Valuesof Hydraulic Conductivities 19 GROUNDWATER MODEL OF PAHUTEMESA TO YUCCA FLAT • • 32 APPENDIX A. HEAD VALUES FORTHEORETICAL BASIN .... 41 APPENDIX C. GROUNDWATER POTENTIALMODEL......................... 48

Page 6

Groundwater systems are often modeled on computers using a finite difference method in two dimensions, height and length. These models do not account for water movement in the lateral extent normal to the model. An empirical correction matrix has been devised that will partially adjust a two-dimensional model for flow in the third dimension. In an isotropic and homogeneous medium, the correction matrix averages a 4.9% difference in potential values. The correction matrix has been applied to the groundwater system from Pahute Mesa to Yucca Flat on the Nevada Test Site, a distance of 24 kilometers. The geology of the area is complex and contains folded and thrust faulted Paleozoic through Tertiary rocks. Groundwater is recharged at Pahute Mesa and flows to Yucca Flat at the rate of 64.8 acre feet/year. The maximum water flow under all of Yucca Flat and into the Ash Meadows groundwater system is 1240 acre feet year. 1.

Page 7

The purpose of this project is to better define and understand the groundwater system from southern Pahute Mesa to Yucca Flat, Nevada. The area studied lies entirely within the Nevada Test Site, located 90 kilometers north of Las Vegas, Nevada. Because of the radioactive contamination of groundwater at the Nevada Test Site, it is important to understand the regional groundwater movement. Funding for this project came from the University of Nevada, Desert Research Institute, Water Resources Center as a part of the Groundwater Study Project of the Department of Energy. All computer work was done at the University of Nevada, Reno, on a Control Data Corporation #6400. FINITE DIFFERENCE EQUATIONS Mathematical Basis For Finite Difference Equations Two types of equations are useful in groundwater modeling problems. First are continuity equations that express the conservation of mass. Second are momentum equations, such as Darcy's Law, that define flow through a saturated isotropic porous medium. These two types of equations are combined to yield a groundwater flow equation that is applicable to the region studied (Huntoon, 1972). The following is developed in the two dimensions Y, Z. GIVEN: Continuity equation for fluid in a porous medium. specific discharge in Z direction (LT 1) specific discharge in Y direction (LT 1)

Page 8

AND: Darcy's Law for flow through a porous medium Under isotropic conditions, Darcy's Law can be written in vector components Substituting the vector components of Darcy's Law into the equation of continuity Which is the equation governing a compressible fluid in an isotropic, heterogeneous porous medium. If we assume that the density is constant, the hydraulic conductivity is constant and steady state conditions prevail, then Eq. 4 reduces to This is the Laplace Equation for groundwater flow and is quite difficult to integrate except for the simplest boundary condition. To find the finite difference equation for an interior node, the second order partial derivatives of Laplace's Equation must be replaced. For any value of (Y,Z) the equation for <j) (Y,Z) can be expanded about any point (Yq,Zq) (Boas, 1966).

Page 9

next node in the positive direction on the Y axis next node in the positive direction on the Z axis point about which the series expansion is taken a small remainder representing the other terms in the series By use of a forward and a backward difference technique (Freeze, 1969, p. 45), which is too extensive to present here, the second order derivatives of the Laplace Equation are approximated by By substituting values for the derivatives in Eq. 5 and letting R equal one nodal space, the equation becomes

Page 10

node 1 = the first node along the forward distance along the Z axis used in Eq. node 2 = the first node along the backward distance along the Z axis used in node 3 = the first node along the forward distance along the Y axis used in Eq. node 4 = the first node along the backward distance along the Y axis used in which is the finite difference approximation of Laplace's Equation for any constant R.

Page 11

lor any interior node in an isotropic, homogeneous medium. Or, by other notations, Where (Y,Z) = the node where head is to be determined i = the number of iteration Eq. 11 is known as the Jacobi Method and is theoretically important. However, other techniques are better suited to computer means of matrix relaxation. This paper uses the Gauss-Seidel Method exclusively, which is about twice as fast in machine time as the Jacobi Equation. In the Gauss-Seidel Method the heads at nodes (Y-1,Z) and (Y,Z-1) have already been computed by iteration i+1 when node (Y,Z) is computed. The Boundary Conditions for Finite Difference Equations Along the water table, <j)= h, where h = the distance above the datum - usually mean sea level. At depth there is an assumed impermeable basement where 3<f>/9y = 0 and 34>/3z 0, which is a zero flux condition. In order to position a zero flux condition procedure is used to place an impermeable boundary along the sides of the matrix in positive Z direction. The location of the boundary is halfway between Z and Z+l. Where the impermeable basement, zero flux in the Y direction, meets the side of the matrix, zero flux in the Z direction, a special boundary condition exists. The conditions are met in the corners of the matrix and the equation ^(Y+l,Z+l) = 4> (Y,Z) applies. The location of this boundary is halfway between Z and Z+l, and is also halfway

Page 12

In an isotropic heterogeneous medium, different hydraulic conductivities are used within the matrix to model different rock types. When using computer models to solve for the head potentials, the absolute values of hydraulic conductivity are not needed as long as the ratios between the values are known. This does not apply to timedependent situations. Within any layer where K does not change, the Gauss-Seidel Model applies. The boundary between two different hydraulic conductivities is modeled in this paper by a harmonic mean technique which places the boundary between two To compute the head at (Y,Z) and (Y+1,Z) the harmonic mean equations are

Page 13

= hydraulic conductivity number 1 at point (Y,Z) = hydraulic conductivity number 2 at point (Y+1,Z) = Harmonic mean hydraulic conductivity which equals DEVELOPMENT OF A QUASI THREE-DIMENSIONAL CORRECTION MATRIX The implementation of finite difference equations is best suited to computer programs using recursive procedures. Finite difference equations are used with recursive procedures to allow the computer solution to converge on the exact solution, represented by the partial differential equation. The difference between the exact solution and the computer solution (convergence criteria) is ±0.000001 in this paper. All recursive programs have been converged to the same tolerances. Although groundwater hydraulics are three-dimensional, they are usually solved by two-dimensional finite difference equations. This is done for several reasons. First, the amount of computer memory needed to solve a three-dimensional problem is much larger than for a two-dimensional problem, although this is less of a problem today than when computers were first used to solve finite difference equations in the 1960’s. mathematical models should be limited to small, simple basins that can be adequately represented by 7500 nodes. The best way to handle large, complex basins is with a series of two-dimensional models representing vertical sections taken parallel to the direction of dip to the water table slope” (Freeze, 1969). Today (1978) computer memory capable of 50,000 nodes is commonplace and will cover most modeling situations.

Page 14

The second reason involves computer time. As the number of nodes increases, the number of computations for a three-dimensional model increases exponentially over the number of steps in a simpler two-dimensional model. It can be seen in Figure 3 that the extra dimension increases the number of computations at a rapid rate. nodes along a unit axis. Figure 3 The machine time needed to compute the extra dimension can turn computer time into a major expense. Third, the number and complexity of the finite difference equations increases significantly when a three-dimensional solution is used. The three-dimensional two layered model used in this paper requires a 350 percent increase in program length For the reasons noted above, the users of two-dimensional models define the basin to be modeled as infinite in the lateral extent. In most cases this assumption allows a simple and fast model that gives acceptable results. However, the assumption carries with it the supposition that the flow paths of the groundwater are parallel in the lateral extent. Under these conditions, areas of recharge or discharge cannot be The major advantage of three-dimensional models, besides the added geological input, is that they allow water to diverge from a recharge area and to converge at a discharge area. This gives the best possible solution to hydraulic groundwater problems. But how much error is introduced by the use of a two-dimensional model that does not account for water movement in a lateral direction, and can this be partially corrected? To examine this problem, a computer program was written to emperically compare the results of a three-dimensional and a two-dimensional model.

Page 15

An imaginary basin was envisioned, with a single recharge area and a single discharge area. The water table was linear between the two areas. The values of head chosen for the high and low areas were arbitrarily set at 20 and 10, respectively; steady state conditions were assumed. The aquifer was assumed to be isotropic and the coding allowed for both homogeneous and heterogeneous cases. The development behind the computer program was first aimed at presenting the most basic case possible, with variations of hydraulic conductivity to be added later. The imaginary basin was modeled by a three-dimensional model called Box and a two-dimensional model known as Mat. Both models were designed so that two different hydraulic conductivities could be examined, an upper horizontal layer of rows Y = 1 to 6.5, and a lower horizontal layer of rows Y = 6.5 to 12. The model Box (see Figure 4) is a (12,21,20) matrix consisting of 5040 nodes. The model Mat (see Figure 5) is a (12,20) matrix having 240 notes. In Figure 6 the space relationship between Box and Mat can be seen. The shaded area is the matrix Mat, centered along the X axis inside the matrix Box. The given water table is initialized along line R-D. The same criteria was used in both models for depth and length of basin, position and value of the water table, and hydraulic conductivity. Since all parameters are the same, except that the model Box also computes head along the X axis, it can be concluded that any difference in the results are due to water movement in the X direction. This water movement is the expansion of water moving away from the recharge area and the contraction of water approaching the discharge area. The flowline divergence and convergence along the flow paths at the water table is shown by the dotted lines in Figure 6. Homogeneous conditions were assumed and both matrices were relaxed by computer. The results of Mat and the area of Box that corresponds to Mat (see Figure 6) were compared. Both matrices have the same mean value of The range of the computed values is greater in the two-dimensional solution (see Appendix A, Tables 2 and 3). Head values in the recharge area of Mat are too large because there is no divergence of flow in the X direction. Similarly, in the discharge area of Mat there is no convergence of flow in the X direction and the head values are too small. Since the center of the two-dimensional solution lies midway between values that are too large and those that are too small, it is the most accurate

Page 17

(see Appendix B, Table 4). The water table values have a correction factor of zero because the water table is the same in both models. The zero values are excluded from all statistical analysis of the matrix. In order to observe and correct for the differences between the two solutions, a new matrix called Corrmat (for correction matrix) was devised. Both the threeFigure 6 Spacial Relationship of Box and Mat dimensional and two-dimensional solutions were divided by their mean value of 15.0. This presents each solution as a field of percentages about the mean. Since the values of head in the two-dimensional solution had a greater range, the percentiles also had a greater range. Then, by simply subtracting the percentiles at each node of Mat from the corresponding node of Box, one produces a matrix (12,20) that shows the degree of difference at each point. The Corrmat matrix can be used to correct for flow in the X direction when applied to a two-dimensional solution of a three-dimensional system. To use this correction method:

Page 18

1. Relax the matrix to be corrected by finite difference means; 2. Either expand or contract the correction matrix by linear methods along the Y and Z axis to fit the matrix to be corrected; 3. Multiply the value at each node in the correction matrix by the mean value of head within the matrix to be corrected; 4. Add the values from step 1 to each corresponding node or group of nodes in matrix to be corrected. To examine the nature of the correction factor, two different example matrices were computed. A "normal" case of isotropic homogeneous conditions were also computed to see the effect in the ideal case (see …