ML20148P071
| ML20148P071 | |
| Person / Time | |
|---|---|
| Issue date: | 03/31/1988 |
| From: | NRC OFFICE OF NUCLEAR REGULATORY RESEARCH (RES) |
| To: | |
| References | |
| NUREG-1249, NUREG-1249-V01, NUREG-1249-V1, NUDOCS 8804080227 | |
| Download: ML20148P071 (184) | |
Text
{{#Wiki_filter:. . . - ._. .. NUREG-1249 Vol.1 l i i I t NRC Model Simulations in Support of the Hydrologic Code 4 Intercomparison Study (HYDROCOl\l) 4 f Level 1-Code Verification , I i s f U.S. Nuclear Regulatory ! Commission
- Office of Nuclear Regulatory Research on anc o
_e 0 , l ['Y 8 . , c ff ,'g
;. p. e f $[*~
l 5~*- * \ l l l r f 1249 R p g)g
NOTICE Availability of Reference Materials Cited in NRC Publicctions Most documents cited in N RC publications will be available from one of the following sources:
- 1. The NRC Public Document Room,1717 H Street, N.W.
Washington, DC 20555
- 2. The Superintendunt of Documents, U.S. Government Printing Office, Post Office Box 37082, Washington, DC 20013 7082
- 3. The National Technical Information Service, Springfield, VA 22161 Although the listing that follows represents the majority of documents cited in NRC publications, it is not intended to be exhaustive.
Referenced documents available for inspection and copying for a fee from the NRC Public Docu-ment Room include NRC correspondence and internal NRC memoranda; NRC Office of lospection and Enforcement bulletins, circulars, information notices, inspection and investigation notices; Licensee Event Reports; vendor reports and correspondence; Commission papers; and applicant and licensee documents and correspondence. The following documents in the NUREG series are available for purchase from the GPO Sales Program: formal NRC staff and contractor reports, NRC-sponsored conference proceedings, and NRC booklets and brochures. Also available are Regulatory Guides, NRC regulations in the Code of Federal Regulations, and Nuclear Regulatory Commission Issuances. Documents available from the National Technical Infonnation Service include NUREG series reports and technical reports prepared by other fedetal agencies and reports prepared by the Atomic
- Energy Commission, forerunner agency to the Nuclear Regulatory Commission.
Documents available from public and special technical libraries include all open literature items, such as books, journal and periodical articles, and transactions. Federal Register notices, federal and state legislation, and congressional reports can usually be obtained from these libraries. Documents such as theses, dissertations, foreign reports and translations, and non NRC conference p oceedings are available for purchase from the organization sponsoring the publication cited. Sing'e copies of NRC draf t reports are available free to the extent of supply, upon written request to the Division of Information Support Services, Distribution Section, U.S. Nuclear Regulatory Commission. Washington, DC 20S55. Copies of hdustry codes and standards ussd in a substantWe manner in the NRC regulatory process are maintained at the NRC Library, 7920 Norfolk Avenue, Bethesda, Maryland, and are available there for referer'ce use by the public. Codes and standards are usually copyrighted end may be purchased from the originating organization or, if they are American National Standards. from the American National Standards Institute,1430 Broadway, New York, NY 10018.
NUREG-1249 Vol.1
~
- :: T :z :_-_ ___._ _ _ _ _ _ __-- _ ___ _ - _ _ =
NRC Model S!mulations in Support of the Hydrologic Code Intercomparison Study (HYDROCOIN? Level 1-Code Verification Manuscript Completed: February 1988 Date Published: March 1988 Division of Engineering Office of Nuclear Regulatory Research U.S. Nuclear Regulatory Commission W2shington, DC 20555
,j "", ,, '
[ i, i e ,
l
)
i ABSTRACT l HYDROCOIN, an international study for examining ground-water flow modeling strategies and their influence on safety assessments of geologic repositories for nuclear waste, has been divided into three levels of effort. Level 1 evaluates the accuracy of the computer codes by comparing code results to analytical solutions and intercomparison of code results. Level 2 tests the codes' capac411 ties for simulating specific laboratory or field experiments. 1 Level 3 addresses the uncertainty of model results and the sensitivity of those results to changes in either the modeling approach or input parameters. A greater variety of analytic approaches and computer codes was compared in a systematic fashion than would have been possible because of the large number of parti-ipants involved and their corporate efforts. This report summarizes only the combined NRC project teams' simulation efforts on the Level 1 bench-marking problems. The codes used to simulate these seven problems were SWIFT II, FEMWATER, UNSAT2, USGS-30, and TOUGH. In general, linear problems involving scalars such as hydraulic head were accurately simulated oy both finite-difference and finite-element solution algorithms. Both types of codes produced accurate results even for complex geometries such as intersecting fractures. Difficulties were encountered in solving problems that involved nonlinear effects such as density-driven flow and unsaturated flow. In order to fully evaluate the accuracy of these codes, post processing of results using particle tracking algorithms and calculating fluxes were examined. This proved very valuable by uncovering disagreements among code results even though the hydraulic-head solutions had been in agreement. ! I I i f I
)
! iii l l
i t TABLE OF CONTENTS ! Pa21 I A B S T RA C T . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii L ,i FOREWORD ............................................................... xi ! r
- 1. INTRODUCTION ...................................................... i
+
1-1 ! s 4
- 2. DESCRIPTION OF COMPUTER CODES ..................................... 2-1 1
2.1 USGS-3D (Modular Three-Dimensional Finite-Difference I 1 Ground-Water Flow Model) ..................................... 2-1 : 2.2 TOUGH ...................................................... 2-2 l 2.3 UNSAT2 ....................................................... .. 2-3 2.4 SWIFT II ................................................. 2-4 ! 2.5 FEMWATER ..................................................... . . .. 2-9 4
- 3. HYDROCOIN LEVEL 1 PROBLEMS ........................................ 3-1 l 3.1 Case 1 ........................................................ 3-1 3.2 Case 2 ..... ................................................. 3-11 3.3 Case 3 ............................. .....................
i 3-45 3.4 Case 4 ....................................................... . 3-59 3.5 Case 5 .................................................... .. i 3-68 r 3.6 Case 6 ..................................................... . 3-78 3.7 Case 7 ...................................................... 3-114 j
- 4.
SUMMARY
AND CONCLUSIONS .........................................., 4-1 ! 4.1 Research Findings ............................................ 4-1 4.2 Regul a to ry Impl i cati on s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ( , 4-2
- I j
b REFERENCES ............................................................. R-1 ; 4 Appendix A, Corrections to Storativity, Leakance, and Transmissivity in Order to Simulate Radial Flow Conditions ................ A-1 Appendix B, Fracture Correction ........................................ r B-1 l Appendix C, Particle Tracking Algorithm ................................ C-1 f i f i t ) l $ i ; l ! i i m a ! i
- v !
i 1 [
. - . . - . - - - - . - - - . - - - m, ._ _ ,, , . _- --_.--__ , , , . , n . _ m _ , , , . - , - . - - . .c
LIST OF FIGURES r Figure Page ( 3.1.1 Conceptual model for Case 1 ..................................... 3-2 3.1.2 Assumptions for Case 1 .......................................... 3-3 3.1.3 Specifications for model output for Case 1 ...................... 3-5 3.1.4 Grid used to simulate Case 1 .................................... 3-7 3.1.5 Hydraulic-head distribution at r = Sm and Z = 2.5m for Case 1 ..................................................... 3-8 , 3.1.6 Hydraulic-head distribution at r = Sm and Z = Sm ! for Case 1 ...................................................... 3-9 3.1.7 Hydraulic-head distribution at 2 = 4m and t = 100s for Case 1 ...................................................... 3-10 3.2.1 Conceptual model and assumptions for Case 2 ..................... 3-12 3.2.2 Geometry of modeled :iomain for Case 2 ........................... 3-13 3.2.3 Specifications for model output for Case 2 ...................... 3-14 3.2.4 Fine grid used for SWIFT II simulation of Case 2 ................ 3-15 3.2.5 Medium grid useo for SWIFT II simulation of Case 2a . . . . . . . . . . . . . 3-16 3.2.6 Coarse grid used for SWIFT II simulation of Case 2 . . . . . . . . . . . . . . 3-17 3.2.7 Medium grid used for SWIFT II simulation of Case 2 .............. 3-18 3.2.8 Hydraulic-head profile for Z=0, Case 2 - SWIFT II ............... 3-19 . 3.2.9 Hydraulic-head profile for Z=-200, Case 2 - SWIFT II ............ 3-20 3.2.10 Hydraulic-head profile for Z=-400, Case 2 - SWIF1 II ............ 3-21 3.2.11 Hydraulic-head profile for Z=-600, Case 2 - SWIFT II ............ 3-22 3.2.12 Hydraulic-head profile for Z=-800, Case 2 - SWIFT II ............ 3-23 3.2.13 Hydraulic-head profile for Z=0, Case 2a - SWIFT II . . . . . . . . . . . . . . 3-24 : 3.2.14 Hydraulic-head profile for Z=-200, Case 2a - SWIFT II . . . . . . . . . . . 3-25 3.2.15 Hydraulic-head profile for Z=-400, Case 2a - SWIFT II ........... 3-26 3.2.16 Hydraulic-head profile for Z=-600, Case 2a - SWIFT II ........... 3-27 3,2.17 Hydraulic-head profile for Z=-800, Case 2a - SWIFT II . . . . . . . . . . . 3-28 3.2.18 Pathlines for fine grid, Case 2 - SWIFT II ...................... 3-29 3.2.19 Pathlines for medium grid, Case 2 - SWIFT II .................... 3-30
- 3. 2. 20 Pathlines for coarse grid, Case 2 - SWIFT II . . . . . . . . . . . . . . . . . . . . 3-31 '
3.2.21 Pathlines for Case 2 from GWHRT model ........................... 3-33 3.2.22 Case 2 - USGS grid ...................................... ....... 3-36 3.2.23 Hydraulic-head distribution at Z = Om for Cast 2 - USGS ......... 3-37 3.2.24 Hydraulic-head distribution at Z = -200m for Case 2 - USGS ............................................................ 3-38 3.2.25 Hydraulic-head distribution at Z = -400m for Case 2 - t USGS ............................................................ 3-39 3.2.26 Hydraulic-head distribution at Z = -600m for Case 2 - USGS ...................... ..................................... 3-40 3.2.27 Hyoraulic-head distribution at Z = -800m for Case 2 - , USGS ........ ................................................... 3-41 1 3.2.28 Hydraulic-head contours for Case 2 - USGS ....................... 3-42 ! 3.2.29 Model predicted pathlines for Case 2 - USGS ..................... 3-43 l 3.2.30 Coarse grid used for FEMWATER simulation of Case 2 .............. 3-46 3.2.31 Fine grid used for FEMWATER simulation of Case 2 ................ 3-47 3.2.32 Hydraulic-head distribution at Z=0m for Case 2-FEMWATER ......... 3-48 3.2.33 Hydraulic-head distribution at Z=-200m for Case 2-FEMWATER ...... 3-49 l 3.2.34 Hydraulic-head distribution at Z=-400m for Case 2-FEMWATER ...... 3-50 i vi I _ .._ _..__ _s_~ __
LIST OF FIGURES (Continued) Figure Page 3.2.35 Hydraulic-head distribution at Z=-600m for Case 2-FEMWATER ...... 3-51 3.2.36 Hydraulic-head distribution at Z=-800m for Case 2-FEMWATER ...... 3-52
- 3. 2. 3 7 Pathlines for coarse grid of Case 2 - FEMWATER . . . . . . . . . . . . . . . . . . 3-53 3.2.3t1 Pathlines for fine grid of Case 2 - FEMWATER .................... 3-54 3.3.1 Conceptual model, assumptions, and specification of output for Case 3 ........................................ ............. 3-56 3.3.2 Grid for Case 3 - UNSAT2 ........................................ 3-58 3.3.3 Model predicted water-table elevation for Case 3 . . . . . . . . . . . . . . . . 3-60 3.4.1 Conceptual model and assumptions for Case 4 ..................... 3-61 3.4.2 Grid used to simulate Case 4 - SWIFT II ......................... 3-63 3.4.3 Temperature rise vs. time, r = 0 meters, for Case 4 - SWIFT II .. 3-64 >
3.4.4 Temperature rise vs. distance, r = 0 meters, for Case 4 - SWIFT II ............................................... 3-65 3.4.5 Pressure rise vs. time, r = 0 meters, for Case 4 - SWIFT II ..... 3-66 3.4.6 Pressure rise vs. distance, r = 0 meters, for Case 4 - SWIFT II . 3-67 3.4.7 Pathlines for 100 year start time of Case 4 - SWIFT II .......... 3-69 3.4.8 Pathlines for 1,000-year start time of Case 4 - SWIFT II ........ 3-70 3.4.9 Pathlines for 10,000 year start time of Case 4 - SWIFT II ....... 3-71 3.5.1 Conceptual model and assumptions for Case 5 ..................... 3-73 3.5.2 Specifications for mouc1 output for Case 5 ...................... 3-74 3.5.3 Grid used to simulate Case 5 .................................... 3-76 3.5.4 Grid used to simulate Case 5a ................................... 3-77 3.5.5 Pressure profiles along Z = -50m and -100m for Case 5 ........... 3-79 3.5.6 Pressure profiles along Z = -150m and -200m for Case 6 .......... 3-80 3.5.7 Pressure profiles along Z = -250m and -300m for Case 5 .......... 3-81 3.5.8 Salt concentration profiles for Case 5 .......................... 3-82 3.5.9 Model predicted pathlines for Case 5 ............................ 3-83 3.5.10 Pressure profiles for Z = -50m and -100m for Case 5a ............ 3-84 3.5.11 Pressure profiles for Z = -150m and -200m for Case 5a ........... 3-85 3.5.12 Pressure profiles along Z = -250m and -300m for Case !.a ......... 3-86 3.5.13 Salt concentration profiles for Case 5a ......................... 3-87 3.5.14 Model predicted pathlines for Case 5a ........................... 3-88 3.5.15 Salt concentration contours for Case 5 .......................... 3-89 3.5.16 Salt concentration contours for Case 5a ....... ................. 3-90 3.6.1 Diagram of Case 6 including problem as.umptions ................. 3-91 3.6.2 Plan viewr, of Case 6 including hydraulic parameters ............. 3-92 3.6.3 Specifications for model output for Case 6 ...................... 3-93 3.6.4 Grid used for Case 6 - SWIFT II ................................. 3-95 3.6.5 Hydraulic-head contours for layer 1 of Case 6 - SWIFT II . . . . . . . . 3-96 3.6.6 Hydraulic-head contours for layer 2 of Case 6 - SWIFT II ........ 3-97 3.6.7 Hydraulic-head contours for layer 3 of Case 6 - SWIFT II ........ 3-98 3.6.8 Hydraulic-head contours for layer 4 of Case 6 - SWIFT II ........ 3-99 3.6.9 Hydraulic-head contours for layer 5 of Case 6 - SWIFT II ........ 3-100 3.6.10 Hydraulic-head profiles 1-4 for Case 6 - SWIFT II ............... 3-101 3.6.11 Hydraulic-head profiles 5-8 for Case 6 - SWIFT II ............... 3-102 vii
i i LIST OF FIGURES (Continued)
.F_i.gu re Page 3.6.12 Pathlines 1-4 for Case 6 - SWIFT II ............................. 3-103 3.6.13 Pathlines 5-8 for Case 6 - SWIFT II ............................. 3-104 3.6.14 Hydraulic-head contours, layer 1 of Case 6 - USGS ............... 3-107 3.6.15 Hydraulic-head contours, layer 2 of Case 6 - USGS ............... 3-108 3.6.16 Hydraulic-head contours, layer 3 of Case 6 - USGS ............... 3-109 3.6.17 Hydraulic-head contours, layer 4 of Case 6 - USGS ............... 3-110 3.6.18 Hydraulic-head contours, layer 5 of Case 6 - USGS ............... 3-111 3.6.19 Hydraulic-head profiles 1-4 for Case 6 - USGS ................... 3-112 3.6.20 Hydraulic-head profiles 5-8 for Case 6 - USGS ................... 3-113 3.7.1 Conceptual model for Case 7 ..................................... 3-115 3.7.2 Model grid for Case 7 - UNSAT2 .................................. 3-117 3.7.3 Hydraulic-head contours for variant 1 of Case 7 - UNSAT2 ........ 3-119 3.7.4 Streamlines for variant 1 of Case 7 - UNSAT2 .................... 3-120 3.7.5 Vertical velocities through center of blocks (Z = 12.5m) of Case 7 - UNSAT2 ................................................. 3-121 3.7.6 Vertica; selocities through top of blocks (Z = 15.0m) of Case 7 - UNSAT2 ................................................. 3-122 3.7.7 Vertical velocities through bottom of blocks (Z = 10.0m) of Case 7 - UNSAT2 ,................................................ 3-123 3.7.8 Hydraulic-head contours for variant 2 of Case 7 - UNSAT2 ........ 3-125 3.7.9 Streamlines for variant 2 of Case 7 - UNSAT2 .................... 3-126 3.7.10 Initial finite-element discretization for Case 7 -
FEMWATER ......................................................., 3-127 3.7.11 Final finite-element discretization for Case 7 - F E MW AT E R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-128 3.7.12 Hydraulic-head contours for variant 1 of Lase 7 - FEMWATER ...... 3-129 3.7.13 Hydraulic-head contours for variant 2 of Case 7 - FEMWATER . . . . . . 3-130 3.7.14 Vertical velocities along a horizontal profile located at bottom of trenches (Z=10m) of Case 7 - FEMWATER ................. 3-131 3.7.15 Vertical velocities along i norizontal profile located at midpoint of t renches (Z=12. 5m) of Case 7 - FEMWATER . . . . . . . . . . . . . 3-132 3.7.16 Vertical velocities along a horizontal profile located at top of trenches (Z=15m) of Case 7 - FEMWATER .................... 3-133 3.7.17 Particle trajectories for variant 1 of Case 7 - FEMWATER ........ 3-134 3.7.18 Particle trajectories for variant 2 of Case 7 - FEMWATER ........ 3-135 viii l
I l LIST OF TABLES ' Tabl e. Pg ; t 3.1.1 Input parameters for Case 1 ..................................... 3-4 i i , 3.2.1 Travel times and distances for pathlines predicted ' for Case 2 - SWIFT II ........................................... 3-34 3.2.2 Travel times and distances for pathlines predicted for Case 2 - USGS ............................................... 3-44 3.2.3 Travel times and distances for pathlines predicted for Case 2 - FEMWATER ........................................... 3-55 3.4.1 SWIFT II and HYOROCOIN pathline comparisons ..................... 3-72 3.6.1 Accumulated travel times and distances for individual paths predicted for Case 6 - SWIFT II ................................. 3-105 t i t 1 l l r 4 1 i i i ix . j
- i l
LIST OF APPENDIX A FIGURES Figure Pa21 t
- A.1 Block dimensions in radial and Cartesian coordinates used
. in adjusting storativity ........................................ A-5 A.2 Block dimensions in radial and Cartesian coordinates used in adjusting leakance .............,............................. A-6 A.3 Block dimensions in radial and Cartesian courainates used in adjusting transmissivity ..................................... A-7 LIST OF APPENDIX B FIGURES 8.1 Terms used in adjusting model parameters for fracture simulation ...................................................... B-3 t LIST OF APPENDIX C FIGURES : C.1 Path of integration of Equation 3 through a discontinuous velocity field .................................................. C-7 C.2 Domains of integrations of Equation 4 ........................... C-8 C.3 Hypgtheticc1 pathline calculated by periodic evaluation j of V ............................................................ C-9 , C4 Linear velocity interpolation ................................... C-10 C5 Finite-difference representation of diagonal flow ............... C-11 i i 4 i t t J l J
l FOREWORD t Tnis NUREG report is the first volume of a three-volume report documenting the '
! work of U.S. Nuclear Regulatory Commission staff and contractors as part of an -
1 international ground-water modeling project entitled "Hydrologic Code Intercom-parison Study" (HYDROCOIN). This first volume presents the results of the model verification phase (Level 1) of HYDROCOIN. It will be followed by Volume 2 dealing with model validation (Level 2) and Volume 3 dealing with uncertainty and sensitivity analysis (Level 3). 1 The computer codes used to date are SWIFT II, UNSAT2, TOUGH, FEMWATER, and
; USGS-30. NUREG/CR reports written as user's manuals for these codes are avail-able for purchase from the NRC and are listed in the references.
The work described in this NUREG report documents the combined efforts of the NRC project team led by Timothy J. McCartin, RES, and the Sandia National Laboratories project team composed of Dave Updegraff, Walt Beyeler, Rodger Miller, and Ken Brinster and led by Paul Davis. Ms. Louise Gallagher, NRC, served as the technical editor. Questions and comments on this report or the NRC activities in the HYDROCOIN effort should be directed to: Thomas J. Nicholson Waste Management Branch (Mail Stop NL-007) , Office of Nuclear Regulatory Research U.S. Nuclear Regulatory Commission j Washington, DC 20555 4 1 1 i 1 l Xi L
I l < l 1. INTRODUCTION , This report describes the computer simulations and analyses performed by the i HYOROCOIN project teams at the Sandia National Laboratories and at the NRC on the Level 1 HYOROCOIN test cases. HYOR0 COIN (Hydrologic Code Intercomparison Study) is an international project for studying ground-water modeling i strategies. The HYDROCOIN study is divided into three levels: Level 1--code verification (i.e., the subject of this report), Level 2--model validation, and 4 Level 3--uncertainty and sensitivity analyses. The Level 1 problems, in parti- , cular, were designed to test the numerical accuracy of the codes by using well-defined problems for which there are analytical solutions and/or facility for intercomparison of model results. These Level 1 test c;ses cover a wide range of hydrogeologic situations, including saturated and unsaturated flow, flow through media of highly contrasting permeabilities (e.g., fracture zones and ! argillaceous strata at low-level waste sites), and flow systems affected by
- thermal sources and the presence of variable ground-water densities. Model simulations reported in this document were previously transmitted to the ,
, HYOROCOIN Project Secretariat via magnetic tape for their use in preparing
, a final report of all 22 project teams' modeling results and their
- intercomparisons.
1 i The following sections provide a description of the computer codes used, a detailed discussion of each test case simulated, explanations of how these codes were used to model these test cases, and an analysis of the modeling results. In addition, appendices are provided to clarify the processing of , model inputs for Test Cases 1 and 2, and to describe the particle tracking algorithm used to produce the HYDROCOIN-required model outputs for many of the l cases. I , i i 4 t l 1-1 l ( i
- 2. DESCRIPTION OF COMPUTER CODES The NRC staff expects that a variety of ground-water models may be used in high-level and low-level waste disposal performance assessments. For example, modeling capabilities in the areas of saturated and unsaturated flow, fracture flow, and non-isothermal flow will be important. Therefore, a number of codes that are currently in use or expected to be used for performance assessments in these subject areas were selected for use in the Level 1 simulations. A discussion of each of the selected codes follows.
2.1 USGS-30 (Modular Three-Dimensional Finite-Difference Ground-Water Flow Model) HYDROCOIN Cases 1, 2, and 6 were simulated with the U.S. Geological Survey (USGS) Modular Three-Dimensional Finite-Difference Ground-Water Flow Model. (Ref. 1). The code solves the governing equation: (a/8x)(Kxxah/ax) + (a/8y)(Kyyah/ay) + (0/az)(K zz ah/0z) - W = S 38h/at where b = potentiometric head (L), K yx
= hydraulic conductivity 'L/T) in the x direction, K
yy
= hydraulic conductivity (L/T) in the y direction, K,7 = hydraulic conductivity (L/T) in the z direction, S, = specific storage (1/L),
t = time, W = source-sink terms (1/T), and x,y,z = Cartesian coordinates, aligned with the principal axes of the hydraulic conductivity tensor. The governing equation is discretized on an orthogonal, variably spaced, block-centered grid. The code provides no-flow and constant potential boundary conditions and six different source-sink terms, including specified and potential dependent condi-tions. The number of source-sink terms applied in one block is not limited by the code. Additionally, the discretized equation may be used in a linear form (for simulation of confined ground-water conditions) or in a nonlinear form (for simulation of phreatic ground-water conditions). Different forms may be combined within a single simulation, and the code may be allowed to convert from confined to phreatic conditions, or from phreatic to confined conditiores at any point or any time as conditions change. Matrix solution may be accomplished with either slice-successive overrelaxation (SSOR) on a strongly implicit procedure (SIP). 2-1
The USGS code is intended primarily as a water-resource model. It is limited to conditions of single phase, constant.-density flow and contains input / output structures designed specifically for representation of layered aquifer systems. The x Cartesian coordinate is, by convention, parallel to grid rows, the y coordinate parallel to grid columns, and the z coordinate normal to grid layers. The z coordinate is not necessarily vertical but the x, y, and z axes most correspond to the principal axes of the hydraulic conductivity tensor, which may take any orientation. The code is written with a modular structure and is available in Fortran IV and Fortran 77 versions. The modular structure of the code encourages modification of the existing code and addition of new code modules. Many new modules pro- - viding different source-sink functions hhve been proposed, and several (not employed in these applicauons) are currently available. The work completed with the USGS code was performed on an IBM Personal Computer AT using Microsoft Fortran, version 3.2. The code was obtained from the Holcomb Research Institute end differed slightly from the originally published version. Microsoft Fortran, version 3.2, is a subset of Fortran 77. Implementation of the code in the subset language required modification (completed by personnel of the Holcomb Research Institute) of eight format statements. 2.2 TOUGH The TOUGH code (Ref. 2) solves the two phase flow of air and water in vapor and liquid phaqes and heat flow in a fully coupled way. The formulation used in TOUGH is analogous to that used in multiphase, multicomponent geothermal or steam-flooded hydror:arbon reservoirs. The governing equations account for gaseous diffusion, Darcy flow, and capillary pressure effects. Vaporization and condensation with latent heat effects and conduction and convection of heat flow are included in the energy balance. Water, air, and rock are assumed to be in thermodynamic equilibrium at all times. The flow domain can include liquid, vapor, and two phase regions, indicating that the code can handle both saturated or unsaturated flow problems either individually or simultaneously. The thermophysical properties of water are represented by the International Formulation Committee's steam tables. Air is approximated as an ideal gas, and additivity of partial pressures is assu' _d for air-vapor mixtures. TOUGH solves three nonlinear partial differential equations simultaneously. These consist of the conservation equatioris for air, water, and heat. Air and water can be transported in both the liquid anJ vapor phases. Flow is governed by Darcy's law. ihe co:1e can simulate flow in one, two, or three dimensions because the method of solution is based on a very general integrated finite-difference method. Time stepping is accomplished by either the Crank-Nicolson or fully implicit procedure. The resulting nonlinear difference equations are linearized by the Newton Raphson technique. The resulting equations are solved by the Harwell matrix solver, which stores only the nonzero elements of a matrix. This reduces core storage. Attempts at using TOUGH for the Level 1 problems were made on the Sandia CYBER 70 Model 76 and CRAY-1 computers. 2-2
r 2.3 UNSAT2 The UNSAT2 code (Refs. 3 and 4) was used to solve HYDROCOIN Level 1 Cases 3 and 7. UNSAT2 is formulated to solve two-dimensional transient, saturated-unsaturated flow problems in porous media. The governing partial differential f equation used by UNSAT2 is: C($) + pS, a +S=0 5x 4 r(*)Ks (i J) + Kr (*)Es (i J) where C = specific water capacity, K = relative hydraulic conductivity, 7 K,(4 ) = saturated hydraulic conductivity, S = sink terms, S, = specific storage, t = time, xj = spatial coordinates, h = 1 for saturated conditions and 0 for unsaturated conditions, and
$ = pressure head (tension).
The fact that S is zero when the porous medium is unsaturated assumes that the soil is incompressible. - UNSAT2 can handle various types of boundary conditions. These include constant pressuni head, prescribed flux, impermeable or zero flur., seapage faces (either set at the initiation of the problem or calculated at '.ater time steps), and fluxes determined by atmospheric conditions such as evaporation or infiltration. Several internal sources or sinks can ta accommodated by UNSAT2. These include discharging or recharging wells and water uptake by the roots of plants. UNSAT2 is probably unique with respect to water uptake by plants. However, none of the features mentioned above was used for the HYDROCOIN modeling. The equations describing unsaturated-saturated flow are discretized by the Galerkin finite-element method. Elements can consist of either quadrilaterals or triangles. Triangles are very useful for modeling irregularly shaped regions. However, before performing integrations of the shape functions ever the elements, UNSAT2 converts all quadrilateral elements to triangles. The time derivative is discretized by the finite-difference method. However, UNSAT2 does not use the finite difference of pressure head at each node of the grid but uses weighted averages of the time derivative over the entire flow region. The purpose for this is to improve the convergence properties of numerical solution. The equa-tions are integrated over t *me in one of either two methods. A time-centered 2-3
scheme is recommended for most problems. However, when S, equals zero and when nodes change from unsaturated to saturated conditions, a backward-in-time discretization is used. The code, as supplied to Sandia, was set up to rv'. on IBM computers. With .ne removal of REAL*8 cards from the Fortran code, the code ran on Sandia's CDC CYBER 155 Model 855 and CRAY-1 computers. 2.4 SWIFT II Cases 2, 4, 5, and 7 were simulated with the SWIFT 11 computer code (Refs. 5 and 6). SWIFT II (Sandia Waste Isolation, Flow and Transport) is a fully tran-sient, three-dimensional finite-difference model, which simultaneously solves coupled equations for transport in geologic media. It is designed to simulate flow and transport processes in both single and double porous media. The pro-cesses considered are: Fluid flow Heat transport Brine migration Radionuclide-chain transport The first three processes are coupled via fluid density, viscosity, and porosity. Together they provide the velocity field on which the radionuclide transport depends. The transport equations are obtained by combining the appropriate continuity and constitutive relations. Sink terms (F) are included for fracture zones in which losses to the rock matrix are significant. The resulting relations may be stated as follows:1 Fluid:
-V-(pu) -
q qq conduction production sink / source 2
+
Rg ry = h(4p) (1) salt loss to dissolution matrix 1All terms are defined after Equation 16. 2This term refers to a sink / source other than a well. A positive sign Jenotes a sink, and a negative sign denotes a source. 2-4
Heat: ;
- V-(pHu) + ) ~
NI 9 ~ ~ t' ([H N9 4H ' convection conduction / injected produced sink / dispersion enthalpy 5 enthalpy source r
- = O (Hry + r g) g [$pU + (1 - 9)pR RU3 loss to matrix accumulation in fluid and rock (2)
Brine: I A A A A
-V-(pCu) +
v.(pgCVC) - Cgy Cq t convection dispersion / injected produced 4 diffusion brine brine .' ^ a
+ R C
(Cry + FC)
- E I3) salt loss to matrix accumulation 4 dissolution Radionuclide r:
j - U (pCpu) + - r Cr4 ~ 9r
- Qwr convection v.(pgC[VC) disper ion / produced sink / wasa diffusion component source la .ch '
r N :
+ 1 k rsAs[$pC, + (1 - $)p gW,] i
- (C IrW+P) r
; loss to 8*henerationofcomponent matrix r by decay of s l ,.
l Ap [$pCp + (1 - $)pRW r] = g [9pCr * (1 ~ *)P W3 Rr decay of component r accumulation (4) ,. j Several quantities in Equations 1 through 4 require further definition in terms of the basic parameters. The tensors in Equations 2, 3, and 4 are defined as
. F sums of dispersion and molecular terms:
r 4 E =D+Dy (5) and Eg =Spcp +s m (6) 3This is a source term since, by the adopted sign convention, the rate of fluid injection is inherently negative. 2-5 i i
- - _ . _ .._,.m.
where D jj=aVT ij 64 ) + (at - a T)Vi Vj /IVI (7) in a Cartesian system. Also, sorption of radionuclides is it.cluded via an assumption of a nonlinear Freundlich equilibrium isotherm: W r * *r(pCp ) r (8) Equations 1 through 4 are coupled by auxiliary relations for: Darcy Flux: u=-(k/p)-(vp-fvz) c (9) Porosity: 0 = og[1 + c p(p p )] g (10) Fluid Density: p=pg [1 + C,(p - og ) - cT (T - Tg ) + cC 3 (11) Fluid Viscosity: p = pR( ) exp (B(C)(T 1 - T R 1)] (12) Fluid Enthalpy: H = U, + U + p/p (13) Fluid Internal Energy: U = cp (T - Tg ) (14) (This equation assumespc is constant with regard to temperature.) Rock Internal Energy: U=cpR(T R - T,) (15) 2-6
where parametercc in Equation 11 is defined in terms of an input density range (p y PN )!Po* c c * (PI ~P)PN o (10) Symbols B Viscosity parameter C Concentration of inert contaminant c Specific heat of the fluid p c pR Specific heat of the rock (single porosity) or of the fracture-fill material (dual porosity) Cp Concentration of radioactive (trace) components C R Compressibility of the pores (single porosity) or of the fractures (dual porosity) c T Coefficient of thermal expansion C, Compressibility of the fluid Dj Dispersion / diffusion tensor D, Molecular diffusion E X Dispersion or conduction / dispersion tensor for heat (X = H), brine
^
(X = C), or radionuclide (X = Cr ) within the global system g Acceleration of gravity g Units conversion factor equal to g for the English system and equal c to unity for the SI system H Fluid enthalpy Hg Enthalpy of injected fluid I Unit tensor k Permeability tensor for the global system k rs Product of branching radionuclide and daughter parent mass fraction K, Heat conductivity tensor for fluid and rock (single porosity) or fluid and fracture-fill material (dual porosity) p Pressure 4 P, Reference pressure for system, initial pressure for the unsteady-state aquifer model and aquifer boundary condition for the steady-state aquifer model q Rate of fluid withdrawal 2-7 ;
7 q, Radionuclide source due to waste leaching qX Sink / source other than a well for fluid (X = W), heat (X = H), and radionuclide (X = r) R C Brine source rate due to salt dissolution Rf Fluid source rate due to salt dissolution t Time T Temperature T R Ambient temperature of rock surrounding well bore T, Reference temperature of system, interface temperature between system and over/underburden, and surface temperature for radiation model U Mass specific internal energy of the fluid U R Mar,s-specific internal energy of the rock (single porosity) or of the fracture-fill material (dual porosity) U n Me.ss-specific internal energy of the fluid at reference fluid conditions lvl Magnitude of v4 ) Darcy velocity v$) W r Solid phase concentration of radionuclide r W s Solid phase concentration of radionuclide s x,y,z Cartesian coordinates a L ngitudinal dispersivity l a i Transverse dispersivity F Total loss to the rock matrix for fluid (X = W), heat (X = H), brine X (X = C), or radionuclide (X = r) 6 4,) Dirac delta function r1 Freundlich isotherm parameter x Freundlich isotherm parameter K Modified Freundlich isotherm parameter r x R Thermal diffusivity of rock surrounding the well bore A Decay constant p Viscosity pR Viscosity parameter p Fluid density pg Fluid density at reference temperature and pressure and unit brine concentration 2-8
pH Fluid density at reference temperature and pressure and zero brine concentration p, Fluid density for the initial conditions pg Formation density
$ Porosity t f
(g Porosity at the reference pressure Furthermore, an internal energy U n is included in Equation 13 to account for the difference in reference conditions as specified by the analyst and the reference conditions specified internally for the enthalpy. SWIFT II allows the option of using either a direct solution technique or the iterative two-line successivo-overrelaxation technique. All the SWIFT II simulations were performed on the CRAY-1 computer at Sandia National Laboratories. 2.5 FEWATER Cases 2 and 7 were simulated with the FEWATER computer program (Ref. 7). FEWATER solves two-dimensional transient saturated-unsaturated flow problems under isothermal conditions in porous media. The governing partial differential equation solved for in the FEWATER program is:
- 9. ($ . (vh + VZ)]= a+esfh where 4 h = pressure head, K = hydraulic conductivity tensor,
] n, = effective porosity, c t = time, ' x,z = spatial coordinates, Z = elevation above datum, o = modified compressibility of the medium (compressibility of consoli-dation times the density of water times gravity), ' 4 S = modified compressibility of water (cort.pressibility of water times the density of water times gravity), and 0 = moisture content. The FEWATER program can apply to one or more of the following types of boundary conditions: constant or time-dependent pressure head, prescribed flux, and a i [ l 2-9 1 i
~ - - -
rainfall or seepage flow. Any portion of the boundary not assigned a boundary condition is automatically assigned a zero flux condition by FEMWATER. , The governing partial differential equation is solved using the Galerkin finite-element method. A quadrilateral bilinear finite-element is used. Additionally, the Darcy velocities are solved for by applying the Galerkin finite-element method to Darcy's law. This technique is used to avoid discontinuities of the flow field at the nodes. The time derivative is approximated using a finite-cifference method. There are three possible choices within the FEMWATER program for this approximation. These are Crank-Nicolson formulation (central differ-ence), backward difference formulation, and mid-difference formulation. In addition, FEMWATER assumes that the off-diagonal elements of the hydraulic conductivity tensor are zero, implying that the grid should be oriented along the direction of the principal axes of the hydraulic conductivity tensor. The computer program supplied to the NRC was set up to run on a CDC computer. How-ever, with minor modifications (removal of LEVEL and PROGRAM cards), the code wasconvertedtoanIB$$p70machineattheNationalInstitutesofHealth(the fact that the program was originally developed on an IBM computer obviously made this transition easier). I i I l t i ! I l t s 2-10 ! t i !
1
- 3. HYOR0 COIN LEVEL 1 PROBLEMS The following section contains a brief description of each of the problems along with a discussion of which codes were applied to the problem and how well the results matched analytical solutions and results from other codes.
3.1 Case 1 Problem Oescription Case 1 involves transient radial flow from a well in a saturated porous medium underlain by a fracture (see Figs. 3.1.1 and 3.1.2 and Table 3.1.1). The pur-pose of the problem is to test a code's ability to simulate transient flow from a well having a finite radius into a region with a finite boundary. Code results are to be compared with an analytical solution for the times and locations given in Figure 3.1.3. Simulations Performed This problem was simulated with the USGS code. Completion of Case 1 with this code initially presented two problems: i
- 1. The code employs a rectangular finite-difference grid and contains no options for representation of radial flow.
- 2. The code, in order to accommodate heterogeneous conditions, employs the harmonic mean of transmissivities specified in adjacent blocks (rather j than the actual value specified) in solving the governing equation.
The model may be corrected to represent the conditions of radial flow by spe-cifying values of transmissivity, storativity, and leakance that vary as func-tions of the radius. The corrections are derived in Appendix A; only the results are listed here. S'g,) = S3rg0z) L'g,) = 2r gK'/(Oz) + Oz) j) f' j ,j = T)(r y - r )/1n(r 3/r ) g g g g i 3-1
i l I I ha(t)
" OREHOLE h=0 ,
/N////////////// ' HHHHHHHHHHHHHHH/HHH//HH/// A r V
.Z ROCK MATRIX a
d
= b =
4 HH/H/HHHHHHHHHHHHH/HHHH/HH/HHHHf d//HH/HHHH//// FRACTURE a RADIUS OF BOREHOLE (m) b RADIAL DISTANCE TO CONSTANT HEAD BOUNDARY (m) d THICKNESS OF ROCK MATRIX (m) h ASYMPTOTIC VALUE OF ha(t) (m) ha(t) IMPOSED HEAD IN BOREHOLE (m) r CYLINDRICAL RADIAL COORDINATE (m) z AXI AL COO.:31NATE (m) Figure 3.1.1 Conceptual model for Case 1. 3-2
ha(t) V h=0 l 3 v
//H/H/HHHHJ I /HH HHHHHHHHHHHHHHHHH, d -
U l* d 1 l
= b =
n
/H/HHHHHH/HHHHHHHH/HH/HHHHHHH/HH/H/H///////H/H, ASSUHPTIONS - The rock is saturated - Isotropic and homogeneous permeable medium - Flow can be described by Darcy's law - A prescribed time-dependent head is maintained in borehole relative to a fixed head at a radial distance b - Rock is confined between impermeable horizontal boundaries - Transmissivity of fracture = 10.s,2fs - Hydraulic conductivity of matrix = 10 8m/s Figure 3.1.2 Assumptions for Case 1, 3-3
n s Table 3.1.1 Input parameters for Case 1. Symbol Parameter Value T Transmissivity of fracture 10.s ,2 s 1 S Storage coefficient of fracture 10 10 K Hydraulic conductivity of matrix 10 'as 1 S, Specific storage of matrix 10 7m1 a Radius of borehole 0.1m b Radial distance to boundary 10m d Thickness of rock matrix Sm t Time constant for borehole head 0.ls where ,
= thickness of layer j, Dz)
Dz),y= thickness of layer j+1, K' = vertical hydraulic conductivity (constant), 1.'j,3 = corrected leakance for column i of layer j,
= radius to the block center, column i, rg rj,1 = radius to the block center, column i+1, S'g,) = corrected storativity for column i of layer j, S, = storativity of the layer (constant),
T = transmissivity of layer j (constant), and T f g,) = corrected transmissivity for flow from column i to column i+1 inlayerj. The harmonic mean transmissivities employed in the calculations, however, do l not correctly preserve the radial variation in transmissivity. The specified transmissivity must be calculated so that the harmonic mean of transmissiv-111es specified in adjacent nodes provide T', the correct radial variation. The correction, derived in Appendix A, leads to: Tg,y = Drg,yT g T/2Tg in(rg,y/r )g - OrgT 3-4
ha(t) v h=0 i '
//NNN/NN/M I
/NN, NNNNNNNNN/NN/N/N/NN/
d r Z d I
= b =
n
/HHHHHHHH/HHHHHHHHH/HHHH/H/HHHHH/HHHHHH/Hi REQUIRED OUTPUT - HYDRAULIC 4
HEAD IN THE MATRIX AS A FUNCTION OF TIME UP TO 10 s FOR r=5 m AND z =2.5 m
- HYDRAULIC HEAD IN THE MATRIX AS A FUNCTION OF TIME UP TO 10's FOR r=5 m AND z=5 m - HYDRAULIC HEAD IN THE MATRIX AS A FUNCTION OF DISTANCE BETWEEN r = a = 0.1 m AND r = b = 10m AT A TIME OF 100 s Figure 3.1.3 Specifications for model output for Case 1.
3-5
where Or g = block dimension of column i, measured along a row, Drj y = block dimension of column i+1, measured along a row, Tj = transmissivity to be specified in column i, and T$ y = transmissivity to be specified in column i+1. Other terms are defined above. Application of this equation requires the modeler to supply Tg. Different T, but all combine values of T$ lead to different values for later terms of to produce the same values of T'. The Case 1 problem is represented with a grid (shown in Fig. 3.1.4) consisting of 1 row, 50 layers, and 81 columns. The grid represents a 15-degree slice of the aquifer, with the well radius at the block center in column 1 and the held-potential boundary at the block center in column 81. Block spacing and time increments were selected to provide computed values at each point used in the problem description. The time-dependent hydraulic head at the well is represented with the General Head Module (a source-sink function) provided as an option in the code. This option calculates the volumetric flux to a block from a point outside the model grid, assuming that the flux is proportional to the difference between the head calculated in the node and a reference head existing outside the node. The reference head and the factor of proportionality must be specified for each node where the function is applied and may vary as a discrete function of time. The reference value is specified as the average value of the het.d in the well over each discretized time period. The factor of proportionality is set at an arbitrarily high, constant value in order to fix the head in column 1 at the reference value. The problem was solved to a convergence criterion of 1.0x10 8 The computer hydraulic head compared reasonably with the analytic solution (Figs. 3.1.5, 3.1.6, and 3.1.7), but problems arose in calculating the flux from the well (represented by flows to column 1 from the general hcad boundary). The flows from the well were calculated by the model (after completion of the solution in a double precision calculation) as a single precision variable. The difference between the reference heads and the model-calculated heads in column 1 was zero (single precision) at all points. The flows calculated by the model were also, therefore, zero. Modifying the code to complete the cal-culations at the double precision level alleviated the problem but led to discovery of a second problem. That is, the convergence criterion of 1.0x10 8 was too large. The heads calculated in column 1 were nearly equal to the reference heads, and variations of 1.0x10 8 led to errors in the calculated flow that not only varied either side of zero but were several orders of mag-nitude larger than the other mass-balance terms. The convergence criterion 3-6
5 o s 3, 5 2
, s s
4 0 1 l 0 8 I 3 5 l l 5 7 l I 7 0 7 l - - _ _ - _ _L I 0 7 e l 6 e l I 6 0 l I I 0 9 I 8 5 5 l l I 5 5 1 e 0 5 I
- _ _L I 0 5
C s a t e 5 I I I 5 a 4 I 4 l u
- - _ - _ _ R i m
s 0 0 4 I I l I 4 o - t _ d 5 l
- - _ 5 e
s 3 l l I 3 u d i r 0 - - _ G
- _ 0 l
3
- - _ _ _ L_ I 3
4 1 5 l 5 . 2 I I 2 3
- - _ _ e
- _ _ _ r u
0 g 2 l l l I 0 i 2 F - 5 1 l
- - _ - I
_ _ _ l _ I 5 1 oI - o t
- _ 1_
I
- - _ _ _ 1_ t 5 I l I l I 5 _
Z Ty
L O.18 . O.16 -
- ANALYTICAL SOLUTION O
0 USGS MODULAR CODE - k O.14 W O - I O.12 - , S 0.1 - 0 - w J O.08 - h O.06
>g- 0.04 l -
0.02 -
- ' I ' I ' ' '
O 10 30 100 300 1000 3000 10,000 TIME (seconds) Figure 3.1.5 flydraulic-head distribution at r = Sm and Z = 2.5m for Case 1. i I
0.18 - 0.1s - ANALYTICAL SOLUTION -
'" ~ -
USGS MODULAR CODE o O
< 0.12 - -
W
- I O o,1 - -
s w a E m 0.08 - o 1 0.06 - 0.04 - 0.02 - o .......! . . ......I . . . . ....I . . . . . . . . 10 30 100 300 1000 3000 10.000 TIME (seconds) Figure 3.1.6 Hydraulic-iiead distribution at r = Sm and Z = Sm for Case 1. r
^
1 5' I I I I I I I I I
.o. . .
ANALYTICAL SOLUTION
.06 -
O USGS MODULAR CODE
.or O
y .oe - y 9 _ w J .os - O 3 E - O .04 - I ~
.03 -
~
.02 ~
~
.01 -
0 0 1 2 3 4 5 6 7 8 9 to r (meters) Figure 3.1.7 Hydraulic-head distribution at z = 4m and t = 100s for Case 1.
was lowered to 1.0x10 8 and the model failed to converge below 6.7x10 8, where the errors were still large enough to dominate the mass balance. Therefore, this report does not contain flow r.alculations for Case 1. 3.2 Case 2 ' Problem Description A vertical plane cut by two intersecting fractures is to be simulated in this case (see Figs. 3.2.1 and 3.2.2). The problem is designed to represent an ' idealized fractured crystalline rock mass. This rock mass is assumed to be saturated, with the ground-water flow system being in steady-state conditionc. The purpose of this problem is to test the abilities of different codes to handle large permeability contrasts. In addition, the effect of discretiza-tion is to be investigated by requiring each participant to solve the problert with at least two grids. Listed in Figure 3.2.3 are the specifications for the output that are to be used to compare model results. ; Simulations Performed Three sets of simulations were performed using the SWIFT II code, the USGS code, and the FEMWATER code. The geometry of the intersecting fractures of r this case is more easily represented with an irregular grid of the finite- ' element method than with a regular grid of the finite-difference method. The finite-difference codes (SWIFT II and USGS) required more grid manipulation to l represent the fractures than the finite-element code (FEMWATER). Although the finite-difference grids were manipulated for both the SWIFT II and USGS codes, a much greater effort to accurately represent the fracturas was made in the ; simulation with the USGS code. Following is a description of the application ! of these codes to Case 2. i SWIFT II Three different grids were eLployed in solving this problem (see Figs. 3.2.4, 3.2.5, and 3.2.6). The number of nodes used was 4136, 2925, and 1148, respec-tively. Modifications were made to the elevation and thickness of the top row of blocks in an attempt to simulate the irregular top surface of the problem.
', Fractures were treated simply by assigning the hydraulic conductivity of the fracture to any block that was crossed by a fracture (Figs. 3.2.4, 3.2.5, and <
3.2.6). Initially, in the case of the medium grid, fracture G (Fig. 3.2.1) crossed only blocks lying on a diagonal (Fig. 3.2.5). Because of the finite-difference scheme, this resulted in fracture blocks not being directly coupled to each other. Therefore, additional blocks were designated as fractures in order to allow for a direct connection of blocks lying along the fracture (Fig. 3.2.7). Results of both representations are included herein. All the model results are shown in Figures 3.2.8-3.2.20. ! An analytical solution for this case does not exist. However, two bases for comparison exist. One is to compare the SWIFT II results for the different grids, and the other is to compare the SWIFT II results to those from other models. Shown in Figures 3.2.8-3.2.12 are the hydraulic-head profiles for the l 3-11 : l t
F A g D E ROCK MATRIX FRACTURE ZONES C B ASSUMPTIONS
- Both matrix and fracture zones are isotropic and homogeneous permeable media
- Flow can be describeo by Darcy's law
- Matrix and fracture zones are saturated
- Steaoy-state flow
- Hydraulic conductivity in matrix (10.s m/s) 100 times lower than in fracture zone (10 8m/s)
- No-flow boundaries along A-B, B-C, and C-D
- Ground-water level assumed to coincide with top boundary (A-G-F-E-D)
- Effective porosity for fracture, and matrix is .03 Figure 3.2.1 Conceptual model and assumptions for Case 2.
3-12
h 12-5 89 34 67 l 16 19 I 18 ' 1 15 14 13 121110 POINT- X Z I 1 0 150 !
; 2 10 ' 150 ;
3 395 100 ( 4 405 100 5 800 150 6 1192.5 100 7 1207.5 100 ) 8 1590 150 1 9 1600 150 ' 10 1600 -1000 11 1505 -1000 i 12 1495 -1000 13 1007.5 -1000 ' 14 992.5 -1000 i 15 0 -1000 l 16 1071.3462 -566.3459 i 17 1084.0385 -579.0383 ) 18 1082.5 -587.5 l 19 1069.8077 -574.8077 ! l l Figure 3.2.2 Geometry of modeled domain for Case 2. i 3-13 i
A 1 - A- -
/ 4- A, B -
- - B' C ~ ~
/ x C' '
3 D - - - D' E - -
/ E'
/
Required output (from at least two levels of discretization)
- Distribution of hydraulic head along five horizontal lines: A-A',
B-B', C-C', D-D', and E-E'
- Particle tracking with starting points 1-4 according to figure above Figure 3.2.3 Specifications for model output for Case 2.
l 3-14 l
~~~' ~~~~~ ~~
__m==*~"
~___ _ ____
a ri m Il iI II IIIE I I
~~
m 11 1111IllE I I J I II IIIIII ! I I m II iI II11II EaI IIIIIIEI I I I I m 11 I I I I I I EJ l I u II IIIIIIEI i 1 m 11 IliII uI I 1 m Ii iiiII II I I
. 11 IIIII II I I m Il ll ll l !I I I I m ll IIIIIEII I I s II IIII hMiI I I
. II IIII-.II I i:
m Il iIII EIII I II -~~
. II IllI EIII I II
. II II11EIII I II 11 IIII EI11 1 11
_ lI iIIIEIII I II 11 1II _IiI I II 11 IllMIIII I II 11 IIIEIIII I 1I 11 IIFEIIII I II II IIIEIIII I.II
--_, I E II Il IIII 11 11 Iin iIE_IIIII IIII 11III IIII l IIIII II11 l IIIII IIII w 11111 lill II I I I I I II 5 .. - ,_,_._ li 111111 llIIII IIII IIII II 1 I II III II II Il1 - 11III 11II
-_.__-._.--._ iIIN MIIII Iii1 1 I I R-_ _1 I I II II 11 i E_ _iI IIII Iii II II III
[_I liII Ill I l l i __ IIII II u 1 III i 1 II II I I IIIII il II I IIIIII I II I IIIIII 1 l'_I+I+1!!ll I I II 111111 IIIIII l'>f 1 I I II ~~~~~ I II I11111 I II l l 11 I IIIII I II
._.- 1 II I111II I 11 i IIiII I II
]II ..IIIIIIi
.._ _ _. I II
.. IIIIII1 I i1 a 1 E_____
FRACTURE BLocxs Figure 3.2.4 Fine grid used for SWIFT II simulation of Case 2.
9 i i l 1 l I E I I E - 1 E e
- -. g
.- - O l E O a
I G I E y 9 e N 1 m 0 l >= I E O $
- u
- - .m.--------- --------- . - --- --- ------ E 6
,_____________________ s,.
, C i
--+--- -------- -------------- O M ---------------
M
------------------ -, m g -
WM um.m.m.m men,name,mue,n m..mn.um.um m.mwwe.>_m=-->=.=-n=ms.- na -.mamus .umm - .am. e .u men,in __,_-._-_._ _ --_-._._-- ._ -_ . - . _._- .-_ e
----._------------------ ---------- N I
_--------- m -,.m.- 1 i
. i I i
. .i -
m I h 3 I I iE i 1 e i . L 6 O I i I l _
- i. i V
{ i 1 . g I I I E !I e 3 I i i i ! I l !
. .! .V.
; i
~
i @
]
l i .
! f?
3 1 i E~ i i. .
. 4, 1 1 i i i x
1 1 1 m 1 !! i i W l E ! . I i m I ! I i i 4 I - L I 3 l : . . . e I I + a= 1 I
" l 1
l I I i f I t I I I I I I I I 3-16
I I I I I I IIIIIII
~I I I I I l IIIIIII l I I I I I ___
l _ _ __ _ __
\
\,
_ _ __ _ __ am ox . so~ uo g __ _ _ _L .m l s
,-+
/j s
.b
=
___ _ _ g it e
- ___ _ _ a 7
7 __-___ ___ _ ._ e _____ ___ _ _ s W s o
! I I II I l iI h I II I l I f ; I I II I I I ~~
l I I II I I I l i II I I l . 1 I II I I l
; I II I 1 1 I II I I I ! I II I I I
4 I I II I I I I 3-17 .i
I I II iiIII III I _. ..._.LIJ l I I l_.III I I II IIIII III I ITIII II ~I I I I I II IIlII III I IIIII II Ii1 1 1 II IIIII II1 1 11III lI III I II IIIII III I I11II II II ~~- 1 1 11 I11II III I IlIll II ygg__Qy___ggg 7 p ____++p.H+H---+H-+----H-H+- .___- u 2: ___ _- -__ _ O . i ____ __ ___ _ ___g N
----- -- --- --_-W 0
____ g j I II I I III 1 II I ,_ _,, , _ , y ___ W . EE __ ___ i __ . . _ = e t g y ____ 7, E y e i
, i
___ EEI - i _ --__-_- y _ __- e __ _ - __._,_- e
~~ ~
l ___ i 11 11 1 y
=
i _____ _ _ ____ L__._ i i ! r _y __
~
i _,_.-m i - 7 z 7___7__ 7 i I m s
!II1l! I II Ii m L._ _ .-- _I _ - ._ . _ . _ .._.._,_
4 g fl I J lIi l I III I I II II L" II I I Il! II l!II l lI II II I I I1i I I -
! II i II II II I I III IJ ! III iI II fII I i1 1I i III II II I IIi i II II ~ lII II 1I I ! I i I I I I III II II I IIl I I II II1 1I II I i 1 1 I I I III II II II I I II III II II I I I I I I I I I II II I I I 1 II III II II ~~
I I I I II III Il I1~ I II I I II III II II I I I I I I I 'II II 1[ 3-18
HYDROCOIN LEVEL 1 - CASE 2 HYDRAULIC HEAD DISTRIBUTION AT Z = 0 METRES I I I I I I I 14o - _
~ ~
GWHRT MODEL
. . . . . . . sWirT II MODEL 81 m 56 _
~ -
-- - SWIFT 18 MODEL 65 m 45 g 130 g SWIFT II MODEL 41 x 28 /
y g -- o - \ /* _ g - N / - I -
\ _
w o 120 - \ - h h -
\ \ /
g \ - N \ -
\ g -
s N
~
(%/ \
* ~ -
I I I I I i 1 0 200 400 soo soo iooo 1200 1400 1600 DISTANCE (m) Figure 3.2.8 Hydraulic-head profile 'ior Z=0, Case 2-SWIFT Is.
HYDROCOIN LEVEL 1 - CASE 2 HYDRAULIC HEAD DISTRIBUTION AT Z = -200 METRES I I I I I I I - 140 - - GWHRT MODEL [ _ . . . . . . SWFT N MODEL 81 m 56 _
-- SWIFT 18 MODEL 65 x 45 y 130
-- SWIFT 88 MODEL 41 x 26 ]
! O - ~ til x - - T 120 -'N % B
~O
~ s ' , .~==
s s - 1 % % / r % / I tiO
- N [4 g
4 s g',, *'.s Ng,/
,f- --
\/ -
1 l 1 1 I I I 0 200 400 600 600 1000 1200 1400 1600 DISTANCE (m) Figure 3.2.9 Hydraulic-head profile for Z=-200, Case 2-SWIFT II.
HYDROCOIN LEVEL 1 - CASE 2 HYDRAULIC HEAD DISTRIBUTION AT Z=-400m i I I I I I i 140 -
~
GwHRT MODEL - _ ...... SwtFT H MODEL 81 m 56 -
~ ~
_ - == - Sw!FT 88 MODEL 65 a 45
~
w E 130 -
-- swift 88 MODEL 41 a 28 _
o - 4 us - m _ E
~
9s 120 - -
~
3 ;r % K
am. W-o "**** % % -,,,,,#
# ~
1
~ s *= _ _. _
~
1 o " " a et u n gg n % su ,N p , _ { -
~
'00 i 1 1 I I I 1 o 20c 400 soo 800 1000 1200 1400 1600 DISTANCE (m)
Figure 3.2.10 Hydraulic-head profile for Z=-400. Case 2-5WIFT II.
HYDROCOIN LEVEL 1 - CASE 2 HYDRAUL'. HEAD DISTRIBUTION AT Z=-600m l 1 1 I I I I _ iso -
~
~
c wtety m
****** SwwT H tdODEL 81 m 64 _
n
- ' N O WL 65 a 45 5 iso - -
NL 41 x 2s ~
~
o g - E y o
~
N - '
" 3 tro -
3 a- * % "~ m - O -
'- - ~ .
I ?.****r.n o " ' n r. *,,,%
** *v . r , , ~~_~
i8o - . ,
*n n n n r. . . . . . . . . . . . . . . . . . . . . .
~
I I 1 . I I I o roo 400 eoo soo sooo 1200 1400 isoo DISTANCE (m) Figure 3.2.11 Hydraulic-tead profile for Z=-600. Case 2-SWIFT II. l.
- - - - - _ - "--"w w --4 w- w , . - - , - - - . .w e
HYDROCOIN LEVEL 1 - CASE 2 HYDRAULIC HEAD DISTRit3UTION AT Z=-800m { l i i i I I I _ 140 - ~
~
OWHRT MODEL [ _ * * * * *
- SWIFT II MODEL 81 x 56 _
m SWIFT II MOLEL 65 x 45 -
*E iso
-- SWIFT 11 MODEL 41 x 28
- o - ~
l < - i m _ _ l I l o - - w ] 120 4 w o
-~
~
cr _
~ -
~ ____ _-
? n====~~- ____
l z ,,o _
~ _ _ _ _ _ ., ., r.v va
'O ~ -
l l I I I ! I o 200 400 coo 800 1000 1200 1400 1800 DISTANCE (m) Figure 3.2.12 Hydraulic-head profile for Z=-800, Case 2-SWIFT II.
HYDROCOIN LEVEL 1 - CASE 2a HYDRAULIC HEAD DISTRIBUTION AT Z=Om I I I l l l 1 _
~
140 -
~
GWHRT MODEL ,
= * * * *
- SWIFT 11 MODEL 81 x 56 _
) - m
- --- SWIFT II MODEL 65 x 45 (s s E 130 l
gN /,_ \
-- SWIFT 11 MODEL 41 x 28
[y..
< kN y
% / \ / -
/
T c
- 120
\g f
/ I j -
g o - N f f i
< \% f
./. / -
4gN' x / 7 g / : I 220 _ .N z? s / _
.N
, N .
.,..,= 9 _
~
' I I l i I I I 600 800 1000 146. 1400 1600 l 0 200 400 DISTANCE (m) i Figure 3.2.13 Hydraulic-head profile for Z=0, Case 2a-SWIFT II.
HYDROCOIN LEVEL 1 - CASE 2a HYDRAULIC HEAD DISTRIBUTION AT Z=-200m I I I I i l I . 140 - - GWHRT f e JDEL ,
. . . . . . SWIFT 11 MODEL 81 x 56 -
- - . - SWIFT il MODEL 65 x 45 130 -- SWIFT 11 MODEL 41 x 28 o -
4 _ m _ I -~ g l
. h,120 m
.**s ..*.N
_--,N s /
/
. * *N l iiO
. h...**..... 4
*.-.c.....*K ...,
% J I
'00 --
l l I I l l 1 800 1000 1200 1400 1600 I O 200 400 600 DISTANCE (m) Figure 3.2.14 Hydraulic-head profile for Z=-200, Case 2a-SWIFT II.
HYDROCOIN LEVEL 1 - CASE 2a HYDRAULIC HEAD DISTRIBUTION AT Z=-400m l l l l l l l 140 -
~
GWHRT MODEL ~
***** SWIFT 2 itODEL 81 x 56 ~
n
==== SWIFT 2 MODEL 65 x 45 -
E iso -
- - SWIFT 2 MODEL 41 x 28 _
o - i m -
- I -
w o - 4 e 3 120 2
- ___ ~~
tr
% N --~ ~ ~ ~ ,~ ~ ~ ~ --
l " ' '
-Q' ~~ _ ,~~
- _ - Q, .__: .
l tw -
-=............. Q _ # "'f _
i _ 100 - 0 200 400 600 800 1000 1200 1400 1600 DISTANCE (m) i I i i ! Figure 3.2.15 Hydraulic-head profile for Z=-400, Case 2a-SWIFT II. I J
i i HYDROCOIN LEVEL 1 - CASE 2a HYDRAULIC HEAD DISTRIBUTION AT Z=-600m i I I I I I i 140 -
~
GWHRT MODEL [ eeoe** SWIFT 11 MODEL 81 x 56 _ m
- --- SWIFT 11 MODEL 65 x 45 w
E 130 -
-- - SWIFT 11 MODEL 41 x 28 g
un _ I _ O.- Y g 120 - y < cc o _ ~ w g. _~~p#"" _ y I i 110 -
...............'....... ~~
~
I I I I I I I O 200 400 600 800 1000 1200 1400 1600 DISTANCE (m) l Figure 3.2.16 Hydraulic-head profile for Z=-600, Case 2a-SWIFT II.
HYDROCOIN LEVEL 1 - CASE 2a HYDRAULIC HEAD DISTRIBUTION AT Z=-800m t I I I I I I l 140 - -
~ ~
GWHRT MODEL
,, ...... SWIFT 11 MOuEL 81 x 56 _
SWIFT II MODEL 65 x 45
= = = -
w E 130 -
-- SWIFT 11 MODEL 41 x 28
-~
o . us w I - - l 4 e 0 120 - _
- 3
--- - - - . . , '_'.::: m g --~~~ - .
~
- ................----.....,=-'...---.,,,-
E _ o _ I 110 - ---...................a=***,, i i _ _ l l 100 - - I I I I I I I 0 200 400 soo 800 1000 1200 1400 1800 DISTANCE (m) Figure 3.2.17 Hydraulic-head profile for Z=-800, Case 2a-5WIFT II.
~,,,...m ac"
-~cceee
\ -
___J, , k(h 4 [f r > 3 1 jt g
% r 1 ___%
j 1
/
/
/.___.__.-k -
l, /
/ ".
% J %. i E .-w. ,- -
-g ,
.___ 2 3< 4 .,
.l J h 'd-g ; r
% 's 1 i 'r I.
k-s (-. 14 4. r L-_--.Qk j x e ,- y N
- j 4 g ja
' s- '
y e ( j ff y .9- ----
, a
't __. f .
j
;- i
%-- y L f I'G /
I (' '.N /
--_'a r >
77 1 i i I '3't3 A f l i f 'Q /~ i \
;[I j D(s i J g
'$ W
- 11 '? E n
r
--f w ri '
B. j f' 'D
'd l{ %
w j . Si _ _ _ . _ _ _ . _ _ _-_- ---- ------ - '~~- Figure 3.2.18 Pathlines for fine grid, Case 2-SWIFT II.
l l I I I I n v w % A N A // U \ \ // V N \ #
\ \ \ //
't h 2 9 1 \ \ // ---~
1 1 h // l 1 //
\
9
\' . //
s } // . m 3- I . // [
~c' g s- 1 // ,___ g
~~Cr %
t ' I // w k I C']~_
~~ ~ J ~
~a_
///
.r/
3 m
+ R1 *
// ~~
%-ct
~ - - m ,
l // N-22 + 4 e i .9 ~'~ ~ : 0 h
!/ b
// 'C i D c
// E
// o
!/ kW l l t //
i l //'
! ff h
! // -
\_ *
// ! +
l l ! // I i l // l
! O 1
\ i.// "g
\ Y/J '% 0 i //) \
\ DJ \ m i //J \ 1 l BI i \ N 1 h// h. V d
/47 A i o rp iA L o
! h , c I ! L 1 w I t I I I
I ! l .1 / i I 1 / 1 / / I f f ~~~ V /" / 1 I / / I r
/ / '-
l n f 1 I I . i l ( l 3-30
x -
/
/
~
)
l
/ A A.
' h '
- ^ / A '
A ' (~
/ A A
lI _A f
;IT l l A ' I I
f T _ I) ! HF .h'A' F
)l j I
_ l fir " l / W S
;/ l
- .l /[
l l 2 h) f l )
/
/ / /
/ e
! l
/ s N / / a W I /
/
/
/.
C W i d K r g N a, W e s gN r a o
, Nx c r
NN f o N s h \. e n i
- l h
L %P / ,/ t a
- P LN / 0
. LN 2
- LNi / 2 3 7 N /
' e r
q
/ l I'
i F u g
% .\
N N'
\ \
i N \ y"
three SWIFT Il simulations and the finite-element simulation performed by KEMAKTA consultants using the GWHRT model. As indicated by these figures, the results of all medels are very similar. Note also in the three SWIFT II simula-tions the decreasing hydraulic heads with each finer discretization. Also, at several locations, SWIFT II results indicate lower hydraulic heads than those of the GWHRT model. This is a result of the way each model treats the fractures. With the finite-element model, elements that correspond exactly with the frac-tures were constructed. However, with SWIFT II, large rectangular blocks are t used to represent different segments of the fracture. Thus, in the SWIFT II simulations, the effective width of the fracture is much greater than its real width. This in turn causes a greater reduction in hydraulic head along the l fractures. In addition to the hydraulic-head profiles associated with the grids in l Figures 3.2.4, 3.2.6, and 3.2.7, results are also given for the grid shown in ' Figure 3.2.5 (see Figs. 3.2.13-3.2.17). This is the mesh where fracture G (on Fig. 3.2.1) was represented by blocks lying along a diagonal that were not directly connected to each other. Blocks representing the ren aining fracture, however, are connected. The overall effect as shown in Figures 3.2.13-3.2.17 is that the medium grid results lie outside both the coarse and the fine grid results near fracture G. Near the other fracture, the results from this grid lie between the coarse and fine grid results. This discrepancy is caused by the lack of ir.terconnection of fracture G, which results in less hydraulic-head loss along the fracture. The remaining results are the pathlines as specified in Figure 3.2.3. These are displayed in Figures 3.2.18-3.2.20. Shown in Figure 3.2.21 are the pathlines predicted by the GWHRT model. In general, the location of the paths are in agreement with the GWHRT model results; that is, with the exception of pathline 2 on the medium and coarse grids. Pathline 2 ends up at the intersection of the surface and fracture G for the medium and coarse grids, but it ends at the intersection of the surface and fracture E for the fine grid and for the GWHRT model. The reason for this is found in the nature of the model predicted hydraulic-head surface in the region where pathline 2 intersects fracture G. This surface indicates a ground-water divide in this region. For the medium and coarse grids, the resolution of the model results is such that the pathline falls above the divide, but the higher the resolution of the fine grid reveals the divide to lie more toward the top of the model. This, in turn, causes the particle to travel down fracture G to fracture E and then to the surface. Another curious feature of the pathlines is their jagged appearance near and in the fracture block. This effect is caused by the way the finite-difference 4 model connects the fracture blocks and is discussed in Appendix B. Overall l lengths of the pathlines and their associated travel times are listed in i l
- Table 3.2.1. With the exception of pathline 2, which was discussed previously, l the results generally agree with the results from the other codes.
One final point is that even though many more nodes were used to simulate Case 2 t with SWIFT II than are necessary with a finite-element model, the actual com- ! puter time may be less because of the efficiency of finite-difference matrix l solvers. For instance, a comparison was made between the SWIFT II coarse grid i for Case 2 (1148 nodes) and a Case 2 coarse grid used for the SAGURO finite-element code (416 nodes) on the same computer, the CRAY-1. For this problem, the SAGURO code used 10 seconds of CPU time while SWIFT II used 8.5 seconds. l t 3-32
.,y-- .,.,,.p- -
r -.,-,---.r .w,------- - - - - ,-r--v---
200 , , , , , , , , , , , , , , , 100
' "' "'~5 1
3 0 - N /--
-100 N
g -200 2 \ -
- N g -300 -
\
w 0 O O
-400 -
\ 4 w
s O N o -500 - - t 8 x -600 - -
-Too - -
/ \
-800 N
-900 -
\ -
I I I I I I I I I L I I I I I Nl _,ogo o 200 400 soo 800 1000 1200 1400 1800 100 300 500 700 900 1100 1300 1500 x - COORDINATE Figure 3.2.21 Pathlines for Case 2 from GWHRT model.
i a Table 3.2.1 Travel times and distances for pathlines .; predicted for Case 2-5WIFT II. PATH 1 PATH 2 PATH 3 PATH 4 .1 GRID Distance
- Time ** Distance Time Distance Time -Distance Time 'l C0 ARSE 572 1,037 1,099 5,388 469 709 1,197 9,245 MEDIUM 604 959 1,068 5,602 587 706 1,234 7,518 618 1,065 1,757 12,590 631 754 1,394 13,605 i FINE
- Distance is in meters.
! ** Time is in years. 4 OJ 3 I i
- v wn~ a w -
USGS The construction of a model for Case 2 presented special problems to the USGS and the SWIFT II codes. The problem required that two intersecting fractures, which are narrow relative to the scale of the problem and likely to largely control its hydrologic characteristics, receive the most detailed representation possible. The representation for both codes was obtained by maximizing the correspondence between the fractures and the location of the block centers. For the USGS code, however, three additional steps were taken:
- 1. The grid (Fig. 3.2.22) was rotated relative to problem coordinates so that a column of block centers bisected the angle between the intersecting fractures.
- 2. One block was located at the intersection of the two fractures, and the grid spacing was defined so that the fractures corresponded to diagonal lines in the grid.
- 3. The hydraulic parameters in the model were adjusted to provide a correction for flow along the diagonal lines representing the fractures and to ensure that the harmonic mean transmissivity calculated by the code produced the hydraulic parameters :;pecified in the case description.
The correction is derived in Appendix B and simply summarized here: ; K = wff K (Dx2 + Dy2 )) K K,(Dx2 + Dy2 )/(4K,0yDx - wff 4 where Dx = block dimension measured along a row, Dy = block dimension measured along a column, K = corrected fracture conductivity value (one value for each fracture), K f
= fracture conductivity, K, = matrix conductivity, and w = fracture width, f
i The problem is represented with a grid consisting of one layer, with 29 rows and 65 columns, in which 1195 blocks were active. The block dimension measured along rows is 28.7710 meters. The block dimension nieasured along columns is i 54.9115 meters. The block at row 20, column 40, corresponds to the location of the fracture intersection. The boundaries of the active grid are located as near as possible to the actual boundaries of the problem. Blocks representing the held potential boundary at the top of the profile are assigned potentials equal to the elevation of the model boundary at its closest approach to the block centers. Results of the simulation are presented in terms of hydraulic-head profiles (Figs. 3.2.23-3.2.27), hydraulic-head contours (Fig. 3.2.28), and model predicted pathlines (Fig. 3.2.29 and Table 3.2.2). All the results compare well with the finest grid used for SWIF1 II. 3-35
>=
k a I ' e E cn tS w w a M a c; t a k 1 3-36
i 1 I ) 129 - g i I
~
i I , 126 - OWHRT MOOR ~ _ m*3-30 MODEL / i 123 - f - i j a E
/
! 12e - i O [ - i < W
- 117 -
y 9 4 w y 114 - > " < / ,t 1 C 111 [ - t j r 3 tos - 3 l - ! \ / 105 - ! \ /
- 102 -
j I f i 1 1 I i 1 j o 200 400 soo 800 1000 1200 1400 1800 j DISTANCE (m) f l l i j 1 Figure 3.2.23 Hydraulic-head distribution at Z=0m for Case 2-USGS. l l
0 0 6 _ # - - - _ 1 7 . S 0 G S 1 / I 0 4 U L L [ 1 2 E E D / i e s D O a OM C M 0 0 T 3- 0 r R o H S I l 2 G 1 f WS m G U 0 0 2
- =
0 Z l
- I 0 0
1 t a
) ,
( m i c E t
- u 0 C b l
I 0 N i 8 A r T t S s I i D d d a 0 e g I 0 h _ 6 - c _ i l u a _ x d r i x I 0 0 4 H y
\ 4
\ 2 2
3 0 g I 0 e _ 2 r _ u _ g i F _ - - - - - O 0 , 6 4 2 0 s , 2 , 1 1 1 1 o i 1 , 1 1 1 1 5 ob12d<g>r I-
. ;! , I j l:I! !
i I 116 g g g l l l l 1 115 -
\ _
GWWIT MODEL
\
-- usos-30 MODEL 114 -
\ -
3 \ o 113
\ ~
$ \
Y 52 112 -
\
w 8 N. " Q N / e N. o sii - g / _ I N 110 -
\ ~
\
i 109 - _ l l 1 I I I I I O 200 400 600 800 1000 1200 1400 1600 DISTANCE (m) 4 Figure 3.2.25 Hydraulic-head distribution at Z=-400m for Case 2-USGS.
0 0
- 6
- - - 1 a
- S L
f G L S E E D ' 0 U D O i 0 - I OM 4 1 2 M 0 e T 3- s mSG C a WS r G U 3 o I = I C f
=
2 1 m
= 0 0
- - 6
- =
Z 0 I 0 t 1 0 1 a n
) o i
N ( m t N u E b
% 0 C i r
1 \ I 0 8 N t N A s T i N S d I N D d a N e
% 0 h 1 0 l
N 6 i c N l u N a N r d N y 0 H l N 1 0 4 N 6 g 2 _ N 2 g 3
\ 0 e l
0 r I
% 2 u g
i F
- - - - - O # 3 2 1 0 9 s 0 " 1 1
1 1 1 1 1 1 0 1 1 E* O5r QjOmo$ = wi
' ; !,.! ! i ;{! j ,
1 I I I I I i 113.5 - - l GWHRT MODEL 113 >==. == = = = USGS-3 D MODEL _
's 112.5 - %
g - O '% g
< s S
DJ 112 - g - I sg O
- s l w
] 111.5 3
's %
i E < 's I #
$ 111 -
's -
>- 's l %
i 110.5 - S g - s 110 -
'g -
s 109.5 - s...,_ ... - - io, I I I I I I I O 200 400 600 800 1000 1200 1400 1600 DISTANCE (m) Figure 3.2.27 Hydraulic-head distribution at Z=-800m for Case 2-USGS.
A
?~ -
))
J
>/
(o~ 4' Y)! f m
,1 n
m; L
, (
,m ,l I
Figure 3.2.28 Hydraulic-5ead contours for Case 2-USGS. 3-42
s 1
\
i (
\
I k 1
\
\
\ '
\
\
\
A
\.. -
( Figure 3.2.29 Model predicted pathlines for Case 2-USGS. 3-43
Table 3.2.2 Travel times and distances for pathlines predicted for Case 2-USGS. PATH 1 PATH 2 PATH 3 PATH 4 Distance
- Time ^^ Distance Time Distance Time Distance Time 424 677 1,680 12,929 659 947 1,232 12,189
- Distance is in meters.
** Time is in years.
Y t
i FEWATER Two grids were used to simulate this problem with the FEWATER code (Figs. 3.2.30 and 3.2.31). These meshes contain 112 and 378 elements, respectively. The finite-element code represents the geometry of the fractures more accurately than with the finite-difference codes (SWIFT II and USGS). In addition, fewer total calculational points are needed to represent the whole system with FEWATER (112 and 378 elements) than with SWIFT II (1148, 2925, and 4136 grid blocks) or . the USGS code (1885 grid blocks). However, this does not necessarily mean that ' less time is required to set up and run the FEWATER simulation. Results of the FEWATER solutions are presented in Figs. 3.2.32-3.2.38 and in Table 3.2.3. i. i The head distributions for the two grids show an expected increase in deviation between the two with an increase in depth. The resolution of head values in regions with a large head gradient will be significantly affected by the dis-t cretization density (see Figs. 3.2.30 and 3.2.31 for grid density). Head gra-l dients in the low conductivity matrix will be larger than gradients in the high conductivity fracture zones. Therefore, differences in head values for the two ' i discretizations should decrease near the fracture zones. The head values for i the two grids show good agreement near the wider fracture zone (see Figs. 3.2.32 ! and 3.2.34). The reason for the largest deviation occurring at the lowest depth 1 and the smallest deviation occurrir.g at the highest depth is that the head values . ! on the surface are fixed by the imposed boundary condition. Head distributions i i for both discretizations are influenced by the fracture zones and currently j show a lower head value over the wider fracture zone. Particle trajectories, as shown in Figures 3.2.37 and 3.2.38, were produced ; using a post processor with the FEWATER output. The post processor ettimates a particle trajectory through a finite-element grid by the following process:
- (1) determine the specific element the particle is in, (2) linearly interpolate the velocity at the particle location from the four nodal values of the element, '
(3) determine the time increment needed to move the particle a distance equal , to 2 percent of the smallest side of the element, (4) move the particle accord- . l ing to the previously calculated velocity and time step, and (5) print the loca- ! tion at every fifth time step. The particle paths can be further analyzed by ! comparing the overall length and travel time for each path (see Table 3.2.3). l Particle paths 1, 3, and 4 generally compare with t.e results from the other codes. For the coarse grid, the particle trajectory for path 2 never exited ! the grid and remains trapped within the middle area of the grid. This result is due to the numerical inaccuracies introduced by fhe use of too coarse a dis- l cretization. A correct pathline was ;btained for tie finer discretization ; (see Fig. 3.2.38). ; l 3.3 Case 3 Problem Description t 1 The Case 3 problem involves unsaturated flow ;n layered sedimentary rock. The f different layers have highly contrasting hydraulic conductivities such that the grid network in the vicinity of the contrasts becomes very important. The shape , and boundary conditions of the modeled domain are shown in Figure 3.3.1. ; 3-45 l t
o I- i. : 4, e
/ o i 8 3
l N f o A
% /~ i
% ~
N
- / ! o g
$5
, / e ,
/ ,5.
g
- / 8<I e !
1 g u O o 4 o '
- a. o x Y.
- =
! 8 u
- rcn
! / e i . i e a 8< 1
> e M I i A l l M i e .
l $ 1 + i .i o o o o o 0
- d i i 8 8 i
i
? ? ? ?
j (eJelow) 31VN10WOOO - Z 4 j 3-46 a y._,y_. -,m_.- ,-, - - ..----,.,__---,,ym,-,,y,
- , - -,e w r--m--e+ewmm+--+ev ---,nw --+emww-em er w e--~ ey-m r- t w , -v--e e ewwee f c- wvw- --w -----=v
o . 8 l ,
- //
~
~ ~ /
/
/ //
/
//
h :
/ /
- -- %~ /
j:c :~ ::::: d' x
' ' // //
/ := :=
o jd/- ' ~ 5
/ ::::: :~- -
Y :::::-
~~ ~
- ' ' /
y _
i? 3
/
/
/ / / ~~~'
' ----*i
. :. i
/
/,/ h,// ' -
- 5 i
/// // "
i kl
/// / // -
ti !
/ o
#//// -
/
/
/
/
/
/
m T [ / .h
/ / ' '
8 '
/ , /
/
F o o o o o o S ~ 8 8 8 8
? ? ? ?
(sJenew) 31YN10HOOO - Z 3-47
132.0 g g l l l 1 I FNE ORIO
* * * *
- COARSE GRID 12&O 124.0 -
g 120.0 - 5 O
- y 116.0 - -
1 : J . 4' . H . O . H 112.0 - 1o8.0b . I 104.0 - - 3oo o i I I l l l l i 0.0 200.0 400.0 600.0 800.0 1000.0 1200.0 I t 00.0 1600.0 X - COORDINATE (meters) Figure 3.2.32 Hydraulic-head distribution at Z=0m for Case 2-FEKWATER. 3-48 _ _ . ..-.-..__.-.,.m_ _ . , . - , , _ - , , . . - _ , - . - 3 - - - - - ,
l i 132.0 g g l l l l l FWEGRO 128.0 -
. . . . . . COARSE ORID -
124.0 - - q 120.0 - _ w 8 o < 116.0 - W lI: H O . H 112.0 - - 108.0 - 104 0 - - I I I ! ! I I 100.0 0.0 200.0 400.0 600.0 800.0 1000.0 1200.0 1400.0 1600.0 X - COORDINATE (meters) Figure 3.2.33 Hydraulic-head distribution at Z=-200m for Case 2-FEffdATER. 3-49
1
'50 l I I I I I I 114.0 - * " "U -
., . . . . . . COARSE ORIO 113.0 -
e 9 112.0 '. j e . O *
< 111.0 W
I J H .'. O p 110.0 - ..**...' ,
. e 109.0 -
- e, .
108.0 - - 1 I I I I I I 107.0 O.0 200.0 400.0 600.0 o00.0 1000.0 1200.0 1400.0 1600.0 X - COORDIN ATE (meters) Figure 3.2.34 Hydrault :-head distribution at Z=-400m for Case 2-FEHWATER, 3-50 4
115.0 g y l l l l l FINE GRID 114.0 -
* * . * *
- CO ARSE GRID -
113.0 -
+..,*.'. ,g 112.0 - .4 .
2 . 5 . o *
< 111.0 - .
E . a . ( *.
- O
- r 110.0 -
., ,..***J 100.0 -
104.0 - - I ! ! ! ! ! I 107.0 0.0 200.0 400.0 600.0 600.0 1000.0 1200.0 1400.0 1600.0 X - COORDINATE (meters) l l Figure 3.2.35 Hydraulic-head distribution at Z=-600m for Case 2-FEMWATER. 3-51
115.0 g l l l l l l FINE ORIO 114.0 -
. . . . .
- COARSE ORID -
113.0 - i..**
- iiz0 - .. _
.a .,, ..* g w o ...** 111.0 - z ., J
- O s 110.0 - ' '. , -
10s.o - 10&O - -
! I ! ! ! ! !
107.0 0.0 200.0 400.0 600.0 800.0 1000.0 1200.0 1400.0 1600.0 X - COORDINATE (meters) Figure 3.2.36 Hydraulic-head distribution at Z=-800m for Case 2-FEMWATER. 3-52
\ \ \ \
\ \
s s ' Q
\ N N s N I
- #s'/ --
a
$g
______. e
$S 3 E -
o
. a
, q .
\ ~
Q
\ o
% v O
\
b 8 y - e" 1
\ s
\ :
/^\
1 m i i g l / i 5 I Q l / -
/ /
/ /
\ \ \ \ 1
, o p , o o
8 e %, 1 , (,,,,,, st * * *~* 3-53
e N g N s o N \
% \
.i \
\ /
\ ,
5 g _ _- _/ e 2 1 ______ E k a m g, e g a
# C e
s ';
; e q7 E
'T gg ;
a s
\ o
- 0
- 5
\
1* 1
,S g
/m\ \
M a s \ - 1 a
\ I e
\ s i / C ,
I / -
/ / 1 l /
l /
/
\ \ \ \ \
Q o O Q O o O e $t b i ' '
,,,,,m uw$u ~*
3-54
Table 3.2.3 Travel times anc' distances for pathlines predicted for Case 2-FEMWATER. PATH 1 PATH 2 PATH 3 PATH 4 GRID Distance
- Time ** Distance Time Distance Time Distance Time a
FINE 690 980 2,000 13,000 680 790' 1,340 7,980 COARSE 730 760 *** *** 690 890 1,340 7,930
- Distance is in meters.
l ** Time is in years. ! *** Particle never left the modeled region. I T i E ) i l 1 i I 1 4 'l I
q = 100 mm/yr F t Y I I I Y J
+500 m A B q=0 f
C E O - l q=0 0 h=0 -
-500 m O q=0 2020 m G H ASSUMPTIONS
- All media (A-E) are considered homogeneous and isotropic
- Media A and C are permeable compared to media B and D
- Medium E is much more permeable than media A-D
- Flow through the system is considered to be saturated-unsaturated
- Initially hydrostatic equilibrium
- Infiltration along F-J (100 mm/yr)
- No-flow boundaries along F-G and G-H
- Hydrostatic pressure along I-H
- Seepage face along I-J l
- No precipitation along I-J I REQUIRED OUTPUT
- Location of ground-water table at successive time steps and at steady state Figure 3.3.1 Conceptual model, assumptions, and specification of output for Case 3.
3-56
In this problem a two-dimensional vertical section with five media (denoted A-E in Fig. 3.3.1) is to be modeled. Medium D is assumed to be the host rock for shallow land disposal. The media A and C are permeable relative to the host rock, whereas medium B is to be considered a confining bed between media A and C. Medium E is a very permeable layer of loose material on the right flank of the model. Saturated and unsaturated flow conditions are to be simulated. The system is initially in hydrostatic equilibrium. At a given point in time, a specified flux is imposed across the top surface (medium A) of the modeled domain. The development of the ground-water table is followed as a function of time. The imposed ground-water recharge rate is 100 mm/yr. Hydrostatic pressure is imposed along the boundary I-H. All other boundaries except J-E are assumed to be no-flow. Boundary J-E is specified as a seepage face. All media (A-E) are considered homogeneous and isotropic. The required output for this problem consists of the location of the water table at steady-state conditions and at five times that are evenly spaced between the start time and steady-state conditions. Simulations Performed This problem was simulated with the finite element code, UNSAT2. The model grid is shown in Figure 3.3.2. In the horizontal direction, grid spacing ranged between 10 meters on the right boundary and 50 meters on the left boundary. At points where grid spacing was expanded or contracted, the ratio of one block length to the next was kept at 1.35 or less. In the vertical direction, grid spacing ranged between 7 meters and 50 meters. The larger grid spacing was used in the 0 rock group, below the water table. This part of the modeled area would stay saturated during the simulation 50 that the fine gridding necessary for nonlinear unsaturated flow problems was not necessary in that area. All grid blocks had a rectangular shape except those just below the E rock unit where triangles were used. Simulation Results Two types of simulations were made. The first type of simulation attempted to progress from initial conditions to steady-state conditions by means of a tran-sient solution. The solution would progress normally through time, but after several time steps the solution would not converge to within the specified tolerance. Usually decreasing the time step solved the problem, but after a while the time steps were so small that the solution to the problem hardly changed at all between time steps. It then became a trial and error method to try to find a proper time step. Computer runs with the smaller time steps with UNSAT2 were becoming increasingly expensive, and thus transient solutions were stopped. The second type of simulation was an attempt to obtain a direct steady-state solution so that an idea of the water rise could be obtained. This was done by setting both the specific storage and moisture capacity to zero. This simula-tion yielded steady-state conditions. However, upon examination of the output 3-57
s , E1 fr l
}
// l B
/
}/
R 1
! 4 D .
D A D ff f4 IJ fa e
/l N *
< l:
,,n m ,
l' E
/ H T m
-j W
; # l!
o 1 // !i
// D L- .
g In f f' I
? l T
I I 1 g i b
, L3 i
== t I I N !
I i i
.. M. ;
I l l U L .i II 3 j I y
.- i k 6 l'
- I !
j i ! o t 4 m [ ] I 4 I i . 1 !} 1 , 1 J ! , ", I _U_ l l_ o f i- e i. j 3-58 l , l
-- - - - _ _ , - , , - - - . - . . . . . - , - - _ - _ , . . - - ~._ , _ _ . ,_ - - - _
data, it was noticed that the entire modeled domain stayed saturated in nega-tive pressure head zones. The cause for this was based on setting the moisture capacity to zero. This was done by setting the allowable negative pressure head change to an extremely large value (10 50m) for moisture content changes between 0.00 and 0.03. The calculated pressure head changes were on the order of several hundred meters. When the model-calculated moisture contents were based on the pressure head calculations, the moisture contents always remained very near saturation. Recause the hydraulic conductivity characteristic curves used by UNSAT2 are based on moisture content, the unsaturated hydraulic con-ductivities remained equal to the saturated hydraulic conductivities. There-I fore, even though pressure heads were negative above the water table, conditions there remain saturated. In any event, this did not affect the model predicted location of the water table. The water-table location is shown in Figure 3.3.3. The water-table rise was approximately 67 meters at the left boundary and 0 meters at the right. An attempt was made to compare the UNSAT2 solution to a Dupuit-Forcheimer solution of tne water table. With conditions specified in the problem, the 1 Dupuit model predicted a water-table rise of about 9,858 meters, significantly different from the UNSAT2 results. This indicated that much of the total infiltration at the top of the rock was not reaching the water table. It was found that the flow rate of water in the E rock type along the slanted coundary j was slightly less than the total infiltration to the system. Therefore, it
- appears that the B rock type effectively prevents most infiltration from reach-ing the water table.
3.4 Case 4 Problem Description i Case 4 is a thermal transport problem, it involves a spherical heat source ) located in an infinitely large porous mediuni (Fig. 3.4.1). The heat source, which represents a high-level radioactive waste repository, has in initial power output of 10 megawatts, which decays exponentially with a half-life of 30 years. The radius of the heat source is 250 meters. The porous media are saturated with an isotropic, homogeneous permeability of 10 " m2 and a constant porosity of 0.0001. The properties ol water are assumed to be equal to the values that ! would occur at 40 C. ! Several types of output are required by HYOR0 COIN for this problem and they are to be compared to the analytical solutions available for this problem. First, , the temperature and pressure changes as a function of time between 1 and 10,000 years are required for elevations of 0 meters, 126 meters, 250 meters, and 375 meters directly above the center of the sphere. Second, the temperature and l pressure rise as a function of distance between 0 meters and 750 meters above l the center of the sphere are required for times of 50, 100, 500, and 1,000 years. Third, pathlines of particles leaving the m;dplane of the sphere at radial dis-tances of 0 meters, 125 meters, and 250 meters and at starting times of 100, 1,000, and 10,000 years are requested. The third type of output consists of the radial coordinate and travel time to the plane located 1,000 meters above the midplane of th( sphere. , t 3-59
l i i! I' .\l1ll l1 !l1 ,iI1 1i m w w r w w 7 - T _ T 3 e s N a C w r W o w f w t n i o _ e t - y . a - v e r y l - e w l e w b
- a t -
r v
% V e t v a w w t w
d * . e t c v _ i d e r p
*9 n
w w l e 7 d n
.- o y M e
_. r g . 3 i m ._ 3 3 4 w e - r g u i F U T o,o
- ! 1!!1 11 > .J, : I,l 1, l11 1
t i ; R E P O SIT O R Y... i f r=250 m ? ASSUMPTIONS : a
- Repository is idealized as a uniform spherical heat source with same physical properties as surrounding rock
- Power output decays exponentially with time
- - Rock-mass is assumed to be a homogeneous, isotropic and saturated permeable medium of infinite extent
- Dominant heat transfer mechanism is conduction through the rock rather than convective transfer due to flow of ground water '
- Flow transients arising from compressibility effects are neglected l - Time-dependence of flow field arises solely from time-dependent heat I output from repository
- Density of water varies linearly with temperature and is independent of i pressure
- Viscosity and volumetric expansion coefficient of water are assumed to be independent of temperature and pressure
- Flow can be described by Darcy's law l
- In she absence of heat source there is no flow Figure 3.4.1 Conceptual model and assumptions for Case 4.
3-61 '
Simulations Performed Both the SWIF1 II code and the TOUGH code were applied to this problem. Following is a discussion of both applications. SWIFT II The SWIFT II grid was designed so that the center of the grid blocks would coin-cide with required output points presented in the tables in the Case 4 problem description (Fig. 3.4.2). Within the sphere, the grid spacing was 25 meters. The center of the heat source grid block was located at the center of the sphere. This kept the center of each succeeding grid block at an even multiple of 25 meters above the center of the sphere. This pattern was kept up until 800 meters above and 250 meters below the midplane of the sphere. Beyond that, the grid spacing was expanded by a ratio of approximately 1.20. The upper and lower edges of the grid were located 3,982.5 meters above and 4,122.5 meters below the midplane of the grid, respectively. In the radial direction, the same 25-meter pattern was repeated until the edge of the spherical heat source was reached. Beyond the sphere, the grid was expanded at a ratio of approxi-mately 1.20 until a distance of 4,500 meters was reached. The upper, lover, and outer radial boundaries were all assigned constant pressures and tempera-tures at each node. Initial boundary pressures were based on hydrostatic conditions. Two simulations were made with the SWIFT II code. The first simulation used maximum time steps that were equal to one-tenth the time interval between the outputs required by HYDROCOIN. For instance, the HYDROCOIN problem descrip-tion for Case 4 indicates that output is required for times of 10 years and then 15 years. Therefore, between 10 and 15 years, the maximum time step could only be 0.5 year. This simulation was run for about 600 years. For the second run, SWIFT II was allowed to internally calculate the maximum time step based on r.umerical dispersion and overshoot criteria. However, time steps were restricted so that they would fall on times for which there was required output. The second simulation was run for 40,000 years and ran to completion in one-third the computer time required for the first problem. Plots of temperature rise versus time for various elevations above the midpoint of the sphere and temperature rise versus elevatien above the midpoint of the sphere for various times are presented on Figures 3.4.3 and 3.4.4, respectively. For both figures, the temperature rise predicted by SWIFT II compares well with the analytical solution. Of interest to note is that heat transport in this problem is governed more by heat conduction than by convection and dispersion. i Similar plots of pressure rise versus time and pressure rise versus distance are presented on Figures 3.4.5 and 3.4.6, respectively. A good comparison between pressure rise predicted by SWIFT II and the analytical solution is ; indicated cn Figure 3.4.5. However, SWIFT II seems to slightly underpredict pressure rise at pressure peaks. At later times, greater than 800 years, SWIFT II seems to slightly overpredict pressures calculated from the analytical solution. On Figure 3.4.6, the SWIFT I4 pressure rise predictions compare very , well with those calculated from the analytical solution for times of 50, 100 3 l and 500 years. The pressure rise for these times was computed from the first SWIFT Il run, which used the small time steps. For 1,000 yeat , SWIFT II l 3-62
1 EOCK O 12.8m,11 KOCKS C tem 3.-3442.5m p= 3.4062 a 107 Pe T= 40*C 16 20 26 26 27 24 29 30 1 J 2 NTIAL CONDITION . 3 p5992 T= 40 *C 4 i 5 : e 7 l j e ii i i e i b=0 Or e 10 I l j t i p m pgz (hydrostatic) ' I T = 4 0 *C OT p= 0 , 40 BLOCKS
@ 25 m /
PER BLOCK l h
@!!t ..
! l ! .
J - 60 ' ' t i 70 '!!! l ! ' I I 71 !I!! ! l !
. 72 1111 1 73 l l lllN.
74 hll l .l 75 i j 77 l t i 78 l L i 79 l ( l a s 4122.Sm p=3.6479a10 7 Pe T= 40*C y E Figure 3.4.2 Grid used to simulate Case 4-5WIFT II. t t 3-63 ! _ _ . ~ . _ _ _ _ . _ _ _ , , __ _ ,_. ._ . - _ _
HYDROCOIN LEVEL 1-CASE 4 TRANSIENT THERM AL CONVECTION IN A SATURATED PERMEABLE MEDIUM 100 .
.......i .
- 90 ANALYTICAL SOLUTK)MS _
M ---- AT Z = 0 m
~
80 * * * - *
- AT Z = 12 5 m J7kg -
- - AT Z = 2 50 m - g E 70 -
AT Z = 375 m .I - E o SwwTu / Q* 'o g 60 - SOLUTIONS / . g g O .s a 50 - I Q.D - Y F 3 -\ T 4 40 - f D\ - f
~
30 - cy D% - L W 2 20 - - I- 10 - c@j ' - O - 0 10 30 100 300 1000 3000 10,000 TIME (years) l l l l Figure 3.4.3 Temperature rise vs. time, r=0 meters, for Case 4-SWIFT II.
HYDROCOIN LEVEL 1-CASE 4 TRANSIENT THERMAL CONVECTION IN A SATURATED PERMEABLE MEDIUM 100 i i i i i i i I ^ 90 - - ! ANALYTICAL SOLUTIONS m 80 ( ---- AT T = 0 YRS - 1 v3 h%o.a,b. . . . . . . AT T = 50 vRS 7g _ q. _ g -- AT T = 100 YRS 60 S W 0000 0 ~ O SWFT N SOLUTION sb g Y D 50 - Q~~ - S F %
< k, -
cc W 4 30 O < ro % ~o.%
- D -
1 2 20 -
% 9,.D -
W l- 10 - - O ' ' ' " O 75 150 225 300 375 450 525 600 675 750 DISTANCE (meters) Figure 3.4.4 Temperature rise vs. dist.aig.e. r=0 meters, for Case 4-SWIF1 II.
HYDROCOIN LEVEL 1-CASE 4 TRANSIENT THERMAL CONVECTION IN A SATURATED PERMEABLE MEDIUM 20,000 . 18,000 - ANALYTICAL SOLUTIONS -
^ 't "U ---- AT z = 0 m g 16,000 -
\ --- -- Ar z - 12s m
, 14,000 - - - ATz=25om_ t l
@ k AT z = 375 m l
g 12,000 - p do \ o sm si - 1
/ P sotuTsoMs w 10,000 g
i - - r w . I b 8000 -
/ O u) P .- h y) 6000 -
g ,. > W O E 4000 - M * - 1 P 9
# 9*
2000 ? - , - O 9 ODD 0 0 O O 10 30 100 300 1000 3000 10,000
- TIME (years) 1 Figure 3.4.5 Pressure rise vs. time, r-0 meters, for Case 4-SWIFT II.
l
HYDROCOIN LEVEL 1-CASE 4 TRANSIENT THERMAL CONVECTION IN A SATURATED PERMEABLE MEDIUM 20,000 , , , , , , , , , 18,000 _ ANAtmAL SOwmN _ j D t7 ---- AT T= 50 YRS 16,000 a * * * * *
- AT T = 100 YRS f "
' Qa 14,000 -
f
.o e O G
- - AT T = 500 YRS s o.*
AT T = 1000 vnS
.*i G O SWIFT 88 SOLUTIONS s Q.
m w 12,000 E
*Of
.r O 6.
g a > 10,000 - dp a6 - 8000 - p B N h - E 6000 -
-o g b E -
a pr g g0 00000000 %. fp
. 4000 o0 -
2000 - 2 ' Ol I ' ' ' ' ' I I I O 75 150 225 300 375 450 525 600 675 750 DISTANCE (meters) l Figure 3.4.6 Pressure rise vs. distance, r=0 meters, for Case 4-SWIFT II.
noticeably overpredicts the analytical solution for heights between 50 meters and 550 meters above the midplane of the sphere. Outside that range, the com-parisun between the SWIFT II and analytical solution is good. It should be pointed out that the SWIFT II solution for 1,000 years is based on the second simulation, which used the large time steps. The poor comparison for the 1,000 year solution may be due to the larger time steps. Pathlines for start times of 100 years, 1,000 years, and 10,000 years are presented on Figures 3.4.7, 3.4.8 and 3.4.9, respectively. The shapes and positioning of these curves compare reasonably well with those presented in the HYOR0 COIN description for Case 4. SWIFT II comparisons with the analytical solution for travel times and radial coordinates for the plane 1,000 meters above the midplane of the sphere are poor. Results of the comparisons are presented in Table 3.4.1. For final travel radius, SWIFT II consistently overpredicts the trivel radius calculated from the analytical solution. The error ranges between 0 percent and 15 per-cent. The travel time comparisons are very poor. SWIFT II can either over-predict or underpredict travel times by as much as 44 percent. The large per-centage error may be due to too coarse space and/or time discretization, to the SWIFT II velocity calculations, or to the pathline generator routine that uses the SWIFT II velocities. TOUGH The TOUGH code was also used in an attempt to solve the Case 4 problem. A grid was set up similar to the one used with SWIFT II. The main difference between the two grids was that the grid block sizes for the TOUGH simulations were twice as large as those for SWIFT 11 simulations. This was done to reduce the number of grid blocks in the TOUGH simulatior.s to 500 blocks. The TOUGH code requires large amounts of computer memory and 500 to 700 blocks is about the maximum number of grid blocks the Sandia computers could handle while running TOUGH. Because TOUGH has its own built-in water density equations, it was necessary to use TOUGH to generate the initial pressure distribution. In 810 seconds of CRAY-1 computer time, TOUGH performed seven iterations before having operand range errors failing to calculate the pressure distribution. The cost of the job was 225 dollars. Based on this result, it was apparent that to use TOUGH to simulate the Case 4 problem for the times required by the HYDROCOIN output would cost possibly thousands of dollars even if there were only one iteration per time step. Therefore, it was decided that the TOUGH code would not be used to simulate this problem. 3.5 Case 5 A two-dimensional steady-state saturated ground-water flow system that has salt along its bottom boundary is the subject of this problem (Fig. 3.5.1). The density of the wate* is assumed to vary with salt content, with a fresh-water density of 1,000 kg/m3 and a salt water density of 1,200 kg/m 3. A onstant viscosity of 10 3 Pascals is assumed. In addition, the assumption is made that the salt concentration of water flowing into the system is 0 and that no verti-cal gradient of the salt concentration exists in the outflow region (i.e., dc/dz = 0). The type of output specified to be used for code comparison is given in Figure 3.5.2. 3-68
li1l11\\ l ,l l , li' ;' i l ;i,f ll l l I _ I T F I S W 4 e f , __ s a
)/ I- -
__ C f f __ o i e m l))7 - I i t _ _ t
. r a
_ t
- s J .
r
- f. a e
,f y 1 -
/f 0 s#p/
0 1 f _ i r
. . / o
- __ 1 f r,1 ~
s
/ fI e n
f i l v~ h t a _ _ y- N P
- ,~
_ __ n 7
/ -
4
. _ / 3
-. Ir e
I~I 2J' r _. iII u g iIIf f m m - I F A1I YmD S
l l ill, Il J lJ! ,I s I I I T F I W 5 _ 4 . e _ s a
- C -
f
- n. o e
m i
,! t t
r . a - t - 4 s - r , _. a e e y 0 0 0, I _ h W
- 1 M
r o
,\I- f ,
/ s
/r1I .
_ i e n
/}- s#
2 l h t a P - i1 2 ny J 8
- - ._ 4 .
i;I~ .
,I
/r
)/ J t ,
3 e J f - _ . r
.i1ili
- u
- 4_
/ .
g
._ _ i
-~
_ ._ F n I l
.i ,. 1o _
K .t 1 .i.I.,1 , { _ 1 L s il' 1i3 a - - A W p . Ii Yt - iL m e A x K, u M, 3 3 I i'
4 4 Duud man -_ _ 4 6.e e I l' l l W I e i 9 v l - 4 O k i l a 4 l l 4 l u i
, i lj e
, 9 L
i j, 9 l l U l >i 8 1 , k
. 8
- 1 I C eg l
i : L i O t w 1 . ; . 1 I I .C N g i
=
Q.
- j -. .
~i
! l i
. I I l l 1 li ! .l .
i i i i i, , i i M
-,-s-, ..
I l b I, i i , .- o
] ! l @
; l
. . - , a-
!, [
8 I I i f
! i 1 i ! I
_ _ _ ..__4.-- -. _ ..--.-_.-. . _ .._ ,~.dl. - .. I L4. .. -_.._ 3-71
Table 3.4.1 SWIFT II and HYDROCOIN pathline comparisons. Start Time Start Distance Travel time (years) Travel Radius (meters) HYDROCOIN SWIFT II % Difference HYDROCOIN SWIFT II % Difference 0 1.5135x10 3 1.5678x103 3.6 0 0 0 100 125 2.1135x104 1.8890x104 -10.6 1.1050x103 1.1510x103 4.2 100 250 2.4586x10* 3.3360x104 35.7 9.4772x103 1.0225x103 7.9 100 43.7 0 0 0 1000 0 2.2972x10 3 3.3015x103 38.9 3.2284x102 3.6944x102 14,4 1000 125 2.9057x10 3 4.0361x103 250 5.5446x10 3 7.4537x103 34.4 6.3361x102 7.2306x103 14.1 1000 6.6390x104 -15.2 0 0 0 10000 0 7.8324x104 125 8.0396x10* 6.7193x104 -16.4 1.3429x102 1.4114x102 4.6 y 10000
-20.7 2.6811x102 2.8187x102 S,1 R$ 10000 250 8.7346x10 4 6.9275x104
p=10000 Pa N (900.0) A - (0,0) p=0 HYDRAULIC CONDUCTIVITY -16'm/s POROSITY .2 DISPERSION LENGTH - 20m LONGITUDINAL l 2m TRANSVERSAL NO DIFFUSION 1 l F E D C (0, -300) c n.=1 (900.-300) l ASSUMPTIONS
- Porous sedimentary rock assumed to be homogeneous and isotropic - Steady-state conditions - Fluid density depending only on brine concentration - Dynamic viscosity is assumed to be independent of brine concentration - Triangular pressure distribution along top boundary - No intake of salt along top boundary - No-flow boundary for both salt and fresh water along vertical boundaries - No-flow boundary for water along bottom boundary - Ground water is assumed to be saturated with salt between E and D Figure 3.5.1 Conceptual model and assumptions for Case 5. ;
l 3-73 : t
1 i i P i l [ i t A B f X X X X X [ 1 23 4 5 . . - _ t l i 1 I ! F E D C I b REOUIRED OUTPUT : i
- PRESSURE DISTRIBUTION ALONG THE SEVEN HORIZONTAL LINES
- STEADY-STATE SALT CONCENTRATION DISTRIBUTION ALONG THE HORIZONTAL LINES
- TRACKING OF FIVE PATHLINES (1-5) WITH STARTING POSITIONS ACCORDING TO FIGURE ABOVE l
- VERTICAL VELOCITIES ALONG TOP BOUNDARY (A-8) l l
l Figure 3.5.2 Specifications for model output for Case 5. 3-74
Simulations performed SWIFT II Only the SWIFT II code was applied to this problem as it was the only code currently being used by the NRC project team that solves the coupled processes of giound-water flow and brine transport. Several difficulties were encountered in applying the SWIFT II code to this problem. The main one arises from the bottom boundary condition of a no flow with respect to water and a salt concentration of 1 along the center of the bottom. While the SWIFT !! code allows for the user to fix the salt concentra-tion with Dirichlet boundary conditions, the user is 31so forced to specify a pressure at the same block. Because the salt affects the flos, the pressure is not known before the solution is performed. Therefore, a different approach to fixing a concentration of 1 along the bottom boundary had to be found. The approach settled on was to employ the use of the salt-dissolution submodel. This was done by putting a group of blocks along the bottom of the grid that were saturated with salt and had a high enough salt dissolution rate to keep the water in those blocks saturated with salt. Adjacent blocks along the bottom were assigned a value of zero for the hydraulic conductivity. The resulting grid is shown in Figure 3.5.3. The hydraulic conductivity of the blocks repre-senting the salt was set to 10 15 m/s. This is due to a strict interpretation of the boundary conditions, which state that the bottom boundary is no-flow with respect to water. Another interpretation can be derived from a statement in the problem description which indicates that salt flow along the bottom boundary can occur by "lateral dispersion." This was simulated by changing the hydraulic conductivity of the blocks representing the salt to that of the remaining model (Fig. 3.5.4). This simulation is herein referred to as Case Sa, The next difficulty encountered is also related to the SWIFT II stipulation that if Dirichlet boundary conditions are to be used, pressure and concentra-l tion must be simultaneously prescribed. That is, along the top boundary, a salt concentration of zero was entered along with the required fixed pressure. Unfortunately, this is not consistent with the specified boundary condition of I de/dz = 0. The specified boundary condition would allow salt to leave the system via water moveme~t but not via diffusion. However, the model input boundary condition does not allow any salt to leave the system. In retrospect, I a more realistic approach may have been to use wells to maintain the pressure along the top of the bottom boundary. This would allow water of any salt con-centration to leave the system since the difference between the specified well pressure and the block pressure would be the only driving force. Af ter setting up the model grid and boundary conditions, steady-state simulations were attempted. However, the models failed to converge to a solution with l either the direct solution algorithm or the two-line successive overrelaxation l iterative procedure. Finally, long-term transient simulations were performed to arrive at steady-state conditions. To achieve this type of solution effi- , ciently requirtd a tradeoff between time-step size (and therefore computer cost) and solution ac curacy. Initially, a relatively large t'me-step size was used which, when comained with the grid spacing and coefficient of dispersion, violated the numerical criteria for solution stability. This simulation was run for almost 32,000 years. The resulting solution produced oscillating values 3-75
l lI 1. i i I ,k o - I ,A i = -
; g -
i% a i
~
~ -
i t s u/m - Q . S 30 3
~ . I 1 e
,. = s K a C
e t
~ l a
u
- m i
~
?
e -- s o _ m - - _ . - - t
. - 6 . .
d
. , . . E e
- 8. I '
s u II, 9 '
~
i _ = d
- . . K i r
. - ' ~ G
.- - ~'
S 1
- K 3 C
_ O 5 L 3 B e N r
. ?_ . -
- _ ~
I I O I u
. . - - g
- n 1I i nTUL f i
O S
~ ,
-- S A I _
~ ._.
" a D _
~
1yII i T L
~ _
A _ S 1' . i a a l
~
. A
-~ _
i e I;IT i
, ,IyIT 1
j I;IT A i a 1i ij;i;I' i 0 a = I!1 1I t1 I K Yw&
I -.e... _ .. . . ...-- .--. .-...-. ..- _ _ _ ..-. _ .--.,
.-.--p,.-- ---- -...-..
4r->M6-.esei amed>+ipe- 1p--g~1p + p-em aiesw-.mp esi>-.s>-,qm .-.p.-.ai s ea-.i 4p.-4PMi+4b ed.-4,4 ,a%i.. .-
.iy-,.-ii.eump.am,.
w .. _ . N w ._._._>._ K = 9.8E-6 m/s
.---, - , .... . .-... -.e. -.. .---_ . .
-4F-i>My-g-4e>% mme.4p+- W-4P-==umsb m>46aIa.a>-q>-4>.p-. -4>>-a>+-
..p.--p 4,,g.-4.-.H5.,. -4%.- , -
.--.> - <p -<-<--. MWb dp--Hl-*1P-1Fe1> J <>-i -q-am>+. --4 -gii.H>- +- u s -mm. eq>4k ,-4 ..-q e4b-1>M+. -4m->=-4P=d>+ 4p..-gi%> -ibembdh-4s==mdm.m=4pe "* - - M** 1Pi.m=4i=4>=-4b+>mip..-gi as4p.-qpegb- <.-i*. enb.-4h*'h+4.e4-- 4luumm HM-4P=H.-g.-=d> 1>m.i-b eel-mus.-abHp-=-1 -e4r- p-4emi- .e-1 iii1 1 1 il . 1 i . lli 2
- w ;I -
K=o s ALT ossson.uTION BLOCKS K = 9.8 E-6 m/s K=0 Figure 3.S.4 Grid used to simulate Case Sa. of salt concentrations. To verify that the solution was oscillating around the correct values, the time-step size was draJtically reduced. Then the simulation was run for a little over 300 years. These value: showed little oscillation and agreed with the values ftom the longer simulation. Model results are presented in terms of pressure pred'les (Figs. 3.5.5-3.5.7 for Case 5 and Figs. 3.5.10-3.5.12 for rase Sa), salt concentration profiles (Fig. 3.5.8 for Case 5 and Fig. 3.5.13 for u se Sa), and model predicted path-lines (Fig. 3.5.9 for Case 5 and Fig. 3.L 14 for Case SP . As indicated by the pressure profiles and the pathlines, Cases 5 and sa produce virtually identical results for ground-water flow. With respect to sait crncentrations, each pro-duces nearly the same shape of profiles and contours (figs. 3.5.15 and 3.5.16). However, concentrations predicted by Case 5a (where water flows through the salt) are almost exactly an order of magnitude higher than for Case 5. Given the predicted pressures of each and the predicted paths, the conclusion can be made that the salt concentration is not large enough to affect the flow in either case. On the other hand, the salt profiles and contours indicate that the flow does affect the salt distribution. That is, both cases are characterized by a shift-ing of salt concentrations in the direction of flow. Another point of interest is the upgradient movement of salt by diffusion seen in both salt profiles for Z = -300m. As no analytical solution for this case is available, the only statements that can be made about model accuracy are those about physical reasonableness. That is, pathlines and salt concentration profiles for both cases appear to be com-pletely consistent with the boundary conditions and the physics of flow and brine transport. 3.6 Case 6 Problem Description Case 6 involves a hypothetical repository site in bedded salt (Fig. 3.6.1). Given in Figure 3.6.1 are the assumptions used in the formulation of the problem. Plan views of the problem layers and their associated hydraulic con-ductivities are shown in Figure 3.6.2. The problem is designed to test a code's ability to simulate steady-state flow ; in a geometrically complex three-dimensional system characterized by anisotropy ' and sharp contrasts in hydraulic conductivity. Additionai features of this problem not found in the other Level 1 problems include areal recharge and recharge from and discharge to rivers. Model olutions to this problem are to be compared with one another for the outr' specified in Figure 3.6.3. I Simulations Performeu Two computer codes were used in simulating this prob'em, SWIFT !! and the USGS-30 code. The description of the application of these codes to this problem fo? lows. 3-78 HYDROCOIN LEVEL 1-CASE 5 SWIFT ll l 1 DISTRIBUTION OF PRESSURE ALONG SPECIFIED LINES l 1,200,000 . . i l Z =-50 METERS i 1,100,000 - - - Z =-100 METERS - ----%_ s__ ~~ l - 1,000,000 - --'--'---- l (O I Q. 900,000 - L1J 4 w g y 800,000 3 i 0) U) 700,000 - LIJ i E -
- n. 600,000 -
500,000 - 400,000 ' ' ' O 300 600 900 DISTANCE (meters) f i Figure 3.5.5 Pressure profiles along Z=-50m and -100m for Case 5. i l I i l 1 f HYDROCOIN LEVEL 1-CASE 5
- SWIFT 11 l DISTRIBUTION OF PRESSURE ALONG SPECIFIED LINES 2,200,000 i i i i i i i i i - - Z = -150 METERS
- 2,100,000 -
--- Z = -200 METERS - l Q 2,000,000 ------ -________'~~~'---'------- ; v a. i 1,900,000 - LLI i 5 K - o O 1,800,000 - 0)
- 0) 1,700,000 -
- g i E i
Q- 1,600,000 - ' ' I ' ' ' ' ' 1,400,000 O 300 600 900 DISTANCE (meters) Figure 3.5.6 Pressure profiles along Z=-150m and -200m for Case 5. )l i HYDROCOIN LEVEL 1-CASE 5 SWIFT 11 DISTRIBUTION OF PRESSURE ALONG SPECIFIED LINES 3,000,000 - - ' . i i i ' 2,900,000 i~~~~~~ '~~'~~'"~-~~'----~--.--.._._._.-.T ^G3 2'800,000 - ---- Z = -250 METERS Q. .--- Z = -300 METERS- " 2,700,000 - i W 4 w g - l s 3 2,600,000 - u) - - - - - - - ' m 2,500,000 .-..-..-..._.._.._..""-"'----"------- _ W
- E -
! Q. 2,400,000 - i i 2,300,000 - 1 i I i i 2'200'000 ' ' i
- O 300 600 900 DISTANCE (meters) i 1
! Figure 3.5.7 Pressure profiles along Z=-250m and -300m for Case 5. i 4 i HYDROCOIN LEVEL 1-CASE 5 SWIFT 11 j DISTRIBUTIONS OF SALT CONCENTRATION i 1
- O.0004 -
i Z = -50 METERS [\ j Z _ - - - Z = -100 METERS ,[ \ _ ! O ----- - Z = -150 METERS / l 5 0.0003 - -- Z = -200 METERS / - j ( ---- Z= -250 METERS f l E - - - Z = -300 METERS
- N,' s .~
! m F / $ z
- w O.0002 -
i 0 - / / , f .-/ z - / /- 7 '.: O * ' o 0.0001 - / ,/ / - / / . * ,- } l l .- / ... -, _ _ ll -.. f . / 4 -..< 1 / l l n . . . V. / . O 300 600 900 DISTANCE (meters) Figure 3.5.8 Salt concentration profiles for Case 5. 'an- - u- 'Munu . '9 k %.'%m'%.% % L, g t, g \, 't ' ' k L L k L T k L 1 'l % 1 L 1 % 1 1 % 'l s 'k,. .._-- kk i k \ 't \ 1. 1 't T \ L 1 'l % 1 1 \ \ \ L 1. 1 % , 'l --}. L h g 1 g i I t L tu \ 'k h 1 m 1 1 L 1 tU 1 1 1 0 1 1 1 L 1 < o 1 L < w 1 1 I i i i e ' 0 .C E a to CL T a U epe V L a e 0 D o I E J l Ch L I . I f m I 1 m I J A C D E } J G I E D ' I I J J 3 F F J F CD I I I ? *" T I V ..'\ J W V I [ ] 1 I I I g -}__II T A .V y J [ l I ng E 1 r s i n i r 1 1 E } f E ~f / / A . --, ; r I / / J f f f / f / F / / / / fn J r # J a #"j/ p en Y a n*K MP' l I 3-83 ! HYDROCOIN LEVEL 1-CASE Sa SWIFT 11 DISTRIBUTION OF PRESSURE ALONG SPECIFIED LINES j l 1,200,000 . . i . . i i 1,100,000 - - 4 --- w__s__ ~~'~ l ^ - ^--- a 1,000,000 i 1 O. v 900,000 - Z =-50 METERS - m W --- Z =-100 METERS E T 3 800,000 - - ! M j M 700,000 - - W ! E j n. 600,000 - - 500,000 - l
- I ' ' ' ' '
400,000 ' O 300 600 900 DISTANCE (meters) Figure 3.5.10 Pressure profiles for Z=-50m and -100m for Case Sa. .l l HYDROCOIN LEVEL 1-CASE Sa SWIFT 11 DISTRIBUTION OF PRESSURE ALONG SPECIFIED LINES 2,200,000 , , 2,100,000 - . . . . . . Z = -250 METERS _ - - Z =-200 METERS f _cc 2,000,000 _---~___,___'--S---w_____s- ! ca. 1 v 3 1,900,000 - -
- IM Y E
! g 3 1,800,000 - - i (D $ 1,700,000 - - ! E l o. 1,600,000 - _ ' ' ' ~ ~ ~ - - - - 1,500,000 - . . . . . . . 1,400,000 ' ' l ' ' I i i ! O 300 600 900 i DISTANCE (meters) 1 Figure 3.5.11 Pressure profiles for Z=-150m and -200m for Case Sa. l HYDROCOIN LEVEL 1-CASE Sa SWIFT 11 DISTRIBUTION OF PRESSURE ALONG SPECIFIED LINES 3,000,000 . . , . - i - . 1 2,900,000 - '-*-----w.--s._._.__.s,_,,,_,,_ { g 2,800,000 - -- - Z =-250 METERS O- --- Z =-3OO METERS
- 2,700,000 - -
uJ i Y E 3 2,600,000 - - ! M j m 2,500,000 _- - - - - m . . .. m..- .. ," - " S - -w.. g .. - - - ! E l n. 2,400,000 - - 1 ) 2,300,000 - - 1 l 2,200,000 ' ' ' ' ' ' ; O 300 600 900 DISTANCE (meters) l i { Figure 3.5.12 Pressure profiles along Z=-250m and -300m for Case Sa. l f I i HYDROCOIN LEVEL 1-CASE Sa SWIFT 11 DISTRIBUTIONS OF SALT CONCENTRATION i i l i 1 O.004 - Z = -50 METERS .O.g ! g _ --- Z = -100 METERS , [ \, , ! O ~~~"Z=-5 METERS j g - -- Z = -200 METERS / \ l I- 0.003 - _..- Z = -250 METERS --- Z = -300 METERS \
- 1 Y
5 i-z . / N.N ~ W O.002 l ! / ** / l 0 . / f ../ /.- l f i o 0.001 - ./ / .. f /.. j. s - ../ / , ,n . . r
- 16. :L? ,
o j O 300 600 900 l DISTANCE (meters) i l j Figure 3.5.13 Salt concentration profiles for Case Sa. i 1 ' ~ J Er ~ J g'r e n,, ,J - I J' J J P, / j sl72r / r 2 r r, r 2 r 7y s 2 fp s. 's jg , r .#" d s 'r s, - r a
- s r
#r - v a s." p , - . s W "s , - . r , , a a, m-
- S m "o - ~
e " s = d - s 's ~ ' - a r, C - ~ - " - ~ r d - - o - - - f s - - = e n i l h t a p d e t c i d e r 1IlIa_ ~ p .- l e d o - M = - 4 * - 1 - - - 5 - % . 3 % : ~ n e IIJ l je' s ~ ~ r u ,% % - w g 4 % -- % - i % % % F fiII ' s w ' s s - L' % - g w' s % 7 w x % N ' , ,5 L u' \x, % I1IIlT A . ,( k % k' k \ 1 3 1 3 1 s\ \_1Ls <L L1 _ w=cc oo VERTICAL CROSS-SECTION OF AQUlFER BRINE CONCENTRATION 111111 1111111. 111111 DEPENDENT VARIABLE RANGE MAP CHARACTER 1 22 0.000 - 0.158E-04 ll'111 22 0.158E 0.359E-04 1 1111111 222
- 0.359E 0.560E-04 3giigi. 22 0.5COE 0.761E-04 2 111111 222
- 0.761E 0.963E-04 1111111 222 0.963E 0.116E-03 3 111111 22 3*
0.116E 0.13 7 E-03 111111 222 333 0.13 7 E 0.157E-03 4 1111111 222 33* 0.157E 0.177E-03 111111 22 333 0.17 7E 0.197 E-03 5 1111111 222 333* 0.197E 0.217E-03 111111 222 333 0.217E 0.237E-03 6 11111 222 33
- 0.237E 0.257E 03 111111 222 33 0.2 57E 0.2 7 7 E-03 7 111111 222 333
- O.2 7 7E 0.2 97 E-03 111*1 222 33 0 297E 0.318E-03 8 111111 222 33
- 0.318E 0.338E-03 11111 222 33 0.3 38 E C.3 58E-03 9 111111 222 333 .
C.3 58E-C 3 - 0.3 7 8E-03 111111 222 33 4 0.3 7 8E 0.39 8 E-03 0 11111 222 33 4* 111111 2222 33 444 111111 2222 333 44+ 1111111 2222 33 444 1111111 2222 33 444 1111111 22222 33 444* 1111111 2222 333 4444 1111111 22222 33 444* 1111111 222 333 444 1111111 2222 333 444
- 1111111 222 333 444 1111111 2222 333 44 * '
1111111 2222 333 444 1110511 2222 3333 444
- 1111111 2222 3333 444 11111111 222 3333 444
- 11111111 22222 33333 4444 5 11111111 22222 3333 4444 5*
11111111 22222 3333 444 555 11111111 2222 3333 4444 55* 11111111 2222 3333 44444 5555 111151111 2222 33333 444444 555* 11111111 22222 33333 444444 55555 111111111 22222 33333 444444 55555* 111111111 2222 33333 44444 55555555 111111111 2222 33333 444444 5555555555* 11111111 222222 3333 444444 55555555555555 1111111 22222 3333 444444 555555555555555555
- 111111 22222 3333 4444 55555555555 111111 2222 3333 4444 55555
- 1111 222 333 444 5555 68366 66666666666666666 111 222 333 444 55 666 77777777777777 66666666666
- 1122 333 44 555 666 77 80888888 77777 6666666666 8808 77777 866666666
.......233 11 44 05 666 77 ...7. 88 .999999 ... .....9........7........6.*****.*** Figure 3.5.15 Salt concentration contours for Case 5. i l l ! 3-89 l VERTICAL CROSS-SECTION OF AQUlFER BRINE CONCENTRATION 1111 11111* 11111 11111
- DEPENDENT VARIABLE RANGE MAP CHARACTER 11111 2 0.000 -CL194E-03 11111 22* I CL194E CL392E-03 1 11111 2222 I 11111 222* l 0.3 92E 0.5 90E-03 111111 222 l CL590E CL788E-03 2 111111 22
- 0.788E 0.98 6 E-03 111111 22 CL986E 0.118E-02 3 111111 222
- 0.118E 0.138E-02 111111 22 3 l
0.13 8E 0.158E-02 4 11111 22 3* 0.158E 0.178E-02 111111 222 333 0.178E CL 198E-02 5 11111.1 222 33* 0.19 8 E 0.218 E-02 11111 22 333 0.218E CL237E-02 6 0.23 7E CL25 7E-02 111111 222 333* 0.25 7E 0.2 7 7E-02 7 111111 222 333 CL277E CL297E-02 111111 222 33
- 0.297E CL317E-02 8 11111 222 33 0.317 E-0 2 - 0.33 6 E-02 11111 222 333
- CL336E CL356E-02 9 11111 222 33 11111 222 33
- 0.3 5 6 E 0.3 76E-02 11111 222 33 0.37 6 E CL396E-02 0 11111 222 333
- 111111 222 33 4 11111 222 33 4*
11111 2222 33 444 111111 22222 333 444 11111 2222 333 44* 1111111 2222 33 444 111111 2222 333 444* 1111111 2222 333 4444 1111111 2222 333 444. 111111 2222 333 444 1111111 222 333 444 . 1111111 2222 333 444 . 111111 2222 333 44
- 111111 2222 3333 444 111111 2222 3333 444
- 1111111 2222 3333 444 11111111 22222 33333 4444 =
1111111 22222 3333 4444 11111111 22222 3333 4444 5* 1111111 2222 3333 44444 555 111111 22222 3333 444444 SD* 11111111 2222 33333 444444 5555 11111111 P2222 33333 444444 5555* 1111111 22222 3333 444444 5G5555 1111111 2222 3333 4444444 5555555* j 1111111 22222 33333 444444 5555555555555 l 1111111 2222 3333 44444 555555555555555555* i 111111 2222 3333 4444 5555555555555555555555555555 ( 11111 22222 333 4444 555555 * ! 1111 222 3333 444 5555 666666666666666666666 111 222 333 4444 555 6666 77777777777 6666666666 e 111 22 33 44 5555 66 777 888888 77777 6666666666 11 22 333 44 55 636 77 888 99999 8888 77777 666666666 .... ...................................9........7........6.......... l Figure 3.5.16 Salt concentration contours for Case Sa. 3-90 RAINFALL RECHARGE R AINF ALL RECH ARGE 2MO m - O O O O O O O . RIVER 1000 m 40 /;FA .[,,d* .L CSR . . ' , *, .,,'
- N
' ^ ' # YSTEM, s .; 's~ SALT SYSTEM O y . Om F (...i.v.o , 's ' ~ - TR AN SITION .... ZONE '.',.,.,, v g# , . - LIMESTONE SYSTEM .J 1000 m 4 HELD POTENTIAL # e, N ' (Dirichtet) Oeg " ' IM EisT d N E-S H AL E- boundary " . SANDSTONE' SYST Soo m NO plow 80uno g 2000 m - i l I l i O km 100 km 200 km 300 km 400 km DISTANCE (km) y 150 km RIVERS ' ' ' = 0 ! O 100 km 200 km 300 km 400 km ASSUMPTIONS - Steady-state conditions - Anisotropic hydraulic conductivity - Constant temperature and density - Flow can be described by Darcy's law - System is confined and fully saturated - Recharge to ground water is uniformly distributed over surf ace area - Discharge to four rivers - Boundary conditions according to figure on top of this sheet i Figure 3.6.1 Diagram of Case 6 including problem assumptions. 3-91 je 176km : l:
- 226Cm rl -
- /::.. .;;;;.. . ::..*. : : . :: :: } .:* ::*.. :{:. n
.....f..... ..* a.. . . . . * . . * . . . . : . . , ll
- j[ ..
/**'. * . . ARE AL E , XTENT . *:. .OF: .L AYER FLUVIAL SYSTEM . , , :.1 *.*.
- , . :. .:.. ,E*;
N ' f,* ,[, ./, *. o ./
- a,.* . [l, 2.0 m/dey Kv/Kh.0.10 $
..:,...~,.::.........:......::./..:.... . '. ..; *, :....,;.a::/*.*.*..***.
- /..,*. . . . . ., .. ,. : //.*. . *. ; *.*:,
,... ., , . . :u . . . . . : , , . . * : . ; .. W100 km rl: 300 km r! $kllb$./. % . 7: <.y$.n%. il$A$W!%$ <*49. .m,v. m Qrl Yt %sl+t. y 9,,y. ~.e,.;q .e$ >;;m*L } r s.j. L 3, 3 tj4,,4.. c. ,4;f .- isi :- . .. . . p ., .. ..,,4, s ,%q1/..gz< -y ,s, jk i /YN/',N ..y E .e.I' s. $*. q% IIN'p . . g.NMx y, .N, ,, s ,a.MIb7 h,m! ff , gi ';,; ;j, .f4Y AREAL EX ENT OF LAYER 2 WQ{51 E N L. 4 5 f t. E d , N# < .c. ,kg;l- LACUS RINE SYSTEM = 0.10 FOR ENTIRE SYSTEM ff. Yg* 4; e O *4 d,*j W*$ 1 i ~Y $,r$'.t.' K v0.01/K m/6ay h d*K*W;.%,y /< C"? * = ;.$ of.f: a% .
- ?
M. WF g937 s .9< %., . "y ,c ..,;;.i;% M.FY-3.'1-(! '?"M M' +[; ;4. s-*,yvy wr.x 4i&p .vi.f'.3 [.G.y;j,@. ip ;e g @*$-% :.r . wihlf
- eYl. .*:o W 14 5 k m ti.OM - ,,
7hMN bl&&S$*$%NW,N,,qEhd&.? i ,'.Y$ik 30 hm h 130 km M 240 km r ], W '% o g *. , ARE AL EXTENT OF LAYER 3 gO Oz y ) SALT SYSTEM )k gSg h3
- 1.0E-10 m/ day Kv/Kh . 0.10 gz1 ~
ARE AL EXTENT OF L AYER 4 * [k* N (> i -~* LIMESTONE SYSTEM $* 1 .a E l <z, 0.08 m/ day Kv/Kh. 0 to 1 gyE ;L; ;fj i ARE AL EXTENT OF L AYER 5 <*e ev g; M 4 LIMESTONE-SHkLE-SANOSTONE SYSTEu 0.S m/ day K v / K h = 0.10 i O, ' t :"*'W z y yn *, , .;b'JW.P%V f ~ ' h- 9 $ NYNNi ' i#5 l -= 100 km :l 240km 2 di ) n REAEk1TUUEO[ g
- . ALLUVIAL /F AN DELT A SYSTEM N d 0.00 m/ day Kv/Kh.0.10 0
- iG
_o_ Figure 3.6.2 Plan views Of Case 6 including hydraulic parameters. 3-92 Z 2000 - STARTING POINTS FOR TRAJECTORIES 1-4 AND 5-8 RESP. **=ne _ f l 1000 - *sa s' ,=,,,,"===.., , "C , 7 ((h r%v T' .""an.. Wum *= , ': :
- .. . Y,'
SEA LEVEL 0 - - ~ \~' beijiE:::: Q'!;:_ ,,
- 5h :im:::
18._
- iW.:::.&
, l.l:4 L x ;;.g... 4M -1000 - - ! [ , EVALUATION LINES 1-4 AND 5-8 RESP. ~ ^5 -- -2000 - FLUX CALCULATION T W' . r <1 l l l l
- x 0 100 km 200 km 300 km 400 km REOUIRED OUTPUT
- HEAD DISTRIBUTION ALONG EIGHT LINES IN TWO VERTICAL CUTS (Y = 15km and52.5km) ACCORDING TO FIGURE ABOVE - TRACKING OF ElGHT PATHLINES WITH STARTING POSITIONS IN THE TWO VERTICAL CUTS (Y : 15km and 52.5km) . - DISTRIBUTION OF FLUX IN HORIZONTAL RECTANGLE AT A LEVEL OF z=200. POSITION OF RECTANGLE IS MARKED IN FIGURE ABOVE Figure 3.6.3 Specifications for model Output for Case 6. 3-93 SWIFT II To facilitate using SWIFT II to simulate the complex geometry of this problem, a preprocessing program was written to adjust grid block sizes and depths. The resulting grid contained 40 blocks in the x direction and 19 blocks in the y direction. A vertical section of the grid is shown in Figure 3.6.4. Two approaches were used to simulate the rivers. For the east-west trending rivers, which are in layer 1, Dirichlet boundary conditions were used. Actually, the middle river cuts into layer 2 but, according to the problem definition, is to be treated as if it were in layer 1. Then, for the north-south trending river, which resides in layer 2, a well submodel was employed. This option permits the user to define a pressure that is maintained at the node center by allowing flow in or out of the "well" according to the hydraulic gradient near the well. An areal recharge option is available in the SWIFT II code and was used for this problem. However, the option assumes all recharge goes into the top layer (layer 1). This caused some difficulties because layer 1 does not exist over all of the top of this model (see Fig. 3.6.1). In fact, even layer 2 does not cover the whole region. To overcome this difficulty and still use the recharge option, very thin blocks (.0001 meter) were defined for the region where layers 1 and 2 do not exist (see Fig. 3.6.4). Hydraulic parameters were assigned to these blocks that were equal to those of the underlying layers. Results from the SWIFT 11 code are shown in Figures 3.6.5-3.6.13. General flow directions are indicated by the contours shown in Hgures 3.6.5-3.6.9. Hydraulic-head profiles required for comparison with other HYDROCOIN results are given in Figures 3.6.10 and 3.6.11. Finally, predicted pathlines, also required for HYDROCOIN comparison, are shown in Figures 3.6.12 and 3.6.13. Vertical fluxes across a plane located in the salt layer were not calculated because of the necessity of creating a post processor. Table 3.6.1 provides a listing of the accumulated travel times and distances for the individual flow paths simulated using the particle tracking algorithm. In general, the results are consistent with the problem definition. However, there is one major exception. In layers 1-3 just east of the north-south corner, model predicted hydraulic heads are too high to be realistic. This is shown both by the contoers in Figures 3.6.5-3.6.7 and the hydraulic-head pro-files in Figures 3.6.10 and 3.6.11. In additiori, these hydraulic heads have a severe effect on the lower two model predicted pathlines shown on Figures 3.6.12 ' and 3.6.13. These very large hydraulic heads are located in the two model j coiumr.s that represent the salt system. This is significant because of the extreme contrast in hydraulic conductivities between the salt system and the surrounding rocks. This contrast is about nine orders of magnitude. An additional complication was that a relatively large amount of recharge was being forced into a unit of low conductivity in this region. The problem of high hydraulic heads resulted from the inability of the SWIFT II code to arrive at an accurate numerical solution in this region. Numerous attempts were made to improve the solution. These included forcing the code to take more iterations to solve the problem, changing the direction 3-94 a A ' \ l I f I 1 I I I I \ l I I I \ l I I I \ l 1 I I \ l ! , I I \ l _I ~ ~ l I \ l I I i 1 I I T I I \ l .1 5 I J \ I 1 ? I I \ l Y I I i Y Y X I I I X X X 3 I I I X X X u X X .X o I I I i I I X X X , I I I X X X w I I I / X Y $ I I I X X X v I I I X X X T y * \l I X X 1 1 X / 7 e \ I y X_ Y J h I I \ Y a 17 I I I , I/ f I y 5-i ! T/ I i I 5 l I/ 1 1 y C / I I / I I Y I I I I 1 l 1 Y Y t I X Y II I YlIy 3-95 HYDROCOIN LEVEL - CASE 6 LAYER 1 '8 ' I r i f I fl i 120.0 - - l = , I uJ 90.0 - - o o % 8 , W F . .I o n r t s W a { S 1 0 $ $ m o ~ u a 'no 'O ; 60.0 - - , 2 g g o - o 2 $ 8 $ o lk o o 1 l 30.0 - - 8 l - S o.o ! I I I I I o.0 40.0 ao.o 120.0 160.0 200.0 240.0 280.0 320.0 360.0 400.0 KILOMETERS l l Figure 3.6.S Hydraulic-head contours for layer 1 of Case 6-SWIFT II. I l l i O l l l l I O v l p - g ~. s-O $ 2 5 00,t ? / - o o u o 1400 o e 0091 * - g J
- E 5-WE W &
-1600-~ o H >W 6 W (m W o " 2 < O 2 J J c 2000 - o O M v O __ 00C3_ ._: y .e ~ h ==rr== 0 0 2 L
- e O
e > 0031 ,o I - O 5 00Ci d O I Q o d 1400 - g g e o t e o 'w" 1500 6 e i f f i o o O O o O o O o o o o e o W N O C M S H 3.L 3 W O"II M 3-97 r C l l l l 8 O _ . O 3 800 Y__ ~ m m 4 $ \ -- _ - / h h \ 006 E o O $ _ , O u 1000 * % W 0011 L (O ~ O N 1200 - N 0 W g J n C O w W * >g 0001- o F g w W> f o 2 W L J $ < 1400 N O a c 2 J J o O O _W . -' C! y u O _ - _ O e m 1400 e j f 0036 {; z ; s O m - 00C L d MX N. e z N O @ 1400 6 M e I o .n O W - 1500 - g e i i I t O O O O O O O b 0 D N O O n w w S W 313 W O 71)I 3-98 l l O. O I I O I I y 009 l O. 1 O I l M i
- 009 O ,
l 6 l - - a y to b 002 3, O. e i @ - - O e , l W 009 m ; . m i N " 4 O o 006 - d (/)
- e wg b i at NW @
t W - 000 L - OF 2 i y $oc - - dW S J4 OE + Za oogg O o 6 O. - J t o O - ~ Og u O D E L --
- j l @ , l i
O fu i > .t o , OOC L - d Z j N 2 1 e v O 2 00# L - g q 1 m o m O E - OO9L - g g 9 w ? I O O O O O O* d d d d d W N O @ CO e- e i SW313WOllM 3-99 l l 0 - 0 0 4 ooW 0 i ' 0 6 3 oo@ 0 - i ' 0 I 2 I ooN 3 T F I 0 W S 6 i ' - 0 6 E oom 8 e S 2 s a A C C oo@ 0. S f - , ' o 1 0 R 5 4 L5 2E r e E R ooo" l T y a ._ V E E Y i ' 0. E l L A 0M 0 f r o NL oo,"- 2O s r I O L u I o C i ' 0. 0 K t n O ooN" 6 o c R 1 d D a Y 0 e h H Ooo" i ' c - 0 i c 2 l u 1 a - ;c 0 y - oov" H 0 8 9 _ 6 3 oom" 0 i ' e 0 r u 4 i g _ F _ i 0 0 0 0 0 0 00 0 0 0 0 0 0 5 2 9 6 3 1 1 )Em 0 l g2OJ_M Tg 1i i il1i j! ! ll ;j1il ,! j lllllll ;, 0 o - 0 )* 3 A - 0 - 4 5 s - - b#bp - 5 \ I I \ .j r\ T F I W *. , S .a # _ - 6 _ ' Q) 6 e E _ ~ 0 g s a S m A 0 \ " \ 0 4 0 C r 0 o C 0 0 _ - 0 1 f _- _ 4 3 1 5 _
- 1 L1 1 E
V 1 T= _ ( g. l i s e EI F Y ,- f \ _ , \ o r LW, _
- NS-1 4 - *
* [ g d p a I e O o h C O S E N Y,,.-{- o 0, i l c u O9 \ x R I # ., a r D L O d b y Y H o s% . :. g4 H 0 / o e :m m e 1
- o6 4 -
*/ -- u O o / .j g "@ o i o /
- I o e e m : n a w
/ / E = a w=m / // / w e >sn / W ' I ./ / - o $. z*m mi f ./ z e a o m - 4 = a m I :l oy :
- o V) ?
o a: w z - I 1 :/ _ o. OO - Y 5 o a ' :/ o* e I f', } y ~~ j{/ = , e M e a o C o I I I I o o o e o 6 o o o o o - o e o m N w p (sJelew) OV3H OI10V80AH 3-102 f /1 1 \ l . ,I fI \ l I 7I I I \ I I II I I \ l I /1 I I 1 I I // I I \ l I /Y I I \ l I II I I \ l I // I I \ I I /1 J ~ LI \ l I /t I I \ l I // C I I , 1 I I // y I I I X X / // ? I I I X X / JX
- I I I X X / /X @
I I I X X / //
- A I I I X X / ff I I I X X / // '
5 I I I X X / ff I I I X X / /X T ~ I I I X X//X \I I X X / /X C 't I X X/ /7 .5 \ l X YJ /Y 7 kI ! 7 27 g TI y T/ I K I\ // N F /i/ 17 d T / Y 1/ ei Y l I Y e f I I f 5 I I I .T Y / Y II I i i Y' Y II Y 3-103 Z 2000 - STARTING POINTS FOR TRAJECTORIES 1-4 AND 5-8 RESP. A s '*4m-- _ [ __ 1000 - hs ; mm...,, "' & - r- _ v * ~- ~ ~-- m,,,, y.y' %'i:(:i:. SEA LEVEL % W ..:s... 0 - ~ *~ ^ - % .;;: ::5 , .;. .- ( N 4,% < 4 4 , ; -1000 - ~- _< i - , i)), EVALUATION LINES 1-4 AND 5-8 RESP. A - FLUX CALCULATION ' M, <l I i l l L X 0 100 km 200 km 300 km 400 km REOUIRED OUTPUT - HEAD DISTRIBUTION ALONG EIGHT LINES IN TWO VERTICAL CUTS ACCORDING TO FIGURE ABOVE - TRACKING OF E!GHT PATHLINES WITH S TARTING POSITIONS IN THE TWO VERTICAL CUTS - DISTRIBUTION OF FLUX IN HORIZONTAL RECTANGLE AT A LEVEL OF z= 200. POSITION OF RECTANGLE IS MARKED IN FIGURE ABOVE Figure 3.6.13 Pathlines 5-8 for Case 6-SWIFT II. 3-104 Table 3.6.1 Accumulated travel times and distances for individual paths predicted for Case 6-SWIFT II. Accumulated Length Accumulated Time Path (m) (s) (yr) 1 45,013 8.01x1010 [ 2,540] 2 47,587 8.82x1010 [ 2,797] [ 3 290,094 6.729x1011 [21,338] 4 290,039 7.16x1011 [22,704] 5 45,016 8.01x1010 [ 2,540] O 47,591 8.821x1010 [ 2,797] 7 290,094 6.729x1011 [21,338] 8 290,039 7.16x1011 [22,704] of the sweep for the two-line successive-overrelaxation solution (L250R) tech-nique used to solve this problem, and finally reducing the grid spacing in the
- x direction by a factor of two. Forcing additional iterations had no effect on the accuracy of the solution. Changing the sweep direction for L250R, in ;
this case from the x direction to the y diration, significantly reduced the ' number of iterations taken before the convergence criterion was met but did j not affect the model results. Finally, reducing the grid spacing did bring , the heeds down by about 100 meters but did not completely solve the probles ! Undoubtedly, the Jroblem could be resolved by further reducing the grid spacing, especially by doing so in both the x and z directions. However, further grid * { reduction was deemed economically unwarranted based on the additional insights
- that would be provided. In any event, these large hydraulic heads appear to i be confined to an area very close to the two columas mentioned and only seem
! to affect layers 1-3. j USGS l Of the three cases addressed here, Case 6 most closely approximates conditions that the authors of the USGS rodel apparently intended the code to simulate. , As a result, representation of Case 6 required few of the parameter adl'Jstments employed in Cases 1 and 2. ' The problem is represented in a 15 row x 41 column x 5 layer grid. The grid l spacing in both the x and y directions is set at a constant value of 10,000 [ meters. The block centers in column 41, in order to correctly represent the ' held potential boundaries, are located at the eastern edge of the simulated area. The western boundary of the problem, because of the constant block size employed, corresponds tn the block cen**rs of column 1. Half of each block in columns 1 and 41, therefore, lies outside the boundaries of the problem. I i 3-105 l .- ,. . -. - - . _ - _ _ _ _ _ ~ _ ... . _ _ _ - . Each layer in the grid is defined to represent one of the geologic strata present at the east edge of the problem. The transition zone in the west half of the model is represented by layers 4 and 5, and the boundaries of layers 3, 4, and 5 are projected from the west side of the transition zone to the west edge of the model. Layers 1 and 2 correspond directly to geologic units and are inactive in the west side of the problem where the geologic units do not i occur. The transmissivity (in most blocks) is based on the hydraulic conductivity of the layer and the vertical thickness of the layer at the block center. The leakance (in most blocks) is based on the harmonic mean vertical hydraulic con-ductivity and vertical spacing between the block and the underlying block. Blocks in columns 1 and 41 provided exceptiens to this description. Because the block-center thickness in layers 3, 4, and 5 in column 1 and layer 1 in column 41 was zero, the transmissivities are based on the average layer thick-ness at the block boundaries. Leakance values are corrected to represent vertical flow through only that half of each block within the area of the problem. All rivers defined in the problem are simulated as held potential boundaries, as are the fixed potentials in layers 4 and 5 at the east edge of the problem. , Areally distributed recharge was represented through option 3 of the Recharge Module, which applies the specified rate of recharge to the uppermost active l block in each vertical grid line. The rate of recharge was reduced in columns , 1 and 41 tn account for the area of the block that lay outside the boundary of the problem. L The discretization used nere may provide an adequate value for the average hydraulic head in each block and for the average ground-water flux across each block face; however, its representation of heterogeneity is extremely crude. The calculated potential field cannot be expected to represent the abrupt gradi-1 ent deflections that should occur at unit boundaries, and particle tracks con- a structed from the results of this model cannot be expected to adequately ' represent the deflection of flow lines at unit boundaries. A more detailed representation of the heterogeneity could have been achieved by employing 13 layers with two thin layers representing each layer boundary and a , variable block size with smaller blocks near unit boundaries. This construc- : i tion, although well within the 80-layer limit provided by the code, would [ require more dynamic memory than can be accessed by Microsof t Fortran version
- 3.2 in an IBM Personal Computer AT.
j Model results are displayed in terms of hydraulic-head contours (Figs. 3.6.14-
- 3.6.18) and hydraulic-head profiles (Figs. 3.6.19 and 3.6.20). These results are very similar to the SWIFT II results, even to the extent of having unreal-istically high hydraulic heads east of the north-south river. However, the USGS code, using the strongly implicit procedure (SIP) as a matrix solver, does not l 5
produce as high heads as SWIFT II does. In addition, the USGS results only show the high heads in layers 1 and 2 and not as much in layer 3. Particle tracking as well as the required vertical fluxes were not calculated since calculation of these outputs would require the development of post processors, 3-106 . ~ - _ - a T 3 e 3 i a t g t a R S a R G S a R U 6 a i e s a C a t f o 1 a t OO,.. . r 4 i e ,i y a t a l a OOe- 0 s R s r u 5 R o t n
- OOe- ? o c
5 R d a s R e h a R 5 R i c l 5 R u a 5 R r d 5 e y H 3 i 4 e 4 1 i 5 6 5 t 3 a t e r 3 i u g a e i F 5 a a f a R - - - ~ - * - _ - - .yoy ._ _-. .. .. ._. . =. . 1 s a i e a i a a a = i s i s e i i i s j s e i r i q [ s a a i s [i i a i v1800 - ~ - ,1700 i 1 ~ - 1200 1500 1100 1000 _ s300 ,,,, o - m i, ~ , - i l I - \ 7 - < I J _ 1200- - 1- l 1400 - . . . . . . . . . t , . .I. . J . . . .,. . . . < . . L , . . . . . . . . . I Figure 3.6.15 Hydraulic-head contours, layer 2 of Case 6-USGS. ' [; 5 i 5 . 5 . 5 . I 5 . S . G 5 S U 5 . 6 e 5 . s , a C 5 . f o 5 . 3 5 . r 5 ooee . e y P l a 5 j ooNe . , s 5 . r u 5 t o _ [ t 5 .uOo . 4 . n o 5 . c 5 . d - abOo a - 5 - . e OO w OOMe h - OONe c i 5 . l u 4
- moo
- a r
5 d y 5 . H *"goo . 6 . 5 . 1 5 . 6 - 8 . 3 . **oO e r 6 . u - - g _ - i , 5 . _ F - 5 . " *OO - 5 . 5 . w* Hoe I\ i i , 1 I lllil I I :lll i , i! I I ll. , L , ~ a M M E - a ~ - . _ _ _ - , a . , a O.. , a , i a e , e . , o.O S a G S u a U a a 6 o.N e s s e a C .
- a f _
o.oe o a a 4 - u e r e . y O.# a - l s a . a a s r _ 5 a u aO.O o a e t - n _. u a o c _ a e d . F aa.O a a . e 5 a h 3 a i c - l 3 aM.O a u _ a 3 a r d 8 a y . H 3 e *U o 7 =- 1 8 a . - 6 . - 3 e . 3 - 3 e - e
- yO r
._ u ' a g i ' a F ' a -- 8 e _ - " = - - _ _ _ - ~ =. n Wepwo ,t I I 1 l.! ll'l l l , , , . . . .ii. . . 009 . d w . m - o . e e 00L - + U M - o 008 - m - u W ~ 006 . - . u ~ 000L . a 8 C o . v . e 00LL - 2 .c t .u- . 00EL - E u . u>, - = 00CL m 9 e = . G e 00tl - s .ee . W 009L - j { l j j jI j I I A A E 3-111
- HYDROCOIN LEVEL 1-CASE 6 USGS-3D MODEL 1 HYDRAULIC HEAD ALONG SPECIFIED LINES 1-4 l l l l l (0
w , 8 ' LINE 1 e2000 - ---- LINE 2 - .' E ..... LINE 3 v -- LINE 4 O 4 m 1500 -
- Y I , A - ~~~
N g / % ; ~ . ~3 1000 - %, .**** ~ ~ ~. N,, 4 % E N O . ~* * * . . . *: I 500 - 100,000 150,000 200,000 250,000 300,000 350,000 400,000 4 DISTANCE IN X-DIRECTION (meters) Figure 3.6.19 Hydraulic-head profiles 1-4 for Case 6-USGS. l HYDROCOIN LEVEL 1 CASE 6 USGS-3D MODEL HYDR AULIC HEAD ALONG SPECIFIED LINES 5-8 ^ l i I I i V) s W LINE 5 ---- LIN E 6 8 - .. . .. LIN E 7 - [2000 - - LINE 8 O <C w LU 1500 - - $ Z ,^~~~~ o NJ/ ' ' - ~ ,~. a . W :..........**** * ~ ~ , ' - ,~ 3 1000 - 4 N ...,***.~~.
- OC -
l O ' ~~ :. . I 500 - i i ; i 100,000 200,000 300,000 400,000 t Y= 52SOO m DISTANCE IN X-DIRECTION (meters) Figure 3.6.20 Hydraulic-head profiles 5-8 for Case 6-USGS. 3.7 Case 7 Problem Description A vertical plane of a low-level radioactive waste site situated in a saturated layered argillaceous medium is to be simulated in this case. The argillaceous medium ranges between 17 and 22 meters thick and overlies a 3-meter-thick, high hydraulic conductivity aquifer. Both the hydraulic conductivity and porosity of the argillaceous media vary spatially in the vertical direction but are con-stant in the horizontal direction. The 3-meter-thick aquifer is both homogene-ous and isotropic. The case is simulated for steady-state conditions. A sche-matic of the problem is presented in Figure 3.7.1. There are two variants for this case. For the first variant, the hydraulic conductivity of the two con-crete bunkers and their associated concrete caps is 10 6 m/s, which is larger than the hydraulic conductivity of the surroundin argillaceous media. For the secondvariant,thehydraulicconductivityis10gO m/s, which is smaller than the hydraulic conductivities of the surrounding argillaceous media. In the vertical direction, the hydraulic conductivity of the aquifer and argillaceous media is described by: log K(Z) = [ log K(2 2 ) - log K(Z )][(Z - Z )/(Z 2 3 3 - Z i )] + log K(Z i ) where K(Z) = hydraulic conductivity at height, Z, 2 and 22= sequentially specified heights above the datum, and 3 Z = height of interest. Values of K and Z are specified below: Z(m) K(m/s) 25 10 6 20 10 7 15 10 8 10 10 9 7.5 10 8 5 10 7 3 10 6 i i 0 10 6 Similarly, porosity is described by: { E = [c(22 ) - c(Z )][Z - Z ]/[Z 2 i 3 2 ] + c(Z ) 3 3 l l f i 3-114 i! ll IlllIl, i11 l ll
- t. .
2i k i o* ._ N e _' / r e h / , V . 0 = / . m r 0 = M / _ 7 U , T D / A , c ~ ~ /RY 7 A e m s = S / DNU C a z - ,. O r m h 1 m l 0 m 1 r /B f o r i e l _ - 0 e , 1 2 / d o m - = l E m o, ' g 5 m e / F t a u _ 4'm + 1 p = e y a' f q / o m c n 1 Wo o t C 5 1 /O L F 1 - 2 . 1 . - c,1 [ = / 7 3 e , i TT / r u - m g 1 i F / , 0 0 I m 5 0 = / ON , I 1 2 M / 2 ,- 7 U . T A D V / _ 0I 0I y E L A 5 m 2 / E C = L S C A S L A ^ h /D _' sl c L T ., A N , C O p_ , Z S Rd I i i T I R R Ene , E O - V H wd& - c l 1)' . !. : . P where c(Z) = porosity at height Z. 4- Z, Z , and Z2 were defined above. i Values of c and Z as specified by: Z(m) c . 25 0.3 '
- 20 0.25
] 5 0.2 3 0.1 0 0.1 Simulations Performed This problem was simulated with two finite-element codes, UNSAT2 and FEMWATER. UNSAT2 The head calculations for both variants were performed with the UNSAT2 code. Head values, output from the UNSAT2 code, were then used to calculate vertical velocities along three horizontal lines running along the top and bottom and through the center of the concrete bunkers. (NOTE: .The UNSAT2 code does not internally calculate velocities.) One of the parameters output by UNSAT2 is ; discharge from or recharge to the modeled system at each boundary node. Using ; the discharge data and the flow net theory presented in Reference 8, the UNSAT2 t code was used to generate streamlines for the modeled domain. The streamline - , cutput from UNSAT2 was input into a Sandia contouring program that generated ' both streamline contours and travel times along the streamlines. The finite-element grid for this case is presented in Figure 3.7.2. The same grid is used for both variants of this case. Grid spacing in the vertical di-rection is 1 meter. However, appropriate adjustments in grid spacing were made , along the upper boundary such that the grid corresponds with the upper boundary. In the horizontal direction, grid spacing varies between 2.5 meters and 7.9 meters. In the horizontal direction, the ratio of one block length to next 3 block length is kept at approximately 1.2 or less, except at the middle block of the grid where the factor is approximately 1.35. The 2.5-meter grid lengths were used at the left and right boundaries and at the four vertical sides of the concrete bunkers where they meet the argillaceous media. The 7.9-meter block lengths are located one-half the distance between the left boundary and the left side of the left bunker, and one-half the distance between the right boundary and the right side of the right bunker. Grid elements consisted of rectangular blocks, except along the upper boundary where elements consisted of 4 either trapezoids or triangles. Values of hydraulic conductivity and porosity ] are assigned to each element based on the elevation of the midpoint of the
- element. Because of the irregular shape of the elements along the upper bound-j ary, 70 sets of hydraulic conductivities and porosities are assigned to the i various elements.
1 i ! 3-116 2 T A G N U 7 e s a C r o , f d i r . g i l e d _. o M 2 7 3 0 2 e r u =0 g i
- 1 F 1 =
L1 A: TL NA OC ZIT I RR OE HV E L A C S Tg ~ i Variant 1 Results !
- Head contours for variant 1 are presented in Figure 3.7.3. Although head contours are not required for HYDROCOIN output, they are presented here to '
j point out some important results. First, the contours of the variant 1 case i at the lef tmost part of the region are almost vertical. They appear to be ; strongly influenced by the constant-head boundary condition at the left bound-ary. This would seem to indicate a near horizontal flow into the region from ; the left. Second, the head contours in the aquifer are very nearly vertical, t ! indicating a horizontal flow path in the aquifer. Third, the heads along the
- right boundary, where the hydraulic gradient has a rather large value of 1, appear to have a strong influence on the flow field below the left side of the lef t bunker to the right boundary. Fourth, an unsaturated zone forms at the r lover right-hand corner of the modeled region just above the aquifer. Because I the hydraulic conductivity was kept constant for all values of moisture content i ard the moisture content was kept equal to the porosity for all values of nega- !
t've pressure head, the unsaturated zone was essentially modeled as a saturated zane. Streamlines for variant 1 are presented in Figure 3.7.4. Streamlines entering . the flow system from the lef t boundary are nearly horizontal near the boundary i but quickly drop downward and into the aquifer. Streamlines entering along the top boundary predominantly drop downward to the aquifer. All streamlines enter-
- ing the system, except streamlines 1 and 2, leave the system through the aquifer.
j Streamline 1 is located in a region where the gradient is upward. It immediately ; leaves the flow system. Streamline 2 enters the system along the left boundary ; and leaves along the top boundary (a particle released along the top boundary l' at the point where streamline 2 intersects the top boundary will also immediately leave the flow system). Along the streamline and the boundary above it, some , numerical error appears to be occurring. UNSAT2 indicates a downward flow gradient but actually calculates flow leaving the system, which would require an upward gradient. , Vertical Darcy velocities through three horizontal lines for variant 1 are shown on Figures 3.7.5, 3.7.6, and 3.7.7. For variant 1, vertical Darcy velocities on i all three figures are negative except at the very right boundary on Figures 3.7.6 and 3.7.7. This indicates a predominantly downward flow. Vertical Darcy velo- ,
- cities through the blocks, are greater than those for the surrounding media.
- Further, vertical Darcy velocities on Figures 3.7.5 and 3.7.6 are almost equiva-W lent except at the right boundary. However, the vertical Darcy velocities on Figure 3.7.7 are lower than those on Figures 3.7.5 and ?.7.6 by approximately ,
, 35 percent. This is attributed to hydraulic conductivity changes in the modeled ' ! system. Along the horizontal line at the base of the bunkers, the hydraulic i conductivity is 10 9 m/s. Both above and below this line, the hydraulic conduc-i tivity is higher than at the line by a factor of 30 percent to 60 percent ] because of the interpolation method defined for this case. Discretizations that
- do not have the depth of the base of the bunkers occurring at the midpoint of a finite element will not have a 10 9 m/s hydraulic conductivity and, hence, head calculations do not account for it. This is the cause of the differences in the vertical Darcy velocity profiles. In addition, vertical velocities at the j edge of the blocks appear peaked. This is attributed to using the hydraulic conductivity of the bunkers in the vertical Darcy velocity calculations rather
! than some type of average hydraulic conductivity between the bunkers and the i
- 3-118 l . . _ _ _ __ . _ _ . .
D E T A R UN T AI O SG NE UR . . 2 . T . A S .,,,,, N ,,, U 5 - 0 iOH s/i ) f l 7 ~ 2 1_ l /l2 (e$?s" e f)/ (_u_ . ._ s _ l a / s - _ . _ . _ _ C Y - z_ __ f f>/nYpO _ o - / (o - 1 P / /0 / -M t t i n a / 2 _ _ .. _ r /ICGff . _ / a / //d;g/f _ _ v /' %WFf f r o s p I- _ _ - _. ._ f /l _ _ .'fh ^ _. . , _ s v : s r yJ u o ( - _ t n ]/ . }j m t, . - . o c / u j " 1 . -'i _ _ _ d . !I a f_1l n4, \ . e , . 3 , F _. _ . . _. h 7_ - ._ _ c 1 f jl/ 1 i l _ _ _ u ~ ..g[, . . a ~ - r T . d ' / l i / )f (. H y - Y - i . _. 3 /.%~[ / _ / - _. ._ 7 -- / L[I _. 3 ' / f . e r ' l / l _. i u g l l / F 1 f f l r _ I. / , lf 1 I Y Cw / h_ d a t / .2 k = t .5 * / , M' b a - i e- . - 2 3-12G O O i I I N ! = . N N E E w O e, y < z 1 i Q ?
- o
/ ,
- 11
/ 11 m $ x / o n O U U O i ~ o I y - m a i g w I 7 m I 9 \ r.-- n ~ I . I - E ! i t l_- W $ { OQ s o _g O g u r- 4 z I F 5 m i - O s= 'y, j \1 2 i l .t i f ~ / 8 -f* l I % I .x t 3 ~. M e ! O 5 c o m o I e I (s/w 3 0g.x) A1lcoagA 3yol1ggA 3-121 o i i o N N a E /s e g 7 /E & o / e o 5 1 - / y- 4 / 11 a
- o 8 u l o
- m a 2 ; - . a / 3 o M g ' I (~ s- ,a - N il w \ '~~- of g o I o a g - o4 -p o , S l ! 3 8 l s~3 " f t t / 'G - I _ o { ' / m 1 0 l I / E / a 9 m e I a o c o e O I (s/w 6 0 8 X) ^110013A lyol183A 3-122 tzt-t I VERTICAL VELOCITY (x10 8 m/s) l d I , o m o E o ' .~ N i a ! c. 9 - I
- I 1 m -
g o / I , ?. z I w ; ~ N ao s- i ,m , aH :- e>a zo - I li ~ 0o 's ,E m 1 ;3 ^v i il T Y I O x I n 3 I o " n n - - un - ' l ' n o ol / of O d y O a 3 / z s 3 1 > m / s M I I o surrounding porous media. The hydraulic conductivity of the bunkers is higher than that of the surrounding media by two to three orders of magnitude, thereby accounting for the peaks in the vertical Darcy velocities. Variant 2 Results Head contours for variant 2 are presented on Figure 3.7.8. For the most part, they show similar shapes and effects as those presented for variant 1. The only difference that appears to occur is that the contours turn more upward in the vicinity of the bunkers than for variant 1. Streamlines for variant 2 are presented in Figure 3.7.9. The streamlines for variant 2 show similar effects as those for variant 1. The main differences are that the contours defurm slightly differently in the vicinity of the bunkers and caps than those in variant 1 and that streamline 2 flows into the system and out the right boundary rather than out the lef t boundary as in variant 1. The vertical Darcy velocities are presented in Figures 3.7.5, 3.7.6, and 3.7.7. Outside the concrete bunkers, the shapes of the Darcy velocity profiles are similar to those for variant 1, but the magnitude of vertical Darcy velocities are slightly greater for variant 1 than for variant 2. Inside the bunkers the velocities are approximately two orders of magnitude less for variant 2 than for variant 1. This is attributable to the lower hydraulic conductivity of the bunkers for variant 2 than for variant 1. Because the hydraulic conductivity of the bunkers is near that of surrounding porous media, the peaked vertical Darcy velocities do not occur at the edge of the bunkers as they do in variant 1. FEHWATER Hydraulic heads and Darcy velocities for this problem were also predicted with the FERWATER code. In general, these results are very similar to the UNSAT2 results (see Figs. 3.7.10-3.7.18). However, two discretizations were used with the FERdVER code. The reason for two different discretizations is given below along wit.. other details of the FEMWATER solutions. An understanding of this simulation can be gained by examining the vertical Darcy velocity profiles (Figs. 3.7.14, 3.7.15, and 3.7.16). This discussion will focus on the velocity trend over the entire modeled region length, includ-ing the velocity profiles across the trenches. The Darcy velocities, along all three lines for both design conditions, tend to increase downward when proceed-ing from left to right except for the area within the trenches and at the extreme right boundary. These results are consistent with the imposed boundary conditions that converge flowlines in the lower right region of the site. Near the right boundary, the Darcy velocity decreases slightly at the 10- and 12.5-meter lines while increasing sharply at the 15-meter line. This difference is due to the right-hand boundary condition and is consistent with the equi-potential plots for the two design conditions (see Figs. 3.7.12 and 3.7.13). , Across the area of the trenches, differences in Darcy velocities exist for the two design conditions. Also, there are substantial differences in the Darcy velocities between the three lines for the high conductivity design condition. For the low conductivity design condition, there are no appreciable differences in the shape and magnitude of the velocity profiles for the three lines con-sidered. These results display tt orinciple of refraction of ground-water 3-124 N }'._..___..-,..,,"s ____. _ _ . J *_ /l11' ' l / l f I I i - ..-._--. .--._. ..L....__.. .._.__,/ r / l / t}~([' N l HYDRAULIC ../ / HEAD / ./ / / l / /J ./ - - ho // s CONTOURS / 4- (& lifQ [7 / j/ , 43 _.4/04kl / / - / e:K UNSATURATED Y g / '_ / / s f' _ A p [,% 3 'n [./ ~-C- s _- / [G b frQ : :nd : 61.5h W, W>==r- gr.i; REGION ._y _/ g c., -..__. . _~ 3~5 _.-- S .. --_h l r ._ ._J [ s ,- f ._ _ [h2$%: hf$b _f 4 -- - .-. . . ..- . . .-_ d-{-//f(// _ _ - T 20 ::: E Z::: r E 15 Yit 10 _u ~T d i 5....._ _y. . -.... i 4. _ -. . _ .. .. -- .-. - _. .-. . .. -.. .--J L Figure 3.7.8 Hydraulic-head contours for variant 2 of Case 7-UNSAT2. ) _b N 3 A / m / a y N ea I / .2 - / L e ,. . e / O i u m e ^ l ed D [ t n m N. m I ? N a .? 4 m a *f : 2 3-126 25.0 % ~ ~ ~ ~ __ 20.0 \\\ ______'T'~__ / / ) y\t I / \\ I ~ E 15.0 til 6-- 1 w E O O m G M O o 10.0 O e N s.o - i I i I 1 1 ,,,, so.o 75.o Soo o 325.0 150 0 # X - COORDINATE (meters) Figure 3.7.10 Initial finite-element discretization for Case 7-FEMWATER. --y m - - -- - , ,,v-- -,-- -- ,-g .7 w- - -----,r -r- - -- -r- g 7--- --- -- - - 0 . 0 _ - 0 - 2 _ - i - l 7 o. 5 7 1 R E . T . A _ W M _
- o. E o F 5
1 7 e s . a . C r y I I 0 5 2 1 ) s r e i f r o n o t t e a m i z i j I \ ( t r E e j \ r T c 0A 1 _ _ s I 0 N i 0 I d 1 l D t T I \ \ l / l R O O n e m e T C l e . ) - _ 0 - e \ t I 5 X i _ _ 7 n \ i . f _ - \ _ i l . _ a - i n - F _ 0 _ 0 1 5 1 _ 7 3 _ e r g u 0 i . - I 2 5 F _ _ 0 0 0
- o. o.
5 o 5 0 5 2 2 1 1
- 3 .Ew D<bOEhoiN t b s
tO$ a 1 !, g '%%,',f *- ' *g '% }" p9(," =* "l N N N g\\t / f =..-- ----- 7- , g g s g \gs 3 ( o \ \ \ \\ \ 4 \ )ig\\ \ % N g \ s' ~~---- o \ t t \\ \ l gI g\ \ \ \ o - o \ o \ \\ \i Ilii Ig O s~~~_.--_ s , \ \ l It o \ g i l\ s s
- 5
) \\\\ s'~~~~_--, o 2 i 1 gi t n 9 - o l Il 1 N s'~ 2 0 1 i \ t' ; \, \ l .~~ - - - - . " \ s' \g\ w \ m - % g\ u , e -= % i I s,\ g%z- = g2 : - :oo 5 Og s a , \ \g E s 5 s, lgi - ----, ; ;- ss N ,i \ o 2: s -
- 15 5
\ o - s % * * *\.N N g* - e o I -- ---., O g 'I \\ o o o Il \
- U a
Y o i g o v / \ - g x y s ~ e I / - - - - - - , 9 2 o v \ N % \ \ Q 5-g o e ~ % $ d j ----w-- % %-- o , m M i %g% R 4 s & % o o g - W h % N g e. W o 5 l l I 1 I o o o o d d. 6 5 u . - (sJenow) 31VNIGBOOO - Z 3-129 i ) I l I i i i ,/ i j r / / ; I I / ) ! / / / , / / e 15.0 - m / / l 5 / / l,, ?, ~ ~~~ ~ ' / / l! l , / Jf %' fK- ~ ' ,4,) e a / / ,/ g s- m m' rl m-N = _. - "s% s
- y 1 10.0 O
/I / j Q _ ., * , / _~~~~,,_ a's 8 l 7 f ' I / / l / s' / s' '/ / f_ Q 5 ~ l I I I I I / / / I l I l i l l l l l l l l \ 1 I l l l l l I I l l l I I I I l l I l i l I I 22.0 I 2 0.C I 18.0 I 16.0 I 14.0 I 12.0 I 10. 8.0 I 6.0 I 4. 2.C f , i 24.0 I i I i i li l i l il i i d s i i 25.0 50.0- 75.0 100.0 125.0 150.0 175.0 X - COORDINATE (meters) Figure 3.7.13 Hydraulic-head contours for variant 2 of Case 7-FEMWATER. i ! . . _ _ _ _ . - - . _ ~ . . _ . - - . . . _ - __ _ _ _ _ _ _ _ _ _ . _ _ _ _ _ . _ _ . . . - - - . - - - - - - , - 40 i i i i i 1 i i I j e t O 7 ---- K = 10-6m/s . . . . . . K = 10- 1 om/s x 2.0 -
- a N
E "m : .,,,, - l* ^ . l*-l\ ~ y- _ ..-~ ~ ~ -. a : . ..... A O - O -2.0 - t t 'c' : si - a w \ ,' i .\ s s- , i , r'. N
- s t g .--
'sj i j y -4.0 - ' f ..~."'s' '**- b ,. [ YO \ 8 % i.- 5m %, s % : f- <t s s Q -6.0 - \gr - _ a N <C O - -8.0 - - b-E W > I I I I I I I I I 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 180.0 200.0 X-COORDINATE (meters) Figure 3.7.14 Vertical velocities along a horizont.a1 profile located at bottom of trenches (Z = 10m) of Case /-FEMWATER. ^ m 4*0 I I I I I I I l I O " K = 10-6m/s x 2,0 _ -- K = 10-10 m/s - sn N E . O.O ~~~ *
- ..~.* l . .* ~ *.*
>- n.~* *~ ~ ':.' ; . .~~ ~,s_1:.,~~ ...?. ,A ..,~. g. O t */ O -2.0 - t 2 ,, ,.4 ., ,^g _.1 1 . / I m \ fi / i [...g+s . 4 > \ \ l s; g . r. , C. ~ 1: -4.0 - v i l' M v I .N v .o a 1 x \ 11 g ->. s = /. . ~ 4 O -6.0 - \ II a I .l v Sl _. v O -8.0 - H E W i i i i i i I I i > 120.0 140.0 160.0 180.0 200.0 20.0 40.0 60.0 80.0 100.0 X-COORDINATE (meters) Figure 3.7.15 Vertical velocities along a horizontal profi;e located at midpoint of trenches (Z = 12.5m) of Case 7-FEMWATER. 0 . ~ 0 _ 0 2 s/ /m s ..{ .$:3.:$- 0 x - I m0 I 0 61 - - - s. 8 1 00 11 y. 0 . I == 2- I 0 d t e KK ,- 6 a c - 1 l o ~.** . R eE 0 l T iA f d I I 0 4 )s oM a\s.8 , t r. 1 r e rE pf l 7 a , :.::- l1 0. te t e ns 1 /. .. I 0 oa zC - 2 (m i s',o - rf 1 oo * ::.:/ h ',v E am ) a,g II sls . 0. T gI n S 1 i - I 0A 0 N l o= - . : . ./. f aZ . 1 I ( s * - D es i e /l 0. R th i c cn 0O i I oe *:... v 8 O l r et v ,'i ,i i i \ n f . l o a
- 0. C - io cp .
I ._ . 1 0X 6 tt r ~ et Va . - 6 0 1 - i I ~. 0 7 . - 4 3 e - r u . g 0 . i ~- 1 0 2 i F n 0 O* 0 . 0 0 O. 4 2 0 2 4 s 8 ^cmow
- eNEv >F ooa_w> yOE4o a4O FEwy m- .
i! i; i ; 'l iI! - !; j ( / / / / // I/ , / i / . ; ' . lI - , 4 f /li I i / lI a n 'II llIl . i-5 E d l 11 A t / . A,-----___s ill o u g = =' /sI e l - QQ s = / - ; o /;l s - /Il E y ,~ _ - - -- > II I o ~ ? s Il g - o i c / l Il 8 o o m
- 5 o
~~~ - -- - II 8 : / ff u Io i - /// o y , """ j , - gx ,, ,.p A,, - - - I g l //j~/s ! I .M // $ o Y n)'j,l I i e i ll /j , % fl ~. / " 1 1 I 9l ,,--# /' 1 _ j {, i i / I / / '/l /I' C i I / o **' ~#/ e s' \ i jf - o,j -f -f N[ i i o o o o ! E d d d j (sJonew) 31YN10HOOO - 2 3-134 ) I / sf/ / / // r / k 9 /I I - I III 5 t / il a ti _ : ! ? )I ? ~ s Ilj// l a ,.-.------- .s ; ,_ y I
- i e
~ lI I - i - -# /ll l .E N c 5 ~~/p l o< 2 s- ,-_____--_d I - g5 / f -eo c I pI 8 3 I t I I , e x / /l l ~ E b o - - - .s I .z / / ,, # - " // I A ~ 4 / ~ / I n s l A a I /\ e k-I /,. -
- ~
I ' , / ',/ l l / , I / - '/ / l / / f/ f / v o o o o k $ Y (sJetew) 31YN108000 - Z 3-135 flowlines at a geologic boundary. Specifically, the velocity profiles indicate vertical downward flow in the trenches for the low conductivity condition. For the high conductivity design condition, the ground-water flow is directed from the upper lef t region of the trenches to the lower right region of the trenches. The dips in the velocity profiles for all three lines of the low conductivity condition, shown in Figures 3.7.14 through 3.7.16, and for the 12.5-meter line of the high conductivity condition, shown in Figure 3.7.15, are believed to be due to the finite-element grid coarseness. Particle trajectories were produced using the same post processor that was described in Case 2 FEMWATER simulations. The predicted particle paths are shown in Figures 3.7.17 and 3.7.18 and are consistent with potential plots (Figs. 3.7.12 and 3.7.13). The initial finite-element grid (Fig. 3.7.10) resulted in some particles becoming trapped in the lower right corner and other particles exiting at the boundary where they are released. Further grid refinement corrected these problems (Fig. 3.7.11). The results presented in this report are those for the final refined grid and are consistent with the problem definition. i l 3-136
- 4.
SUMMARY
AND CONCLUSIONS The NRC staff and contractors participated in the HYDROCOIN study to gain experience and knowledge in using ground-water flow codes against standardized ground-water flow problems relevant to low-level and high-level waste disposal. Specifically, Level 1 of HYDROCOIN was designed to test the accuracy of codes that simulate ground-water flow. The tests of these codes were performed by simulating various types of problems and comparing the results with either anal-ytical results or results from other codes. The Level 1 cases sia.ulated included both linear and nonlinear flow problems. The simulation results of these problems were analyzed to examine (1) the flow theories, (2) applications to nuclear waste problems, and (3) use of post processors to compare velocity fields and contaminant paths. Solution techniques and various gridding schemes used in the numerical codes were also examined to determine how the simulation results are affected by them, 4.1 Research Findings
- 1. A reassuring finding of this study is that the ground-water flow codes used to simulate the Level 1 cases were able to solve the linear problems with a high degree of accuracy. Simulation outputs of the scalar quant-ities of hydraulic head, pressure, and temperature were achieved with little difficulty.
- 2. A finding that was not obvious at the beginning of the study is that, for the types of problems simulated, the finite-difference codes were just as capable as the finite-element codes of simulating ground-water flow, even when the simulations involved complex geometry such as intersecting frac-tures. In Case 1, the project teams were also able to successfully apply a finite-difference code (i.e., USGS-3D), which is based on Cartesian coordinate gridding, to a radial flow problem.
- 3. The only difficulty encountered in simulating linear problems was associated with a region of large contrasts in hydraulic conductivity coupled with an extremely high recharge rate (Case 6). Neither the USGS-3D code nor the SWIFT II code (both finite-difference codes) was able to accurately simulate this region. However, refinement of the SWIFT II grid gave better results, and further refinement may have eliminated this difficulty.
- 4. Simulations of Case 3 (unsaturated flow) and Case 5 (brine flow), which involved numerical solutions using hydraulic parameters that are functions of independent variables (e.g., hydraulic conductivity as a function of saturation), posed great difficulties for all codes tested.
- 5. The nonlinearity in Case 5 (brine flow) arose from the dependence of hydraulic conductivity on ground-water density due to variable salt con-centrations and the pressure distribution. This Sandia project team's 4-1
simulations, which were consistent with the original HYDROCOIN case defini-tion, predicted very little brine movement from the boundary into and through the ground-water flow system. Therefore, the Sandia project team did not encounter difficulties in obtaining a solution. However, other HYDROCOIN project teams, seeking more realistic conditions, modified the case to allow for more brine to enter the system. This modified case, which was a rather simple problem with respect to geometry and boundary conditions, proved to be extremely difficult to simulate. After extensive simulation activities and discussions among the project teams at HYDROCOIN , workshops, solutions to this modified case were finally achieved. How-ever, in the process of arriving at their solutions, other HYDROCOIN pro-ject teams needed to modify their codes and modeling approaches and even f to reformulate the governing equation for coupled brine and fresh-water flow.
- 6. Although the codes provided adequate solutions for linear problems, the subsequent processing of the scalar outputs (e.g., hydraulic head) into vector quantities such as flux or particle trajectories, for purposes of ,
cross-comparisons between project teams' results, caused many difficulties. For example, in Case 1, achieving accurate values for the flux to a well bore proved to be difficult for many codes, among them the USGS-30 code. A refinement in the grid near the well bore helped to reduce the error but did not eliminate it. Analyd s of the simulation results convinced the project teams of the need to examine particle trajectories in detail, which initiated the subsequent particle tracking studies.
- 7. Initially, none of the project teams anticipated any problems with the particle tracking algorithms used in conjunction with the flow codes.
However, it soon became apparent that discrepancies among code results and analytical solutions were caused by these algorithms and not the ground-water flow codes. Some of the problems arose from the manner in which the code itself calculates velocity from the heads that it calculates. In these cases, the error cannot be avoided. For example, the FEMWATER code is designed to produce a smooth velocity field. However, at contacts between zones of different hydraulic conductivity, a discontinuity in the velocity field should be predicted.
- 8. Significant problems were discovered when analyzing the sensitivity of the particle tracking algorithms to discretization. These algorithms are more sensitive to discretization than is the hydraulic-head solution. There-fore, in many simulations, a finer grid is required (both in space and time) to ensure accurate pathlines than is required to ensure accurate hydraulic-head solutions.
- 9. Finally, this report only discusses the NRC and Sandia project teams' results. For a more complete description of all the HYDROCOIN Level 1
results, the reader should refer to the HYDROCOIN secretariat's Level 1 report (Ref. 9). 4.2 Regulatory Implications The research results from the NRC staff and Sandia project teams' simulation activities have shown that the ground-water codes used were successful in code c 4-2 i
verification testing. For linear problems, such as saturated flow conditions, the codes were very successful. This implies that licensing reviews of simula-tions using these codes for saturated sites should concentrate on how the site characterization information was analysed and processed for input into the chosen codes rather than extensive retesting of these codes. However, for non-linear problems, such as unsaturated flow conditions, there were difficulties. Therefore, the codes, as well as the input data, will need to be highly scruti-nized during the license review process. Much more work needs to be done on code verification studies for unsaturated flow models. 4-3 1
REFERENCES
- 1. Michael G. Mcdonald and Arlen W. Harbaugh, "A Modular Three-Dimensional Finite-Difference Ground-Water Flow Model," U.S. Geological Survey Open-File Report 83-875, 528 pp., 1984.
- 2. K. Pruess, "TOUGH User's Guide," Sandia National Laboratories / Lawrence Berkeley Laboratory, NUREG/CR-4645, SAND 86-7104/LBL-20700, 78 pp.,
August 1987.
- 3. S. P. Neuman, R. A. Feddes, and E. Bresler, "Finite Element Simulation of Flow in Saturated-Unsaturated Soils Considering Water Uptake by Plants:
Development of Methods, Tools and Solutions for Unsaturated FlcW," Third Annual Report, Technion-Israel Institute of Technology, Haif a, Israel, 104 pp., 1974.
- 4. L. A. Davis and S. P. Neuman, "Documentation and User's Guide: UNSAT2 -
Variably Saturated Flow Model," Water, Waste and Land, Inc., NUREG/CR-3390, WWL/TM-1791-1, 218 pp., December 1983.
- 5. M. Reeves et al., "Theory and Implementation for SWIFT II, The Sandia Waste-Isolation Flow and Transport Model for Fractured FMdia, Release 4.84,"
Sandia National Laboratories, NUREG/CR-3328, SAND 83-1159, 209 pp., August 1986.
- 6. M. Reeves et al., "Data Input Guide for SWIFT II, The Sandia Waste-Isolation Flow and Transport Model for Fractured Media, Release 4.84,"
Sandia National Laboratories, NUREG/CR-3162, SAND 83-0242, 150 pp., April 1986.
- 7. G. T. Yeh, "Training Course No. 1: The Implementation of FEMWATER (0RNL-5567) Computer Program," Oak Ridge National Laboratory, NUREG/CR-2705, ORNL/TM-8327, 169 pp., June 1982.
- 8. G. E. Fogg and R. K. Senger, "Automatic Generation of Flow Nets with Conventional Ground-water Modeling Algorithms," Ground Water, Vol. 23, No. 3, pp. 336-344, May-June 1985.
- 9. K. Andersson et al., "HYDROCOIN Level 1 Final Report: Verification of Groundwater Flow Models," Swedish Nuclear Power Inspectorate, Box 27106, S-102 52, Stockholm, Sweden.
R-1
l l' APPENDIX A Corrections to Storativity, Leakance, and Transmissivity in Order to Simulate Radial Flow Conditions Storativity Referring to Figure A.1, the volume of the radial block center at r = r g, and extending from r = r g - Dr$/2 to r = r 4 + Drg /2, with thickness Dzj and arc d is: v$,3 = r$Dr.Dz jd The volume of water released from storage in the block by a change in potentio-metric head (Dh) is given by: Dvg g = -S Vs i,j 0h The specific storage value required to correct the model for radial flow condi-tions (S,') is the value required to produce the same release from storage, given the same change in potentiometric head, from the grid block centered at x = x j, and extending from x = x 4 - Dx4 /2 to x = x$ + Dx4/2, with thickness Dzj and width dy. The volume of the block is:
= 0x 4Dzjdy v'$
and the volume of water released from storage in the block in column i by a change in potentiometric head (Dh) is given by: Dv 4 = -S s' V ' i JDh Solving for the value of Ss ' in column i yields: Ss '
- S s"iDr d/Dx j gdy A-1
Setting xj = r g, Dx g = Drg , and with dy numerically equal to d, the corrected specific storage in column i is: S3 ' = S,r$ Where the thickness of layer j is Dz), the storage coefficient is: S'4,j = S,r 40z) (1) Leakance Referring to Figure A.2, the flow from a radial block centered at r = r j, extend-ing from r = r g
- Dr 9/2 to r = r$ + Dr /2, g with thickness Dz) and arc d, in layer j to a similarly defined block of thickness Dz +1 j in layer j+1, may be found from Darcy's law to be:
Q = -r q0r L4jDhd where Dh = potentiometric head change between layer j and layer j+1 in column i, L) = leakance between layers j and j+1, and Q = volumetric flux from layer j to layer j+1 in column i. The leak'ance value required to correct the model for radial flow conditions (L') is the value required to produce the same volumetric flux, given the same change in potentiometric head, between the grid block centered at x = x g, and extending from x = x - Dx j /2 to x = x 4 4 + Dxj /2, with thickness Dzj and width dy in layer j, to a similarly defined block of thickness Dz) 1 in layer j+1. Application of Darcy's law leads to: Q = -Dxj l'j,3Dhdy Solving for L'4,3 yields: L'g,) = r 4Drgl j,j /Dxd dy 9 A-2
i Setting x j = r ,9 Dx) = Dr ,g and with dy.= d:
i,j = r L) 9 since L
j = 2K'/(Dzj + Dz),y) (2) L' 9,j = 2r 4K'/(Dz) + Dz),y) Transmissivity Referring to Figure A.3, the radial transmissivity value that must be approxi-mated is the value for radial flow from rg to r g ). Darcy's law in radial coordinates yields: Q = -KDz j Dh j d/In(r g ,y/rg ) where d = arc of the radial section, Dh = change in potentiometric head between columns i and i+1, Dz) = thickness of layer j, K = hydraulic conductivity, and Q = volumetric flux from column i to column i+1 in layer J. Stated in terms of transmissivity (T): s Q = -Tj 0h 9 d/In(r9 ,y/rg ) The transmissivity value required to correct the model for conditions of radial flow is the value that will produce the same volumetric flux with the same change in potentiometric head. From Darcy's law: Q = -K'9 0Zj Dhdy/(xj ,y - x$) A-3 I l
i where K' is the hydraulic conductivity required by the correction. Stated in terms of transmissivity (T'): Q = -T'$ 3 Dhdy/(x j y - x;) (3) Solving for T' yields: T'4,3 = T)(xg ,y - xj )d/In(r i+1/r g)dy Setting x 4 = r g,x i+1 * "i+1, and with d = dy: T'g,3 = T)(ri +1 - "i)/I"("i+1/ "i) (4) Transmissivity values specified by the modeler (T) are processed by the code to allow for heterogeneous transmissivities. The values actually employed by the code are the harmonic means (M) of the transmissivities specified in adjacent nodes: M4,3 = (Dx 4 + Dxj 7)T$ Tj 7,j/(Dx$Ti+1,j + Dy$ yT g) (5) where Dx, = width, measured in the x direction, of column i, Dx4 7 = width, measured h the x direction, of column i+1, layer j , Mg = harmonic means transmissivity for flow from column i to column i+1 in layer j, T4,) = transmissivity specified by the modeler in column i, layer j, and Tj y = transmissivity specified by the modeler in column i+1, layer J. The modeler must, in order to obtain the correct value for T' (Equation 4 above), specify values of T that have harmonic mean equal to T' Combining Equations 4 and 5 and solving for T$ 1,3: T T' 'i+1,j = Dxj y 'T ' 4,j /(2T' ' 4 31n(r$ 1 /r 4) - Dx$T) (6) A-4
~ _ _ _ _
N~ I D,I .
/ -
l t
,.J~~~,~~~'-~
"i ~>
'e g,g r , 's Dz,
,s' N s
_L s' Ady
, y Dx, ,
Oz, L. s _. f/ s%
#~~
e 8""
,w
~ ~ ~ ~ ~ ~e X j ,z Figure A 1 Block dimensions in radial and Cartesian coordinates used in adjusting storativity.
i A-5 i
rd6 Dri
,f ,
h
/ ~~.
< I
*s ~ Dzj
, ' ,1\ f \ 's n,
> ~~
2 egg,7. g+i
's s Dzj + 3 y ,e s s y
,s ,/ -
/ /
/
/
Dxi I
.L Dz1
'l' ~, jL e$\h,\,
l h
- I Dzj+1
, , L
~~_ ju
~
. s 5,1 \ +.)
K, ~~ ,,-
~'
ta; Figure A.2 Block dimensions in radial and Cartesian coordinates used in adjusting leakance, g A-6 ,
A i
~
rd0 Dri+1 1 I II I /fi+1,zJ~~__ ~s Dzj
- b. .")
Dr, I y y, ? " ' *g ri zj k/'/ '
, M,'~x
,/
v Dxi + 3 NY 3-1 1, I DzI Dxi j ,3 + v[i, .A s,~ y l .s s, i:-
'it Figure A.3 Block dimensions in radial and Cartesian coordinates used in adjusting transmissivity. .
A-7
APPENDIX B Fracture Correction Referring to Figure 8.1, the component of flow parallel to a fracture aligned along a diagonal in the finite-difference grid (from [1,j] to [i+1,j+1]) must j be represented by flow along the two paths [i,j] to [i+1,j] to [i+1,j+1] (path 1) and [i,j] to [i,j+1] to [i+1,j+1] (path 2). I Conservation of mass requires that the total mass flux along paths 1 and 2 l equal the mass flux along the fracture. Assuming constant-density fluid: ' Qf = Qt + Q2 where Q = volumetric flux along the fracture, f Qt = volumetric flux along path 1, and Q = volumetric flux along path 2. Conservation of energy requires that the potential energy change on paths 1 and 2 equal the potential energy change along the fracture: Dhf = Dh t I and , Dhf = Dh 2 ; where Dhf = potential energy change along the fracture, Dht = potential energy change along path 1, and Dhg = potential energy change along path 2. B-1
- _ ._ = _ . - .- - _ - _ . - _
Darcy's law, stated here for the discretized case, with unit dimension normal to the plane of the grid, provides relationships between these terms: 2 Qf = -w ff K (h i+1,P 1) h j,3 /(Dx2 + Dy ), Qi = -DxK'(h9 ,y -h 4 )/Dy = -Dy k' (h )/Dx, and
$ ,y , py - hj ,y Q2 = -Dyk'(h 9 g,7 -h g )/Dx = -DxK'(h4 ,y y y - h$g,7)/Dy.
where Ox = block dimension measured along a row, Dy = block dimension measured along a column, h = potentiometric head, K' = hydraulic conductivity for paths 1 and 2, Kf = hydraulic conductivity of the fracture, and w f = fracture width. Combining these equations and eliminating flux and energy terms yields an expressien for K': K' = w ff K (Dx2 + Dy2)/2DxDy The values of hydraulic conductivity specified by the modeler in the blocks representing the fracture (K) must be calculated to produce K' as the harmonic means of the values specified along paths 1 and 2. The harmonic mean conductivity (assuming constant block size) is given by: K' = 2K' ' K,/(K' ' f + K,) where K,= hydraulic conductivity of the unfractured matrix. Substituting for K' and solving for K yields: K = wff K (Dx 2 + Dy )) 2 K K,(Dx2 + Dy 2)/(4K,0yDx - w ff B-2
b Dx ,l
\s K' Km 2
1,j , 1,j+1 Dy
\ u
\
Km K' 01 y r i + 1,
<J+1 i + 1,j Og Kg Wg Figure B.1 Terms used in adjusting model parameters for fracture simulation.
B-3
APPENDIX C Particle Tracking Algorithm The objective of a particle tracking procedure is to calculate coordinates along a curve passing through a given point where the curve is everywhere parallel to a given velocity field. This calculation can be thought of as a two-stage process:
- 1. Definition of a velocity field over the domain of interest.
- 2. Estimation of a convective transport path from the velocity field.
- 1. Velocity Field Definition An expression for the velocity field valid over the entire domain is not generally available; we more often have samples of velocity components at given points, from which velocities at intermediate positions must be calculated.
Pressures or heads calculated using flow simulation codes can be conveniently used to calculate such velocity samples. The flow models for which the tracking program is designed numerically solve an expression for conservation of volume mass for SWIFT II, subject to given bound-ary and initial conditions. The following schematic represents a typical equa-tion formulated for the steady-state conservation of volume in a discrete system:
"i-1 "i "i+1 X X X R$ .y R (H i-1 - H$)
- R _y = (H$ -H7)+R$
(1) C-1
A series of simultaneous equations similar to Equation 1 relating the physical parameters (R$ 's) and the unknown potentials (H $'s) may be solved by a conven-ient numerical scheme. The resulting potentials will allow expressions such as Equation 1 to be evaluated for internode fluxes that collectively obey volume conservation while satisfying the boundary conditions. These fluxes can then be used to estimate Darcy velocities from the geometrical relationships imbedded in the R$. In 2- and 3-dimensional systems, velocities (Vi) can be calculated by dividing the fluxes by the cross-sectional areas incorporated in the R . $ This operation effectively defines a boundary over which the velocity must have an average value of Vi in order to describe a conservative field. This method of velocity calculation is incorporated in many of the flow simulation codes; a post-processor is required in other cases. Mathematically, flow path calculation is an integration of: ds = Vdt (2) Because the domain of integration is not known a priori, the average velocity, V, must be defined over the entire domain of interest. Interpolation from calculated velocity samples is a convenient way to construct a velocity field for the entire domain. The constructed velocity field should be consistent with the calculations performed in step 1; i.e., V averaged over each face should evaluate to the calculated samples in order to preserve conservation of l volume at the level of discretization. Further, the velocity field should be conservative throughout the model domain; that is, for any closed curve C in the domain, the rate of volumetric accumulation within C, given by: G= 9 n ds (3) ! must be 0, where n is the unit outward normal. By the divergence theorem, G may be written as: C-2
G= ]divVdA=0 (4) provided V is continuous throughout A, where A is the region enclosed by C. For computational simplicity, we would like the velocity within a given region to be a function only of position and the boundary velocities. This method of construction may result in discontinuities of the velocity field acron internal boundaries. As an example, consider the situation shown in Figure C.la. Equa-tion 4 applies in both regions A and B, curve C lies in both regions, and V is discontinuous across the boundary. We wish to constrain V so that G = 0 from Equation 3. Suppose the boundary is used to divide C into two halves (Fig. C.2). Equation 3 applies to both C1 and C2, so that: I(C1) + 1(C2) = 1(Crl) + I(Cb1) + I(Cr2) + I(Cb2) = 0 (5) where I(x) denotes integral (3) performed along curve x. Because V is defined everywhere and Cbl and Cb2 are arbitrarily close to the boundary, we have: 1(C) = 1(Crl) + I(Cr2) (6) Therefore, in order for I(C) = 0 we must have: 1(Cb1) = -I(Cb2) (7) Because C was chosen arbitrarily, V.n =-V.n (8) bl b2 that is, the component of the velocity normal to the boundary must be continuous across the boundary. C-3
In summary, observation of the following restrictions in developing an inter-polation scheme for velocities will provide consistency with the specific mcdel results as well as with the underlying assumption of conservation of volume.
- 1. The average velocity over the boundary separating nodes should equal the j
velocity calculated from model fluxes.
- 2. The velocity should be continuous within a model block.
- 3. The divergence of the velocity should be 0 within a model block.
- 4. The normal component of the velocity should be continuous across block boundaries.
- 2. Trajectory Estimation As mentioned in the previous section, flow path calculation is an integration of:
d5 = V dt When V is generated by interpolation from sampled velocities, it is given as a function of position rather than time: dx = Vx(x,y) dt and (9a) dy = Vy(x,y) dt (9b) from which dx Vx(x,y)
= dt (10a) dy Vy(x,y)
= dt (10b)
Equations 10a and b may be directly integrable, or the integral may have to be approximated by using constant values of V evaluated periodically along the C-4
flow path. These constant values are commonly evaluated at the current position for estimation of a subsequent position. The principal shortcoming of this approach is the inevitable underestimation of flowpath curvature, as illustrated in Figure C.3.
- 3. Discussion of Tracking Procedure Used The average velocities calculated from the model fluxes for each block face q were assumed to apply at each point on the block face. This assumption fulfills condition 1 as well as condition 4. Each velocity component is assumed to vary linearly with displacement along the component direction (Fig. C.4). For a 2-dimensional problem, with the lower lef t corner at (0,0), the velocity is then given by: ;
V = (Vxt + (x/Dx) + (Vx2 - Vxi)) i + (Vyi + (y/Dy) = (Vy2 - Vyi))3 (11) div(V) = [0y(Vx2 - Vxt) + 0x(Vy2 - Vys)]/(0x*Dy) (12) which is identically zero within the block provided Vxt, Vx2, Vyi, and Vy2 satisfy conservation of volume. Equation 10 may be written as:
* = dt (13a)
(Vxt + (x/0 (Vx2 - Vx1)) (Vyi + (y/Dy) (Vy2 - Vyt)) = dt (13b) which are integrated to give: x = (x +Dx(Vxt/(Vx2-Vx3))EXP[Vx2-Vxt)(T-to)/Dx]-Dx(Vxt/(Vx2-Vxt)) o (14a) and y = (yo+0y(Vys/(Vy2-Vys))EXP[Vy2-Vyi)(t-to)/Dy]-Dy(Vy1/Vy2-Vys)) (14b) The analytic expressions for x and y enable the exit point for a block to be calculated directly from the entry point without evaluating intermediate posi-tions. The position within the block may be evaluated periodically to estimate C-5
the path length. The flow path may also be carried partially into the block if the velocities are to ba updated to estimate a transport path for a transient flow field.
'4. Difficulties Encountered in HYDROCOIN Level 1 Problems Flow paths calculated for Case 2 illustrate the difficulty in using a finite-difference approach to simulate systems with a significant flow component along a diagonal of the model grid. Diagonal nodes are not explicitly connected in a finite-difference discret.ization scheme; they are connected indirectly through nodes in the adjacent rows and columns (Fig. C.5). Suppose this arrangement is used to represent a system with diagonal flow components upward and to the right, such as flow along a fracture. Node (i+1,j+1) with head H is downgradient from node (i,j) with head H + D. Nodes (i,j+1) and (i+1,j) must have head values intermediate between H + 0 and H and so must be downgradient from node (i,j).
The gradient in the region between node (1,j) and node (i+1,j) (and between node (i,j) and node (i,j+1)) is away from the idealized fracture, while the gradient in the region between node (i+1,j+1) and node (i+1,j) (and between node (1+1,j+1) and node (i,j+1)) converges on the fracture. Alternating regions of divergence and convergence can be clearly seen in the results for Case 2. l C-6
I Vs VA B A Figure C.1 Path of integration of Equation 3 through a discontinuous velocity field. C-7
Cr, ss NB
%s's
's's
% % C 32 s'ss Cr3 a
*bg ' s *'
B EA A C = C r, + C e, 3 C2= C r, + 0 b 2 l I Figure C.2 Domains of integrations of Equation 4. C-8
^'
/
c. B
\
ftO Aggo O E A Og Figure C.3 Hyoothetical pathline calculated by periodic evaluation of 9. C-9
t 4 x2 V - x1 M X JL Vy, l l e a V x1 Vx2 Al . I _ Vy Vy Vy Vy I I Dx Figure C.4 Linear velocity interpolation. C-10 ' 1
i
\
,2 FRACTURE
[ (1+1,J) t- - - - --- (1 + 1) + 1) > { i IDEALIZED FLOW i/ u l lh H+D _40 _ _ _ _j _
~ 'I+
GRADINNT DIRECTION ('" ~
/ , CALCULATED FLOW PATH Figure C.5 Finite-difference representation of diagonal flow.
C-11
==c n =vetia nea vo , co .. i .po , =wn. m r., ., m,, o s. , v, a. . , -
j , EE E BIBUOGRAPHIC DATA SHEET fiUREG-1249
- u. .. ,.ver.o. v...... Vol. 1 /
,m a.~osa,a a....t. -
,j
"" Model S. ' lations in Support of the Hydrologic Code f
4 rcompari.s (HYDROCOIN) Study ...,,..,o.,co.,,,,,, L .1 1--Code 'erification o,v. f , ..
. .vv oais. February l # 1988
. 0 75 tpf i&5viD 0%,. ....
..u~..oo..uy ,
March fa -n. 1988 n ., o. v, o o.si.e . .o . .. . i.c e eao,e cm.s= = , Division of Engineeri "'"'"'" ^"' Office of Nuclear Regu _ tory Research U.S. Nuclear Regulatory ommission Washington, DC 20555 ,
- i. .,o .o....o....a.,,o,... . o...u.o y u u , i. c , ii. .o..yo.,
Technical Same as 7. above. " " *a c " " o - ~ f August 1985 - January 1988 i, - a.i 1....evi. HYDROCOIN, an international study fo' examining ound-water flow modeling strategies and their influence on safety assessments .f geolo repositories for nuclear waste, has been divided into three levels of effor , Lev 1 evaluates the accuracy of the computer codes by comparing code results to analy 'ca' solutions and intercomparison of code results. Level 2 tests the codes' capabi i es for simulating specific laboratory or field experiments. l.evel 3 addresses the certainty of model results, and the sensi-tivity of those results to changes in ei er\the modeling approach or input parameters. A greater variety of analytic approache and 6cmputer codes were compared in a systematic fashion than would have been possible ue to t large number of participants involved and their corporate efforts. This ortsumm%ar(yesonlytheNRCprojectteams'simu-lation efforts on the Level 1 hen sarking proble . The codes used to sirrulate these seven problems were SWIFT II, FE ATER, UNSAT2, US -30, and TOUGH. In general, linear problems involving scalars suc as hydraulic head w te accurately simulated by both finite-difference and finite ement solution algorithms. Both types of codes produced accurate results even for e .. plex geometries such as intersecting fractures. Difficul-ties were encountered in ving problems that involvedmonlinear effects such as density-driven flow and . saturated flow. Inordertoihilyevaluatetheaccuracyof these codes, post-proc -ing of results using particle tracking algorities and calcu-lating fluxes were er .ir:ed. This proved very valuable byhncovering d'.sagreements among code results even th gh the hydraulic-head solutions had bFen in agreesent.
,. oocw .....s .... ..o oi.c. . on g . . .,.,ge ,yv .
algorithms hydrogeology 4 computer codes performance assessment 1 ground-water Unlimited
, numerical simulations % ,,,,,,,,n,,,,,,,,,,,,
ground-water deling particle trajectories g ,,,,,,,,
. ,en,.. .in w. n 2,.., verification i Unclassified y g...~
HYDROCOIt. Unclassified N , , ,, . . . c. . . c n M(i l(i( l
< s Ul \ * *a u
. . .. u . w .. a n : m u , i . . . . m . , n . , -
UNITED STATES ,,,c,,,, ,,,,,a ,,, ,,,, NUCLEAR REGULATORY COMMISSION PoSTAufbFEESPAIO WASHINGTON, D.C. 20555 "' PERU:7 No. C 47 OFFICIAL BUSINESS PENALTY FOR PRIVATE USE, $300 cllCI\ i i..;1001
- s. .
y $ v. ) 7 T, V ,..k () J E '
- 33. ).n: yemc"q. , a , r', 'e-A -
1 s s
.. c T e 0 9 <; .
Vh r'
,t ,n ,r - * ("p p s
nC 1 s. o% - l ,',y
) 3' .- i Q
- e a,., p41 '4 e}}