Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking Rev 00, ICN 00 ANL-NBS-HS-000001 March 2000 1. PURPOSE The purpose of this Analysis/Model Report (AMR) is to compare transport simulations utilizing particle-tracking methods with simulations using the more rigorous fully coupled advectivedispersive (A-D) approach. This is in accordance with AMR Development Plan for U0155 Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking (CRWMS 1999a). The fully coupled A-D flow and transport simulations incorporate advection, dispersion, sorption, and decay processes. These are compared with results from particle-tracking methods including the method used for the Total System Performance Assessment (TSPA) for the Viability Assessment (VA). This AMR supports the Unsaturated Zone (UZ) Flow and Transport Process Model Report (PMR) as well as other AMRs. In this AMR, two particle-tracking methods are compared with the A-D approach. The results of (1) the Finite Element Heat and Mass (FEHM) particle-tracking code (FEHM, Software Tracking Number (STN): 10031-2.00-00, Version 2.0), which was used for TSPA-VA, and (2) the randomwalk particle-tracking code, Dual Continuum Particle Tracker (DCPT, STN: 10078-1.0-00, Version 1.0), are compared to the results from the code T2R3D (T2R3D, STN: 10006-1.4-00, Version 1.4), a fully coupled A-D numerical code. The constraints and limitations of the results presented here are that the radionuclide breakthrough curves presented should not be considered to be predictions of radionuclide transport in the UZ at Yucca Mountain. The results are for comparison purposes only and the input values used in the comparisons are not necessarily the same as those that will be used in TSPA for Site Recommendation (SR) and License Application (LA). The analysis and simulations, though, do utilize inputs representative of the range of conditions at Yucca Mountain, but these are not necessarily the final properties to be used in the UZ PMR and TSPA-SR/LA. Predictions for the radionuclide breakthrough curve for the UZ for TSPA-SR/LA will be provided in future AMRs and the UZ PMR. It should also be noted that because the effect of radioactive decay would be essentially the same for all of the methods being compared here, it was not necessary to include this process in the comparisons presented here. Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 12 March 2000 INTENTIONALLY LEFT BLANK Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 13 March 2000 2. QUALITY ASSURANCE This AMR was developed in accordance with AP-3.10Q, Analyses and Models. Other applicable DOE Administrative Procedures (APs) and YMP-LBNL Quality Implementing Procedures (QIPs) are identified in AMR Development Plan for U0155 Analysis Comparing Advective- Dispersive Transport Solution to Particle Tracking (CRWMS M&O 1999a). This analysis was evaluated with other related activities in accordance with QAP-2-0, Conduct of Activities, and determined to be quality-affecting and subject to the requirements of the QARD, Quality Assurance Requirements and Description (DOE 1998). This evaluation is documented in Activity Evaluation of M&O Site Investigations (CRWMS M&O 1999b,c). The activity evaluation (per QAP-2-0) completed for performance-assessment activities was also determined to be quality affecting and is documented in Conduct of Performance Assessment (CRWMS M&O 1999d). Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 14 March 2000 INTENTIONALLY LEFT BLANK Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 15 March 2000 3. COMPUTER SOFTWARE AND MODEL USAGE The software codes and routines used in this study are listed in Table 1. These are appropriate for the intended application and were used only within their range of software validation in accordance with AP-SI.1Q, Rev. 1, ICN 0, Software Management. The DCPT (DCPT, STN: 10078-1.0-00, Version 1.0) and FEHM (FEHM, STN: 10031-2.00-00, Version 2.0) codes are used to simulate transport of radionuclides using particle-tracking techniques. T2R3D (T2R3D, STN: 10006-1.4-00, Version 1.4) is used to perform numerical simulations for comparison to the particle-tracking code results. The software code TOUGH2 (TOUGH2, STN: 10007-1.4-00, Version 1.4) is used to generate flow fields for input to the transport codes. The Q-status of these codes and macros is listed in Attachment I and discussed below. Table 1. Table of Software Used in This Analysis TOUGH2 (Version 1.4) and T2R3D (Version 1.4) have been qualified under AP-SI.1Q and were obtained from configuration management. The use of TOUGH2 (Version 1.4) and T2R3D (Version 1.4) prior to obtaining them from configuration management is being evaluated under AP-3.17Q, Impact Reviews, but no impact is anticipated. FEHM (Version 2.0) was qualified prior to the effective date of AP-SI.1Q. It has been reverified and was obtained from configuration management per AP-SI.1Q. DCPT (Version 1.0) is being qualified and a Software Activity Plan for use of unqualified software and copy of the code have been submitted to configuration management per Section 5.12 of AP-SI.1Q, Rev. 2, ICN 2. Software Name Version Software Tracking Number (STN) Computer Type FEHM 2.0 10031-2.00-00 UNIX DCPT 1.0 10078-1.0-00 PC w/Windows 95 T2R3D 1.4 10006-1.4-00 Sun Workstation w/UNIX TOUGH2 1.4 10007-1.4-01 Sun Workstation w/UNIX Routines: T2FEHM2 2.0 ACC: MOL. 19990915.0359 UNIX PROCESS1 1.0 MOL. 19990915.0360 UNIX MAKEPTRK 1.0 MOL. 19990915.0361 UNIX PrepareKDfile 1.0 MOL. 20000127.0120 PC ExtractFlow 1.0 MOL. 20000127.0121 PC ExBT 1.0 MOL. 20000127.0122 PC StatSpatial 1.0 MOL. 20000202.0193 PC Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 16 March 2000 T2FEHM2 (Version 2.0), PROCESS1 (Version 1.0), and MAKEPTRK (Version 1.0) are singleuser software routines qualified per AP-SI.1Q, Rev. 1, ICN 0 and the documentation has been submitted to the Records Processing Center (RPC), the TDMS and is included in Attachement III. PrepareKDfile (Version1.0), ExtractFlow (Version 1.0), ExBT (Version 1.0) and StatSpatial (Version 1.0) were qualified per AP-SI.1Q and the documentation has been submitted to the RPC and is included in Attachment III. T2FEHM2 is a routine written to create FEHM-readable files from TOUGH2 output flow fields. PROCESS1 is a software routine that post-processes the results of the FEHM particle-tracking simulation to provide columns of time versus mass flux and cumulative mass at the water table. MAKEPTRK creates a transport parameter data file for FEHM to read in the particle-tracking simulation. PrepareKDfile is a routine written to create a DCPT-readable file from a TOUGH2 mesh file and T2R3D input file. ExtractFlow is a routine written to create a DCPT-readable file from a TOUGH2 output file. ExBT is a routine written to extract a breakthrough curve from the T2R3D output file. StatSpatial is used to calculate the distribution of particles along a user-specified line based on the DCPT output file. Grids from the UZ Flow and Transport Model are used for comparing these transport codes. Input and output files for this AMR are provided in Attachment II. The commercially-available graphics plotting program Tecplot (Version 7.0) and the plotting portion of KaleidaGraph v.3.09 were also used but are not subject to software qualification assurance requirements. Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 17 March 2000 4. INPUTS 4.1 DATA AND PARAMETERS The input data used in this AMR are summarized in Table 2. The Q-status of these data is provided in Attachment I. Table 2. Input Data The transport simulations comparing T2R3D and the FEHM particle-tracking method use the hydrologic base-case parameter set (DTN: LB971212001254.001) that was used for TSPA-VA. The values used for the sorption coefficients, diffusion coefficients, and dispersivities are given in Section 6.4.3. The precise values of these flow and transport parameters are not considered inputs that require additional verification because the purpose of this analysis is not to document specific transport simulation results, but to compare several transport simulation methodologies for the same transport system. The one-dimensional (1-D) computational grid representing borehole USW SD-9, used in transport simulations comparing the DCPT and FEHM particle-tracking methods to T2R3D results, was obtained from the grid used for TSPA-VA and was used for comparison purposes only. The extraction of this 1-D column is documented in the Scientific Notebook YMP-LBNLGSB- 1.6.3 (pp. 39-40). All input files are listed in Attachment II (DTN: LB990901233129.001 & DTN: SN9908T0581699.001). 4.2 CRITERIA At this time, no specific criteria (e.g., System Description Documents) have been identified as applying to this analysis in project requirements documents. However, this AMR provides information required in specific subparts of the proposed U.S. Nuclear Regulatory Commission rule 10 CFR 63 (see Federal Register for February 22, 1999, 64 FR 8640). It supports the technical basis for methodologies used in performance assessment by comparing outputs with other detailed process-level methodologies (Subpart E, Section 114). DTN Description LB971212001254.001 DKM Basecase Parameter Set for UZ Model, FY97 (Used for FEHM and TOUGH2 Input Parameters) LB997141233129.001 Calibrated Basecase Infiltration 1-D Parameter Set for the UZ Model, FY99. (Used for TOUGH2/DCPT and T2R3D Input Parameters) LB990501233129.004 3-D grid (FY99) used for T2R3D Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 18 March 2000 The DOE interim guidance (Dyer 1999), requiring the use of specified subparts of the proposed NRC high-level waste rule, 10 CFR Part 63 (64 FR 8640), was released after completion of the work documented in this AMR; it has no impact on this work activity. 4.3 CODES AND STANDARDS No specific formally established standards have been identified as applying to this analysis. Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 19 March 2000 5. ASSUMPTIONS This AMR evaluates three numerical simulators by comparing their outputs for radionuclide transport problems, using the input data given in Section 4. The results of these simulations are not to be considered as predictions of transport from a potential nuclear waste repository because the input data are not necessarily the final input values that will be used for TSPA-SR/LA, and because radioactive decay is not included in these simulations. Radioactive decay is handled exactly the same by all simulators and has been ignored because this simplifies the comparisons between the simulation outputs. Any numerical simulator is a simplification or approximation of the physical world. This section lists the principal simplifications and approximations that are used by all the simulators tested in this AMR. It is assumed that these simplifications do not significantly distort the outputs. 5.1 INPUT DATA It is assumed that the input data are sufficiently representative of the conditions at Yucca Mountain that the comparison among the simulators and the findings of this AMR would not change if the input data used for TSPA-SR/LA were not identical to those used here. This assumption is based on several years of evaluations by many investigators and considered to be the only available source of the data. This assumption is used throughout this AMR and requires no further justification. 5.2 TRANSPORT PROCESSES The transport processes included in this analysis are those that were used in TSPA-VA, except for radioactive decay. These are: advection, diffusion and dispersion, and equilibrium sorption of solutes. Radioactive decay has been ignored to facilitate comparisons between the simulation outputs. It is assumed that inclusion of radioactive decay would not significantly affect the comparison among the methods. This assumption is justified because radioactive decay is mathematically simple and is handled identically by all the simulators. This assumption is used throughout Section 6 and requires no further justification. 5.3 DISCRETIZATIONS All standard numerical flow and transport simulators, including those used here, rely upon spatial and temporal discretization, and therefore provide spatially and temporal approximations of the natural system (Wu et al. 1999, pp. 190-193). Also, the methods tested here use dualpermeability grids, described in Section 6 (Wu et al. 1999, pp. 187-188, Doughty 1999, pp. 100-104). It is assumed that the spatial and temporal discretizations, and the appropriate use of dual-permeability grids, do not cause significant errors and do not distort the comparisons among the methods. This assumption is justified by the process of grid development, in which various degrees of grid refinement are tested until further refinement yields little improvement. This assumption requires no further justification. Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 20 March 2000 INTENTIONALLY LEFT BLANK Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 21 March 2000 6. ANALYSIS/MODEL Transport calculations are integral parts to the simulation and prediction of the movement of radionuclides in the UZ. The UZ Model is formulated to rigorously solve both the transport conservation equations and the flow equations using finite-difference techniques. However, as the complexity of the model increases, solving the full transport equations becomes computationally intensive. An alternative approach that is generally less computationally intensive is the use of a particle-tracking method. In addition, compared with finite-element or finite-difference methods, particle-tracking methods usually give better spatial resolution, eliminate numerical dispersion effects, and reduce large truncation errors. However, particle-tracking approaches can vary according to the methods for describing the movement of particles and the assumptions used to determine their interaction with the flow field. Particularly, the exchange of mass between the fractures and matrix in the UZ makes the implementation of particle-tracking approaches more complicated. Therefore, it must be demonstrated that the particle-tracking approach yields acceptable results relative to the more rigorous fully coupled advective-dispersive transport approach. For this AMR, transport simulations are performed with two particle-tracking methods. One is the residence-time-transfer function particle-tracking method of Finite Element Heat and Mass (FEHM, STN: 10031-2.00-00, Version 2.0) that was utilized in the TSPA-VA. The other is the random-walk particle-tracking method used in the Dual Continuum Particle-Tracker (DCPT, STN: 10078-1.0-00, Version 1.0). The FEHM particle-tracking method has been described in the FEHM User’s Manual (FEHM, STN: 10031-2.00-00, Version 2.0) while DCPT is described in this AMR as well as in its software qualification package (DCPT, STN: 10078-1.0-00, Version 1.0). Transport simulations are performed to compare the DCPT to transport problems with analytical solutions and advective-dispersive numerical results using T2R3D (T2R3D, STN: 10006-1.4-00, Version 1.4). Other transport simulations are performed to compare the results using the FEHM particle-tracking method to T2R3D results for a 1-D column. The cumulative breakthrough curves of two radionuclides (one sorbing and one nonsorbing) are compared using the different methods. All test cases used for comparisons with T2R3D simulations use the realistic Yucca Mountain geology from the UZ Model. The results are evaluated for differences between the three approaches, and assessments of the impacts of the differences are provided. Radioactive decay is not included in this comparison analysis because the effect of radioactive decay would be essentially the same for all of the methods being compared here. To facilitate simulation of water flow and solute transport in the fractured porous media, dualpermeability grids are used for all methods in this AMR. In a dual-permeability grid, the problem domain is represented by two overlapping grids, respectively representing the matrix continuum and the fracture continuum. Water or solute can flow between adjacent grid cells in one grid (in the same continuum) or between the two grid cells in different grids that overlap each other (between two continua). This mass transfer between fracture and matrix is a unique feature of transport in fractured porous media. Because the pore-water velocities in the fracture and matrix continua can differ by orders of magnitude, correct simulation of mass transfer between the two continua is one of the key factors that determine the success of a numerical model. In this AMR, the same dual-permeability grid is used for each case, but the approaches used to model the mass Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 22 March 2000 transfer between the fracture continuum and the matrix continuum differ among the three methods. The detailed descriptions are provided in relevant sections (Sections 6.1, 6.2, and 6.3). Key scientific notebooks (with relevant page numbers) used for the analysis described in this AMR are listed in Table 3. Table 3. Scientific Notebooks 6.1 DESCRIPTION OF PARTICLE-TRACKING IN FEHM A complete description of the FEHM particle-tracking model can be found in the Models and Methods Summary for the FEHM software qualification package (FEHM, STN: 10031-2.00-00, Version 2.0) and in AMR U0065 (CRWMS M&O 2000b). Only a brief summary from those documents is provided here. The particle-tracking method in FEHM views the computational domain as an interconnected network of fluid storage volumes. The two steps in the particle-tracking approach for steady-state flow fields are: (1) determine the time a particle spends in a cell, and (2) determine which cell the particle travels to next. The domain can consist of a single-continuum or dual-continua (e.g., fracture plus matrix) representation of the flow field. The time that a particle spends in a cell is a function of the mass of liquid in that cell, the mass flow rates out of that cell into neighboring cells, and the diffusive, dispersive, and sorptive processes within that cell. For advective flow only, the residence time is uniquely defined by the ratio of the mass of liquid in a cell to the sum of the mass flow rates out of that cell. However, dispersive, diffusive, and sorptive processes provide distributions of particle “breakthrough” times for each cell, which are used to determine the effective residence time for a particle in each cell. The standard advection-dispersion equation (with sorption) is used to evaluate the breakthrough times for each cell. If diffusion into an adjacent matrix cell occurs, a onedimensional diffusion equation for transport between the fracture cell and the matrix cell is also included. The analytic solution for diffusion into the matrix cell in the current particle-tracking model assumes an infinite domain. The analytic solutions for the advection-dispersion equation with possible diffusion into a matrix cell yield cumulative, normalized breakthrough concentrations for each cell as a function of time. These curves also represent the cumulative distribution functions for the residence time of a particle that experiences advection, dispersion, diffusion, and sorption in each cell. A random number generator is then used to select a value between 0 and 1, which prescribes a particle LBNL Scientific Notebook ID M&O Scientific Notebook ID Page Numbers YMP-LBNL-GSB-LHH-1 SN-LBNL-SCI-035-VI 83 - 89 YMP-LBNL-GSB-LP-3 SN-LBNL-SCI-155-VI 1 - 105 YMP-LBNL-YSW-2 SN-LBNL-SCI-120-VI 106-108 YMP-LBNL-GSB-1.6.3 SN-LBNL-SCI-085-VI 39-40 Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 23 March 2000 residence time from the cumulative distribution functions. The cumulative distribution function for the residence times is accurately represented with a sufficiently large number of particles that pass through the cell. The probability of a particle traveling to a neighboring cell is proportional to the advective mass flow rate to each neighboring cell. Only outflows from a cell are considered; therefore, the probability of traveling to a cell that has mass flow coming into the current cell is zero. The mass flow rate to an adjacent cell divided by the total mass flow rate out of the current cell is equal to the probability that the particle will travel to that cell. A cumulative distribution function is derived from all the probabilities, and a random number selected between 0 and 1 therefore prescribes the cell to which a particle will travel. Again, a sufficiently large number of particles are used to reproduce the appropriate cumulative distribution function. As described above, the FEHM particle-tracker simulates the advective portion and the diffusive portion of the fracture-matrix mass transfer separately. The advective portion of mass flow between the fracture and the matrix is accounted for by calculating the probability of a particle traveling to a neighboring cell (the matrix cell is treated as one of the neighboring cells to the fracture cell, vise verse). Therefore, the probability of a particle traveling from a fracture cell to a matrix cell is proportional to the advective mass flow rate in the same direction. However, the FEHM particle-tracking algorithm yields only additional residence time (a retardation) for the particles in the fracture that experience diffusive mass flow from fracture into matrix, but the particles do not actually get transported into the matrix. The additional residence time is calculated based on an analytical solution for a single fracture system (Tang et al. 1981, pp. 555-564). This model implies that the particles diffusing into the matrix cannot move vertically unless they first diffuse back to the fractures. Though FEHM particle tracker has this capability, radioactive decay is not used in this comparison because the effect of radioactive decay would be essentially the same for all the methods being compared here. 6.2 DESCRIPTION OF PARTICLE-TRACKER DCPT 6.2.1 General Approaches and Overall Structures The random-walk particle tracker DCPT describes the history of individual particles instead of focusing on fixed points of space. It uses the Lagrangian point of view, not the Eulerian point of view. The movement of a plume is described as a sum of the movements of individual particles. The coordinates of a moving particle are represented as functions of time (Bear 1972, p. 70, Equation 4.1.18): (Eq. 1) where X and . are the vectors that describe the positions of the particle at time t and some initial time (e.g., t = 0), respectively. Note that X is the dependent variable (vector) in Equation 1 while the function includes factors such as velocity, dispersion coefficient, and adsorption parameters. ) , ( t . X X = Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 24 March 2000 The velocity, dispersion coefficients, and adsorption parameters are generally functions of space. These data are provided as tables of values on a discretized space (e.g., a grid). DCPT transforms these fixed-space values in the Eulerian point of view into the parameters of Equation (1) in the moving-particle (Lagrangian) point of view. Because the whole domain is discretized into subdomains or grid cells, the velocity field, or other fields of parameters, can also be disassembled in the same way. Cells are the basic units of a domain. Each cell has two sets of parameters, each of them corresponding to one continuum, and one set of parameters that defines the interactions between two continua. In dual-continua media (i.e., fractured-porous rock), a particle will travel either in the fracture continuum or in the matrix continuum, two overlapping continua often with very different velocities and parameters. The random switch between the fracture and the matrix is governed by a particle-transfer probability that should be consistent with the mass flow between two continua within that cell. The object-oriented-program approach is used in developing the DCPT. Two major objects are used in DCPT. One is called CELL, which has all the information of the continua (e.g., the geometry, local velocities, dispersion coefficient tensor, and other parameters for both fracture and matrix). The other is called PARTICLE, which has properties describing the current status of a particle including the current time, the current XYZ position, the current cell, and the current continuum (fracture or matrix). Therefore, the major algorithm of particle tracking for a given particle and a given time step can be summarized as below: 1. Calculate the displacement that the particle will take during the time step based on the current status of the particle (see Section 6.2.2); 2. Determine whether the path of the particle intersects with any face of the current cell; if it does not, go to Step (3); otherwise, use the intersection point as the new location of the particle, reduce the time step accordingly, and get the neighboring cell ID; 3. Determine whether the particle will switch to the other continuum at the next time step (see Section 6.2.4); 4. Update the status of the particle with the results of Steps 2 and 3; 5. Check whether the particle has exited through the domain boundary or the specified maximum time has been passed; if yes, finish the simulation of this particle, otherwise go to Step 1. In short, DCPT simulates the random walk of particles in a continuous space with discretized continua (cell based), but uses the particle-transfer probability to control which continuum a particle will travel in at a particular time. The following are some details of the approaches used in DCPT. Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 25 March 2000 6.2.2 Calculation of Particle Displacement The new location of a particle after a time step .t is a random vector and can be calculated as (LaBolle et al. 1996, Equation 3, p. 584, symbolically replacing Xp and .w with X and , respectively) (Eq. 2) The drift term A (see LaBolle et al. 1996, Equation 10, p. 584) is approximated to be the local pore velocity V. The tensor B and its transpose BT are given by BBT = 2D where D is the local dispersion coefficient tensor. W is a random vector, each component of which observes the N(0,1) distribution. For simplicity, two additional terms in drift term A related to the divergence of D and the gradient of the volumetric water content are neglected. As shown in Sections 6.4.3 and 6.4.4, this approximation is acceptable for the advection-dominant transport processes in a steadystate flow field, such as was used for TSPA-VA. For a particle, the mean displacement vector is V.t while the variance tensor is 2D.t in a given continuum. Whether the properties of the fracture or those of the matrix are used in Equation 2 depends on which continuum the particle travels in. 6.2.3 Sorption and Decay For a reactive solute, only a portion of particles are mobile as described by Equation 2 with the remainder being sorbed. The probability, Pr, of a particle being in fluid can be defined as: (Eq. 3) where Kd , f, and .R are the sorption distribution coefficient (m3/kg), the porosity (m3/m3) and the rock density (kg/m3) of the particular continuum, respectively; and . is the volumetric water content. In terms of implementing the sorption process in particle tracking, we can take Pr in a deterministic way, as the percentage of the total mass of a moving particle. Therefore, the effective displacement of the particle will be Pr times the original displacement, which can be implemented by simply multiplying A and B in Equation 2 by Pr . To simulate the radioactive decay, the mass of each particle, Mp, is calculated as a function of time, t: (Eq. 4) t W . t BW t A t X t t X . + . + = . + ) ( ) ( K P R d r ) 1 ( - + = . f . . 5 . 0 / 2 ) 0 ( ) ( t t p p M t M - = Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 26 March 2000 where t0.5 is the half-life. Though DCPT has this capability, radioactive decay is not used in this comparison because the effect of radioactive decay would be essentially the same for all of the methods being compared here. 6.2.4. Particle-Transfer Probability: Mass Transfer between Fracture and Matrix The mass-transfer process between fractures and matrix is simulated by random particle exchanges between two continua as controlled by the particle-transfer probabilities of either fracture-to-matrix or matrix-to-fracture progression, as described in Step 3 in Section 6.2.1. For other variables such as velocity and dispersion coefficients, grid cells are used in DCPT as the basic units for evaluating the particle-transfer probability. For each net mass flow between two continua in the fixed-space Eulerian point of view, there are two corresponding particle-transfer probabilities in the moving-particle Lagrangian point of view. One is the particle-transfer probability of particles from fracture to matrix, and the other is that from matrix to fracture. The challenge is how to transform correctly the net mass flow in the Eulerian point of view into two separate particle-transfer probabilities in the Lagrangian point of view. In the following derivation, we focus on the particle-transfer probability Pfm from a fracture to the matrix. The other probability Pmf can be similarly derived. If the particles in the fracture continuum of a given grid cell at t = 0 have mass M0, and the fraction of them that enter into the matrix continuum during the time interval (0, t) have mass Mfm, the particle-transfer probability Pfm can be defined as: (Eq. 5) For a single particle in the fracture at t = 0, Pfm is the probability at which it will be in the matrix at time t. Mfm is directly proportional to the mass flow from fracture to matrix. For a given grid cell, the net solute mass J transferred from fracture to matrix during a small time interval dt through a small area of the interface dA is: (Eq. 6) where qfm is the water flux (L/T) between fracture and matrix and is the normal vector of the interface and points from fracture to matrix. This being the case, only one of the two advection terms in Equation 6 will take effect, depending on the direction of water flow. C is concentration and D is the dispersion coefficient specifically for the fracture-matrix interaction. A and t are fracture-matrix interfacial area and time, respectively. The variable s is the distance away from the fracture-matrix interface (s = 0 at the interface). Because in reality the detailed geometry of the interface and those variables defined on the interface are not available, it is not practical to derive a formulation to calculate the total mass flow between fracture and matrix (even in cell-scale) 0 M M P fm fm = dt dA s C D - ) , C q (- - ) , C q ( = J d s m m fm f fm fm .. . .. . . . =0 0 max 0 max n& Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 27 March 2000 without some simplifications. In DCPT, a lumping approach similar to that in T2R3D is used to estimate the net mass transfer between the fracture and the matrix at the grid-cell scale. Assuming (as in Section 5.3) that all dependent variables and the parameter D can be used in a sense of average values within the grid cell or over the interface, we can get the net mass transfer during the time interval (0, t) by integration of Equation 6 over the whole interface area: (Eq. 7) where S is the characteristic distance of the fracture-matrix system proportional to the fracture spacing (e.g., 1/6 of fracture spacing depending on the assumptions of the fracture network). Qfm is the net water flow rate (M/T) between fracture and matrix. Its value is positive if the mass flows from fracture to matrix. Note that t is a particular time, i.e., the end of a time step, while t is the variable of integration. Equation 7 can be rewritten as: (Eq. 8) Ffm is the effective flow rate from the fracture to the matrix while Fmf is the effective flow rate from the matrix to the fracture. Equation 8 states that the net mass transfer from fracture to matrix is the total mass flow from fracture to matrix less the total mass flow from matrix to fracture. The first term on the right hand of Equation 8 is the mass flow from fracture to matrix during the time interval (0, t). However, it is not Mfm in Equation 5 because Cf includes not only the particles that are in the fracture at t = 0, but also those particles that enter the fracture during (0, t) from either the matrix or other neighboring blocks. In other words, if CE is the concentration of the particles that entered the continuum during (0, t) and CNE is the concentration of the particles already in the continuum at t = 0, we can split the first term on the right hand of Equation 8 as: (Eq. 9) t d A] ) C - C ( S D + ) , C Q (- - ) , C (Q [ = J m f m fm f fm t fm 0 max 0 max 0 . t t t t d C F d C F = d C A ] S D + ) , q (- [ - d C A ] S D + ) , q ( [ = J m t 0 mf f t 0 fm m fm t 0 f fm t 0 fm + . . . . 0 max 0 max t t t d C F d C F d C F f E N t 0 fm f E t 0 fm f t 0 fm . . . + = Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 28 March 2000 Equation 9 simply states that only a portion of the mass flow from the fracture to the matrix consists of the particles that were initially in the fracture (the second term). Only this portion is needed to calculate the probability of the particles, which are in the fracture at time zero, being transferred from the fracture to the matrix. CNE decreases with time t monotonically. In what follows, we will only discuss the particle-transfer probability corresponding to mass transfer from the fracture to the matrix and drop the subscript “f” and the superscript “NE” for simplicity. The derivation of the particle-transfer probability corresponding to mass transfer from the matrix to the fracture is similar and will not be repeated. For the particles that are in the fracture of a given grid cell at time = 0, we can write a mass conservation equation for those particular particles as follows: (Eq. 10) where V0 and VT are the volume of water in the fracture and the total volume of the fracture within the cell, respectively. .b is the bulk density. Fout is defined as follows: (Eq. 11) where M is the number of interfaces between the grid cell and other neighboring grid cells. Qi (outward positive), Di , Ai , and Si are water flow rate, dispersion coefficient, area, and distance between the neighboring nodes of the i-th interface, respectively. The left-hand side of Equation 10 is the initial mass of the particles in the fracture of the given grid cell, while the first term of the right-hand side is the mass of the particles that still stay there at time t. The second term and the third term on the right-hand side of Equation 10 are the mass of particles flowing into the matrix of this cell and into the fracture continuum of other neighboring cells, respectively. Taking derivatives on both sides of Equation 10 with respect to t, we have a first-order ordinary differential equation: (Eq. 12) with initial condition , where (Eq. 13) is the characteristic time in the system, which indicates how slowly the mass will be replaced for a given cell. t t t t . . d ) ( C F + d ) ( C F + V K V (t) C = V K V ) ( C 0 out 0 fm T b d 0 T b d 0 . . + + ) ( ) ( 0 .. . .. . . S A D + ) , Q ( = F i i i i M =1 i out 0 max 0 1 = (t) C t + dt (t) dC c C = (0) C 0 out fm T R d c F F V K V t + - + = . f) 1 ( 0 Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 29 March 2000 The solution of Equation 12 is readily obtained as: (Eq. 14) Therefore, the probability of a particle being transferred from fracture to matrix during (0, t) can be calculated as: (Eq. 15) Similarly, the particle-transfer probability corresponding to the mass flow from matrix to fracture, Pmf, can be calculated based on Equation 15 by replacing Ffm with Fmf and using Fout and tc of the matrix continuum. 6.2.5 Adaptive Time Steps Particle-tracking time steps used in DCPT are adaptive to the local flow field, cell size, and other transport parameters. For a given type of solute, each cell has two time steps corresponding to the fracture continuum and the matrix continuum, respectively. For either continuum of a cell, the time step is calculated as follows: (Eq. 16) where .xy and |Vxy| are the lateral size of the cell and the magnitude of the lateral component of the pore velocity within the cell, respectively, and .z and |Vz| are the height of the cell and the magnitude of the vertical component of the pore velocity, respectively. Equation 16 limits the time step so that the Courant Number (defined by for the horizontal or vertical direction) is equal to or less than 0.05. This limit is sufficient for the proper accuracy of the explicit approaches used in DCPT by establishing an adequate temporal resolution regarding particle transfer between fracture and matrix. If sorption exists, effective velocities are used in Equation 16 by multiplying the original pore velocities by the factor Pr (see Section 6.2.2). 6.3 DESCRIPTION OF ADVECTIVE-DISPERSIVE SOLUTIONS WITH T2R3D As a member of the TOUGH2 family of codes, T2R3D (Wu et al. 1996, pp. 8-32) provides a capability for modeling liquid or gas tracer or radionuclide transport in multiphase and nonisothermal flow systems. In particular, T2R3D can be used to simulate tracer transport in a complex, heterogeneous fractured rock using a general, irregular 3-D grid. In addition to incorporation of a full dispersion tensor in evaluating dispersive tracer transport, the code takes into account linear adsorption and first-order decay effects. The model formulation and numerical ) t / (-t C = (t) C c 0 exp [ ] ) t / t (- - F + F F = V K V C d ) ( C F = M M = P c fm out fm T b d 0 0 t 0 fm 0 fm fm exp 1 ) ( . t t + . ) 25 . 0 , | | 05 . 0 , | | 05 . 0 min( c t Vz z Vxy xy t . . = . Co .t Vxy .xy ------------------or.t Vz .z --------------- = Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 30 March 2000 scheme make it easy to include many transport mechanisms, such as nonadsorption, multidecay chains, or thermal/mechanical effects. T2R3D is built on the framework of the TOUGH2 code (Pruess 1991, pp. 5-9). The basic mass and thermal-energy balance equations for threecomponent fluid and heat solved by T2R3D are similar in form to those for the standard TOUGH2 EOS3 module (Pruess 1991, pp. 21-23). The integral finite-difference method and a first-order, backward finite-difference scheme are used for spatial and temporal discretization, respectively. The tracer transport mechanisms include molecular diffusion and hydrodynamic dispersion in the liquid or gaseous phase, in addition to advection terms. First-order decay is taken into account, and adsorption of a tracer on rock matrix and/or fractures is described by an equilibrium isotherm with a constant sorption distribution coefficient. The model formulation considers advection/dispersion transport processes of a liquid or gas tracer with a full-dispersion tensor, in a heterogeneous geological system. The grid can be either regular or irregular. In addition to advection terms for the tracer transport, the dispersive and diffusive mass flux, , is described by: (Eq. 17) where is the combined diffusion-dispersion tensor accounting for both molecular diffusion and hydrodynamic dispersion; is fluid density; and Xß is mass fraction of the tracer in phase ß (ß = liquid or gas) and superscript . represents the solute component. A general dispersion model for 3-D tracer transport in T2R3D is: (for ß=liquid or gas) (Eq. 18) where aT and aL are the transverse and longitudinal dispersivities, respectively; vß is the Darcy velocity vector of phase ß through fractures or matrix; t is the tortuosity of the medium; dm is the molecular diffusion coefficient in phase ß; and dij is the Kronecker delta function (dij = 1 for i = j, and dij = 0 for i . j). One of the key issues in implementing the general 3-D dispersion tensor of Equation 18 is how to interpolate velocity fields for determining the dispersion tensor. The averaging or weighting scheme used to evaluate a velocity vector at the interfaces between cells is called “projected area weighting method” (Wu and Pruess 1998, pp. 139-146). In this method, we calculate a velocity component, vn,i , of the velocity vector of cell n by the summation of the flow components of all local connection vectors in the same direction, weighted by the projected area in that direction: (for i = x,y,z) (Eq. 19) FD . ( ) FD . ( ) = -.ß D • .Xß . ( ) D ß . D = aT vß dij+ aL - aT ( )vßvß vß +fSßtdmdij . . = m i nm m i nm nm i n A v A v ) ( ) ( ) ( , n n ni Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 31 March 2000 where M is the total number of connections between cell n and all its neighboring cells, vnm is the flux along the connection to cell m in the local coordinate system, and ni are the directional cosines of connections. The velocity vector v at the interface of cells n and m is then evaluated by harmonic weighting to preserve total transit time for solute transport traveling between the two cells. The mass fraction gradient of the tracer/radionuclide is evaluated at the interface between cells n and m as: (Eq. 20) with (Eq. 21) The net mass flux of diffusion and dispersion of a tracer/radionuclide along the connection of cells n and m is determined by Equation 17. In the above calculation, the connection to the overlapping cell in the other continuum is excluded because it involves the mass transfer between the fracture continuum and the matrix continuum. This mass transfer is treated as a 1-D advection-dispersion transport process and added to the mass conservation equation of each cell. Though T2R3D has this capability, radioactive decay is not used in this comparison because the effect of radioactive decay would be essentially the same for all of the methods being compared here. 6.4 COMPARISONS OF FEHM AND DCPT WITH T2R3D In this section, DCPT is first compared with analytical solutions for 1-D and 2-D cases in Sections 6.4.1 and 6.4.2. The FEHM particle-tracking code has been previously compared to analytical solutions as part of its qualification (FEHM, STN:10031-2.00-00, Version 2.0). Both the DCPT and FEHM particle tracker are compared with T2R3D for a 1-D case in Section 6.4.3, and then the DCPT is compared with T2R3D for the full 3-D case of the Yucca Mountain UZ Model. Again, the FEHM particle-tracking code was used in TSPA-VA. Radioactive decay is not included in all of the comparisons discussed below because the effect of radioactive decay would be essentially the same for all of the methods being compared here. The numerical values of physical and geometric parameters for the selected test cases were chosen to provide reasonable representations of the real-world scales and properties appropriate to the flow and transport process under consideration. .Xnm . ( ) = nx.Xnm . ( ),ny.Xnm . ( ), nz.Xnm . ( ) ( ) .Xnm . ( ) = Xm . ( ) - Xn . ( ) Dm + Dn Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 32 March 2000 6.4.1 1-D Cases Comparing Analytical Solutions with DCPT The first test case is 1-D solute transport in a fractured-porous medium with parallel fractures, for which the particle-transfer-probability approach and the sorption model of DCPT can be tested against an analytical solution (Sudicky and Frind 1982, pp. 1634-1642). The analytical solution is based on the approximation that solute transport between fractures and matrix occurs through matrix diffusion in the direction perpendicular to the fracture only. Matrix advection and diffusion in the direction along the fracture is ignored. Furthermore, the initial solute concentration is zero in the system, and the concentration at inlets of fractures (z = 0) is constant for time t > 0. The diffusion/dispersion in the fractures is also ignored. The rationale for the parameters shown in Table 4 is documented in Scientific Notebook YMP-LBNL-GSB-LHH-1, pp 83-89. For this case, the integral of the breakthrough curve corresponding to a pulse input is equivalent to the breakthrough curve corresponding to a constant concentration input, which is the solution in Sudicky and Frind (1982, pp. 1634-1642). Table 4. Parameters Used in Transport Problem in a Parallel Fracture System Figure 1 shows the cumulative mass fraction (the integral of the DCPT breakthrough curve divided by the initial mass released) flowing out from the fracture at the outlet as a function of time. The results from DCPT are similar to those for the analytical solution. This implies that the particle-transfer-probability approach (used in DCPT) of diffusive mass exchange between fracture and matrix is representative for this transient case. Note that the fracture spacing is 1.0 meter, which is within the range of the fracture spacing in the unsaturated zone of the Yucca Mountain site. The CPU time used in simulation by DCPT on a PC (Pentium II 300) is about 10 seconds, excluding the time used for reading/writing files. Filenames are given in Attachment II. Parameter Value Molecular diffusion coefficient 2.5×10-11 m2/s Fracture spacing 1.0 m Retardation factor 30 Velocity in fracture 1.1574×10-5 m/s Grid spacing 0.5 m Matrix volume per cell 0.25 m3 Fracture volume per cell 0.5×10-4 Fracture/matrix interface area 0.5 m2 Domain length 36.75 m Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 33 March 2000 Based on data submitted with this AMR under DTN: LB990901233129.001 Figure 1. Comparison between DCPT and the Analytical Solution for a Parallel Fracture System 6.4.2 2-D Cases Comparing Analytical Solutions with DCPT The second test case is 2-D solute transport in a porous medium (no fractures) with a dispersion tensor, for which the advection and dispersion model of DCPT can be tested against an analytical solution. Table 5 shows the case specifications with all parameters dimensionless; the rationale for these parameters is documented in Scientific Notebook YMP-LBNL-GSB-LP-3, pp. 1-105. 0 0.2 0.4 0.6 0.8 1 Relative Mass TIME (years) DCPT Analytical 10 2 10 3 10 4 10 5 Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 34 March 2000 Table 5. Parameters of the 2-D Case As defined in Table 5, the scenario is a simultaneous injection of mass at time zero in a 2-D uniform flow field (flow in z-direction only). If M is the mass of a point source injected at (x0, z0) at t = 0, the concentration distribution in the field at any later time is given by (Bear 1972, p. 633, Equation 10.6.34, symbolically replacing n, ., ., , , and q/n with f, z0, x0, Dz, Dx, and Vz): (Eq. 22) where Dx (= aTVx) and Dz (= aLVz) are dispersion coefficients corresponding to x-direction and z-direction, respectively, and f is porosity. The problem is actually simulated with DCPT as a 3-D transport problem with no discretization in the y-direction. Solutes are released at time zero in the form of a point pulse source (M/f = 1). At t = 10, the relative concentration along x = x0 (= 10.25) is calculated within the specific slice. Figure 2 compares these results with the analytical solution. The concentration distribution simulated by DCPT is consistent with the analytical solution. This consistency indicates that DCPT properly incorporated the dispersion tensor. All values are dimensionless. Two million particles were used in the simulation (Figure 2), and the CPU time used is about 10 minutes on a PC with a Pentium II 300 processor. All input and output filenames are given in Attachment II, Section 1. Parameters Value Domain dimension (x, y, z) 20.5×20.5×30.5 Pore velocity Vx = Vy = 0, Vz = 1 Dispersivity aL = 0.05, aT = 0.01 Diffusion coefficient 0.0 Grid spacing .x = .z = 0.5; .y = 20.5 Plume location at t = 0 x = 10.25, y = 10.25, and z = 0.0 Monitoring location at t = 10 x = 10.25 Monitoring resolutions dx = 0.02 and dz = 0.01 D' D ' .. . .. . - - - - - f = t D x x t D t V z z D D t / t z x C x z z z x 4 ) ( 4 ) ( exp 4 M ) , , ( 2 0 2 0 p Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 35 March 2000 Based on data submitted with this AMR under DTN: LB990901233129.001 Figure 2. Comparison between DCPT and the Analytical Solution for a 2-D Transport Problem with Dispersion Tensors. 6.4.3 Comparison of Numerical Solutions (T2R3D) with FEHM and DCPT for 1-D Cases Analytical solutions are only available for the simplified cases (e.g., no advective flow between fracture and matrix). In those cases, the critical features of the particle-tracking models cannot be fully tested. A more realistic one-dimensional transport problem is thus designed to further test the capabilities of the particle-tracking models against the numerical solutions provided by T2R3D, mainly focusing on simulations of the fracture-matrix mass exchange and sorption processes. The case involves a column near borehole USW SD-9 extracted from the 1997 3-D model of the Yucca Mountain site (DTN: LB971212001254.001). The radionuclides are released at the simulated repository horizon at time zero as a pulse. A steady-state flow field is assumed and determined using TOUGH2 version 1.4. The transport parameters used in simulations are shown in Table 6. A total of 2,000 particles are used in DCPT simulation. The CPU time used is about 10 seconds for both DCPT and T2R3D, with DCPT executed on a Pentium II PC and T2R3D on a DEC ALPHA. A total of 27 cells are used. 0 0.1 0.2 0.3 0.4 C -10 0 10 20 30 Analytical DCPT Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 36 March 2000 Table 6. Parameters Used for 1-D Radionuclide Transport In this case, significant mass flow occurs as a result of advection and dispersion between fracture and matrix. Figure 3 shows the cumulative mass fraction at the water table versus time. The cumulative mass fraction is defined as the cumulative mass flowing out to groundwater divided by the total mass released at the repository horizon. In both cases, the results are very similar except that the DCPT shows fewer numerical mixing effects than the T2R3D. The good agreements between the DCPT and the T2R3D show that the approximation (A = V) in Equation 2 is acceptable for the UZ transport of radionuclides in the Yucca Mountain site. The input and output files and the process for performing DCPT simulations are provided in Attachment II, Section 1. The comparison of the FEHM particle-tracking simulations with the advective-dispersive (A-D) transport simulations of T2R3D consisted of a single 1-D flow simulation along borehole USW SD-9, with four subsequent transport simulations. The details regarding input and ouput files and use of software macros for this part of the analysis are provided in Attachment II, Section 2 (DTN: SN9908T0581699.001). The four transport simulations are detailed in Table 7. Parameter Value Molecular diffusion coefficient of technetium 3.2×10-11 m2/s Molecular diffusion coefficient of neptunium 1.6×10-10 m2/s Fracture longitudinal and transverse dispersivity 20 m and 0 Matrix longitudinal and transverse dispersivity 0 Fracture-matrix dispersivity 0 Fracture and matrix tortuosities 0.7 and 0.7 Temperature 25 °C Sorption distribution coefficients of technetium Zero in both fracture and matrix Sorption distribution coefficients of neptunium Zero in fracture and matrix of TCw, PTn, TSw units; 4.0 × 10-3 m3/kg and 1.0 × 10-3 m3/kg in matrix of zeolitic rock and vitric rock in CHn unit, respectively. Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 37 March 2000 Based on data submitted with this AMR under DTN: LB990901233129.001 Figure 3. Comparison between DCPT and T2R3D for 1-D Radionuclide Transport. (a) Technetium, (b) Neptunium 0 0.2 0.4 0.6 0.8 1 Relative Mass 10-1 1 Time (years) DCPT T2R3D 10 10 2 10 3 10 4 10 5 0 0.2 0.4 0.6 0.8 1 Relative Mass Time (years) DCPT T2R3D 10 10 2 10 3 10 4 10 5 1 0 6 (a) (b) Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 38 March 2000 Table 7. Four Transport Simulations Used in FEHM vs. T2R3D Comparison These four simulations consider the transport of two radionuclides, Tc and Np, under conditions with and without matrix diffusion and fracture dispersivity. The Np is assumed to sorb within the matrix, but the Tc does not, and in no case does sorption occur along the fracture. The sorption distribution coefficients for the matrix of different geological units are given in Table 6. Particularly, the thickness of the vitric rock (Kd=1×10-3 m3/kg) and the zeolithic rock (Kd=4×10-3 m3/kg) in CHn unit is 46.63m and 103.16 m, respectively. A finite amount of radionuclides was released at a cell near the potential repository elevation at 1063 m, and the transport simulation was run for one million years, with the cumulative breakthrough (normalized) of the radionuclide plotted as a function of time for each of the four cases. CPU time used for each simulation using FEHM particle tracker is less than 1 minute on a SUN ULTRA SPARC machine. The results of the simulations are shown in Figures 4 and 5. Figure 4 shows the cumulative normalized breakthrough for technetium. The solid lines are the results of T2R3D; the dashed lines are the results of FEHM. Both cases, with and without matrix diffusion and sorption, are shown. Results for T2R3D and FEHM are very similar for the advection-only case. The FEHM particle-tracking results show sharper breakthrough fronts at the water table. This is reasonable because the particle-tracking method reduces the numerical dispersion associated with finitedifference and finite-element methods as used in T2R3D. The initial breakthrough at around one year is a result of advective transport of technetium through the fractures. Both methods show that over 60% of technetium reaches the water table in the case without matrix diffusion and dispersion. The second major breakthrough in the case without dispersion or matrix diffusion occurs around 10,000 years. This breakthrough represents the transport through the matrix continuum between the repository and water table. Radionuclide Simulated Molecular Diffusion (m2/s) Distribution Coefficient, Kd (m3/kg) in Vitric, Zeolitic Fracture Dispersivity (m) Technetium (Tc) 3.2x10-11 0, 0 20 Technetium (Tc) 0 0, 0 0 Neptunium (Np) 1.6x10-10 1.0 x 10-3, 4.0 x 10-3 20 Neptunium (Np) 0 1.0 x 10-3, 4.0 x 10-3 0 Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 39 March 2000 Based on data submitted with this AMR under DTN: SN9908T0581699.001 NOTE: Diffusion Coefficient = 3.2 x 10-11 m2/s, dispersivity = 20 m (fractures only), Kd = 0 m3/kg Figure 4. Normalized Cumulative Breakthrough of Technetium at the Water Table for FEHM and T2R3D Based on data submitted with this AMR under DTN: SN9908T0581699.001 NOTE: Diffusion Coefficient = 1.6 x10 -10 m2/s, dispersivity = 20 m (fractures only), Kd = 4.0 x 10-3 m3/kg (zeolitic), Kd= 1.0 x 10-3 m3/kg (vitric) Figure 5. Normalized Cumulative Breakthrough of Neptunium at the Water Table for FEHM and T2R3D 0 0.2 0.4 0.6 0.8 1 10-1 10 0 101 102 10 3 104 105 106 Normalized Cumulative Breakthrough at Water Table Time ( years) no diffusion or dispersion with diffusion and dispersio n A-D direct solution (T2R 3D ) Particle- tracking solution (FEHM) 0 0.2 0.4 0.6 0.8 1 10-1 100 101 102 103 104 105 106 Normalized Cumulative Breakthrough at Water Table Ti me (years) no diffusion or dispersion A-D direct solution (T2R3D) Particle tracking solution (FEHM) with diffusion and dispersion Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 40 March 2000 In Figure 4, the technetium transport with matrix diffusion is significantly different in the T2R3D and FEHM simulations. The FEHM results indicate that the initial breakthrough in the fractures is “smeared,” but the asymptotic plateau is the same as the plateau for the case with no matrix diffusion (~65%). The reason is that the implementation of the diffusive mass flow from the fracture to the matrix in the FEHM particle-tracking algorithm yields additional residence time (a retardation) for the particles in the fractures that experience diffusive mass flow from fracture into matrix, but the particles do not actually get transported into the matrix. As a consequence, the shape of the initial breakthrough for the FEHM simulation yields the same plateau as the case with no matrix diffusion. This approach is based on an analytical solution by Tang et al. (1981, pp. 555-564), which assumes that diffusive mass flow within the matrix only occurs in the direction perpendicular to the fracture. Hence, the particles can leave the flow system via the fracture only. However, like the DCPT and the T2R3D, the FEHM particle tracker implements the advective mass flow in both fracture and matrix continua, which allows the particles to transport to the water table through either fracture or matrix. As a result, the FEHM particle tracker gives very similar results to those of T2R3D for the advection-only cases. The T2R3D results, in contrast, show an initial breakthrough that has less than 30% of technetium arriving at the water table through the fractures. Recall that without the diffusion between fracture and matrix (which is controlled by the matrix diffusion), over 60% of the technetium arrived at the water table through the fractures. The balance is transported into the matrix via the diffusive mass flow according to the method used in T2R3D. Once inside the matrix, the radionuclide can be advected through the matrix only at a much slower rate than through the fractures unless it transports into the fracture again at some later time either by advection or dispersion. This approach apparently yields a slower overall transport to the water table than the FEHM simulation, which adds additional residence time to particles in fracture elements rather than allowing them to actually transport into the matrix via the diffusive mass flow. The median breakthrough time of radionuclides with matrix diffusion is about 100 years for the FEHM simulation and several thousand years for the T2R3D simulation. Similar results are obtained in Figure 5, which shows the normalized breakthrough of neptunium at the water table for the FEHM and T2R3D simulations. These simulations include sorption in the vitric and zeolitic matrix elements. The runs with no matrix diffusion or dispersion show little difference in the initial breakthrough of neptunium except for the sharpness of the front, similar to the technetium simulations. For neptunium, however, the secondary breakthrough is delayed past 100,000 years because of sorption in the matrix. The runs with matrix diffusion show a disparity between the results of T2R3D and FEHM that are similar to the technetium runs. The reasons are the same for both radionuclides, but the difference is even more pronounced when sorption occurs in the matrix. The median breakthrough time for neptunium is nearly a thousand years for the FEHM simulation, but it is nearly 100,000 years for the T2R3D simulation. These results indicate that a significant difference exists in representations of the diffusive mass flow between fracture and matrix in FEHM and T2R3D. The diffusive mass flow between the fracture and matrix model in T2R3D allows the radionuclides to diffuse into the matrix, yielding much lower initial breakthrough via the fractures. FEHM results are based on an analytical solution that accounts for transient gradients in the matrix (though not valid for the finite matrix and the flow field here), but the absence of radionuclide transport into the matrix via diffusion is Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 41 March 2000 less consistent with the dual-permeability formulation used in the flow simulations. The input and output filenames associated with these runs are described in Attachment II. 6.4.4 Comparison of Numerical Solutions (T2R3D) with DCPT for Full 3-D Model of Yucca Mountain Site The full 3-D model of the Yucca Mountain unsaturated zone is a comprehensive, mountain-scale model. It includes all known aspects of flow and transport processes in the fractured-porous media, and provides a comprehensive test case for the particle-tracking simulator and other numerical simulators. Comparison of the particle tracker (DCPT) with the numerical solutions (T2R3D) provides insights into these methods for a complex system. A comparison of FEHM particle-tracker with DCPT for full 3-D model of Yucca Mountain Site can be found in AMR U0160 (CRWMS M&O 2000a, Section 6.2.4, pp. 21–22). That comparison shows discrepancies similar to those found in Section 6.4.3 of this report. Figure 6 shows a plan view of the 3-D grid. For these simulations, the radionuclides are released at the simulated repository horizon with time zero as a pulse. Steady-state flow is calculated using TOUGH2 V1.4 (TOUGH2 V1.4, STN: 10007-1.4-00, Version 1.4) with hydraulic properties in DTN: LB997141233129.001. The transport parameters are the same as those in Table 5. A total of 1,680 particles are used in the simulations using DCPT. The corresponding CPU time used for each run is about 20 seconds using DCPT on a Pentium II PC and about 1 hour using T2R3D on a DEC ALPHA. The cumulative mass fractions entering groundwater versus time are depicted on Figures 7 (technetium) and 8 (neptunium). The results of DCPT agree very well with the results of T2R3D. This argument implies that DCPT can provide results nearly identical to those of T2R3D, which rigorously solves the advection-dispersion equation of radionuclide transport in the Yucca Mountain site. Its performance will not diminish as the size of the grid (number of cells) increases, a feature that is particularly important in large-scale models such as the UZ Model of for Yucca Mountain. Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 42 March 2000 Based on data from DTN: LB990501233129.004 Figure 6. Map View of the 3-D Grid 167400 231000 169500 233100 171700 235300 173800 237400 175900 239500 EASTING (m) NORTHING (m) Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 43 March 2000 Based on data submitted with this AMR under DTN: LB990901233129.001 Figure 7. Comparison between DCPT and T2R3D for 3-D Radionuclide Transport of Technetium Based on data submitted with this AMR under DTN: LB990901233129.001 Figure 8. Comparison between DCPT and T2R3D for 3-D Radionuclide Transport of Neptunium 0 0.2 0.4 0.6 0.8 1 Relative Mass 1 10 10 2 10 3 10 4 10 5 10 6 TIME (years) DCPT T2R3D 0 0.2 0.4 0.6 0.8 1 Relative Mass TIME (years) DCPT T2R3D 1 10 10 2 10 3 10 4 10 5 10 6 Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 44 March 2000 INTENTIONALLY LEFT BLANK Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 45 March 2000 7. CONCLUSIONS Different methods for simulating radionuclide transport in unsaturated, fractured media were compared under conditions consistent with those expected at Yucca Mountain. These comparisons utilized 1-D and 3-D flow fields developed using the UZ Model, a dual-continua model calibrated to hydrologic conditions at Yucca Mountain. The methods compared included two particle-tracking methodologies, FEHM and DCPT, and one integral finite-difference method, T2R3D, which utilizes a fully coupled advective-dispersive solution. The latter method is considered to be a more rigorous approach, but is not always appropriate for large-scale problems because of its computational requirements. The modeling results reported in this AMR have been submitted to the TDMS under DTN: LB990901233129.001 and DTN: SN9908T0581699.001. The advantage of using a particle-tracking model (DCPT or FEHM) over a fully coupled advective-dispersive simulator (T2R3D) would be in its computational efficiency and lower CPU requirement, with less numerical diffusion in the case of small physical diffusion coefficients. The comparisons of T2R3D and DCPT revealed that DCPT provides results nearly identical to those of T2R3D for the time frames and scenarios considered. It can effectively simulate complex transport processes of radionuclides in dual-continua media. It is an efficient simulator, in terms of computational requirements, especially when only a cumulative breakthrough curve is required. Its performance will not diminish as the number of the grid cells increases, a feature that is of particular importance in large-scale models. Additionally, the DCPT provides higher spatial resolution since it allows particles to move through a continuous space. One-dimensional comparisons performed using the FEHM particle-tracking method and T2R3D indicated that the two methods agree only if diffusion and dispersion are neglected. For the cases that include diffusion and dispersion, the median breakthrough for FEHM occurred at times more than one to two orders of magnitude earlier than the simulations for T2R3D for the scenarios considered. This difference resulted from the use of a residence-time-transfer function to account for the effects of the diffusive mass flow between the fracture and the matrix in FEHM. Particles advected and dispersed in the fracture continuum are modeled as if they remain along these fast flow paths, and the residence-time-transfer-function algorithm is utilized to adjust particle residence times to reflect the time lag attributed to diffusion into and out of the matrix. This difference between T2R3D and FEHM is more pronounced for radionuclides undergoing sorption in the matrix. Numerical experiments reveal that the diffusive mass flow between fractures and the matrix is one of the key processes that control the travel time of radionuclides to water table in the Yucca Mountain, even though the dispersion processes in either fractures or the matrix have little effect. This notable difference in the results for the particle-tracking methods stems from different implementations of the diffusive mass flow between fractures and the matrix in the two codes. Essentially, as noted above, FEHM utilizes a residence-time-transfer function in accounting for diffusion into matrix, resulting in a formulation less consistent with the dual-permeability approach. As a result, the total mass flow from the fracture into the matrix is underestimated relative to a fully coupled advective-dispersive solution. With DCPT, both advection and dispersion/diffusion are incorporated simultaneously into the particle-transfer probability, providing an approach more consistent with the dual-permeability approach. As such, the DCPT Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 46 March 2000 is a better suited particle-tracking methodology than FEHM for a dual-continua model with a structure similar to that of the UZ Model. For a 10,000-year period, particle tracking using FEHM produces more conservative results by overpredicting the mass of radionuclides that will reach the water table. FEHM has already been used for transport simulations in the TSPA-VA, and past results should be considered conservative given the analysis presented here. Continued use of this code would not underestimate risk and, therefore, would not be invalid from a federal or state regulatory viewpoint. Its use, though, will underestimate the performance of the unsaturated zone as a barrier to radionuclide transport. Utilizing DCPT or T2R3D or similar approaches possibly implemented in FEHM for TSPA calculations would result in better calculated performance of the unsaturated zone, potentially by orders of magnitude compared with FEHM. Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 47 March 2000 8. INPUTS AND REFERENCES 8.1 DOCUMENTS CITED Bear, J. 1972. Dynamics of Fluid in Porous Media. New York, New York: Dover Publications. TIC: 217568. CRWMS M&O (Civilian Radioactive Waste Management System, Management & Operations Contractor) 1999a. Analysis and Modeling Report Development Plan (DP) for U0155 Analysis Comparing Advective-Dispersive Transport Solution to Particle-tracking, Rev. 00. TDP-NBSHS- 000002. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19990826.0106. CRWMS M&O 1999b. M&O Site Investigations. Activity Evaluation. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19990317.0330. CRWMS M&O 1999c. M&O Site Investigations. Activity Evaluation. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19990817.0257. CRWMS M&O 1999d. Performance Assessment Operations. Activity Evaluation. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19990716.0106. CRWMS M&O 2000a. Analysis of Base-Case Particle Tracking Results of the Base-Case Flow Fields. ANL-NBS-HS-000024. Las Vegas, Nevada: CRWMS M&O. CRWMS M&O 2000b. Particle Tracking Model and Abstraction of Transport Processes. ANLNBS- HS-000026. Las Vegas, Nevada: CRWMS M&O. URN-0037 Doughty, C. 1999. "Investigation of Conceptual and Numerical Approaches for Evaluating Moisture, Gas, Chemical, and Heat Transport in Fractured Unsaturated Rock." Journal of Contaminant Hydrology, 38 (1-3), 69-106. Amsterdam, The Netherlands: Elsevier Science Publishers. TIC: 244160. Dyer, J.R. 1999. "Revised Interim Guidance Pending Issuance of New U.S. Nuclear Regulatory Commission (NRC) Regulations (Revision 01, July 22, 1999), for Yucca Mountain, Nevada." Letter from J.R. Dyer (DOE) to D.R. Wilkins (CRWMS M&O), September 9, 1999, OL&RC:SB- 1714, with enclosure, "Interim Guidance Pending Issuance of New U.S. Nuclear Regulatory Commission (NRC) Regulations (Revision 01)." ACC: MOL.19990910.0079. LaBolle, E.M.; Fogg; G.E.; and Tompson, A.F.B. 1996. "Random-Walk Simulation of Transport in Heterogeneous Porous Media: Local Mass-Conservation Problem and Implementation Methods." Water Resources Research, 32 (4), 583-593. Washington, D.C.: American Geophysical Union. TIC: 245563. Pruess, K. 1991. TOUGH2 – A General Purpose Numerical Simulator for Multiphase Fluid and Heat Flow. Report LBL-29400. Berkeley, California: Lawrence Berkeley National Laboratory. ACC: NNA.19940202.0088. Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 48 March 2000 Sudicky, E.A. and Frind, E.O. 1982. "Contaminant Transport in Fractured Porous Media: Analytical Solutions for a System of Parallel Fractures." Water Resources Research, 18 (6), 1634-1642. Washington, D.C.: American Geophysical Union. TIC: 217475. Tang, D.H.; Frind, E.O; and Sudicky, E.A 1981. "Contaminant Transport in Fractured Porous Media: Analytical Solution for a Single Fracture." Water Resources Research, 17 (3), 555-564. Washington, D.C.: American Geophysical Union. TIC: 225358. Wu, Y.S.; Ahlers, C.F.; Fraser, P.; Simmons, A.; and Pruess, K. 1996. Software Qualification of Selected TOUGH2 Modules. Report LBNL-39490. Berkeley, California: Lawrence Berkeley National Laboratory. ACC: MOL.19970219.0104. Wu, Y. S. and K. Pruess. 1998. "A 3-D Hydrodynamic Dispersion Model for Modeling Tracer Transport in Geothermal Reservoirs." Proceedings, Twenty-third Workshop, Geothermal Reservoir Engineering, Stanford, California, January 26-28, 1998, 139-146. Stanford, California: Stanford University. TIC: 245292. Wu, Y.S.; Haukwa, C. and Bodvarsson, G.S. 1999. "A Site-Scale Model for Fluid and Heat Flow in the Unsaturated Zone of Yucca Mountain, Nevada." Journal of Contaminant Hydrology, 38 (1- 3), 185-215. Amsterdam, The Netherlands: Elsevier Science Publishers. TIC: 244160. SOFTWARE CITED Software Code: TOUGH2 V.1.4, STN: 10007-1.4-01. Software Code: T2R3D V.1.4, STN: 10006-1.4-00. Software Code: DCPT V.1.0, STN: 10078-1.0-00. Software Code: FEHM V.2.0, STN: 10031-2.00-00. Macro/Routine: MAKEPTRK V.1.0. ACC: MOL.19990915.0361. Macro/Routine: PROCESS1 V.1.0. ACC: MOL.19990915.0360. Macro/Routine: T2FEHM2 V.2.0. ACC: MOL.19990915.0359 Macro/Routine: PrepareKDfile V1.0. ACC: MOL.20000127.0120. Macro/Routine: ExtractFlow V1.0. ACC: MOL.20000127.0121. Macro/Routine: ExBT V1.0. ACC: MOL.20000127.0122. Macro/Routine: StatSpatial V1.0. ACC: MOL.20000202.0193 8.2 STANDARDS, CODES, REGULATIONS AND PROCEDURES CITED 64 FR (Federal Register) 8640. Disposal of High-Level Radioactive Waste in a Proposed Geologic Repository at Yucca Mountain. Proposed rule 10 CFR 63. Readily available. Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 49 March 2000 DOE (U.S. Department of Energy) 1998. Quality Assurance and Requirements Description. DOE/RW-0333P, REV 8. Washington D.C.: DOE OCRWM. ACC: MOL.19980601.0022. QAP-2-3 Classification of Permanent Items. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19990316.0006. 8.3 SOURCE DATA, LISTED BY DATA TRACKING NUMBER LB971212001254.001. DKM Basecase Parameter Set for UZ Model with Mean Fracture Alpha, Present Day Infiltration, and Estimated Welded, Non-Welded and Zeolitic FMX. Submittal date: 12/12/1997. LB990501233129.004. Mesh Files for 3-D UZ Model Calibration. Submittal date: 09/24/1999. LB997141233129.001. Calibrated Basecase Infiltration 1-D Parameter Set for The UZ Flow and Transport Model, FY99. Submittal date: 07/21/1999. 8.4 AMR OUTPUT DATA LISTED BY DATA TRACKING NUMBER LB99091233123.001. Modeling Results for DCPT Comparisons. Submittal date: 09/27/1999. SN9908T0581699.001. Files to Support 1-D Comparison Between FEHM Particle Tracking and T2R3D Advective-Dispersive Transport Simulations along SD-9. Submittal date: 08/16/1999. Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 50 March 2000 INTENTIONALLY LEFT BLANK Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 51 March 2000 9. ATTACHMENTS ATTACHMENT I - DOCUMENT INPUT REFERENCE SHEET ATTACHMENT II - INPUT AND OUTPUT FILES FOR DCPT AND FEHM 1. FILES FOR DCPT 2. FILES FOR FEHM ATTACHMENT III - SOFTWARE ROUTINES Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 52 March 2000 INTENTIONALLY LEFT BLANK Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 Attachment I-1 March 2000 ATTACHMENT I—DOCUMENT INPUT REFERENCE SHEET DIRS as of the issue date of this AMR. Refer to the DIRS database for the current status of these inputs. OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev ANL-NBS-HS-000001/00 Change: Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 2a 1. DTN: LB971212001254.001. DKM Basecase Parameter Set for UZ Model With Mean Fracture Alpha, Present Day Infiltration, and Estimated Welded, Non- Welded and Zeolitic FMX. Parameters & mesh file N/A – Reference only 6.4 Flow parameters for DKM Basecase, mean permeability, present day infiltration MESH file for present day infiltration N/A N/A N/A N/A 2. DTN: LB990501233129.004. Mesh Files for 3-D UZ Model Calibration. Entire N/ATechnical Product Output 4 3-D Grid Mesh-filename N/A N/A N/A N/A 3. DTN: LB997141233129.001. Calibrated Basecase Infiltration 1-D Parameter Set for the UZ Flow and Transport Model, FY99. 1dbasecaseR1 wodis.xls N/ATechnical Product Output 6.4 Flow parameters for base-case, present infiltration N/A N/A N/A N/A 4. Bear, J. 1972. Dynamics of Fluid in Porous Media. New York, New York: Dover Publications. TIC: 217568. p. 70 p. 633 N/A - Reference only 6.2.1 6.4.2 Equation 4.1.18 Equation 10.6.34 N/A N/A N/A N/A Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 Attachment I-2 March 2000 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev ANL-NBS-HS-000001/00 Change: Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 5. CRWMS M&O (Civilian Radioactive Waste Management System, Management & Operations Contractor) 1999a. Analysis and Modeling Report Development Plan (DP) for U0155 Analysis Comparing Advective- Dispersive Transport Solution to Particletracking, Rev. 00. TDPNBS- HS-000002. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19990826.0106. Entire N/A - Reference only 2 Planning Document N/A N/A N/A N/A 6. CRWMS M&O 1999b. M&O Site Investigations. Activity Evaluation. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19990317.0330. Entire N/A - Reference only 2 Activity Evaluation N/A N/A N/A N/A Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 Attachment I-3 March 2000 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev ANL-NBS-HS-000001/00 Change: Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 7. CRWMS M&O 1999c. M&O Site Investigations. Activity Evaluation. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19990817.0257. Entire N/A - Reference only 2 Activity Evaluation N/A N/A N/A N/A 8. CRWMS M&O 1999d. Performance Assessment Operations. Activity Evaluation. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19990716.0106. Entire N/A - Reference only 2 Activity Evaluation N/A N/A N/A N/A 9. CRWMS M&O 2000a. Analysis of Base-Case Particle Tracking Results of the Base-Case Flow Fields. ANL-NBSHS- 000024. Las Vegas, Nevada: CRWMS M&O. p. 21–22 N/A - Reference only 6.4 Model Comparison N/A N/A N/A N/A Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 Attachment I-4 March 2000 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev ANL-NBS-HS-000001/00 Change: Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 10. CRWMS M&O 2000b. Particle Tracking Model and Abstraction of Transport Processes. ANL-NBS-HS-000026. Las Vegas, Nevada: CRWMS M&O. URN-0037 6.2.4 N/A – Reference only 6.1 Description of Methodology N/A N/A N/A N/A 11. Doughty, C. 1999. “Investigation of Conceptual and Numerical Approaches for Evaluating Moisture, Gas, Chemical, and Heat Transport in Fractured Unsaturated Rock.” Journal Of Contaminant Hydrology, 38 (1–3), 69–106. Amsterdam, The Netherlands: Elsevier Science Publishers. TIC: 244160. pp. 100– 104 N/A – Reference only 5.1 Rationale for dual continua approach N/A N/A N/A N/A Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 Attachment I-5 March 2000 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev ANL-NBS-HS-000001/00 Change: Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 12. Dyer, J.R. 1999. “Revised Interim Guidance Pending Issuance of New U.S. Nuclear Regulatory Commission (NRC) Regulations (Revision 01, July 22, 1999), for Yucca Mountain, Nevada.” Letter from J.R. Dyer (DOE) to D.R. Wilkins (CRWMS M&O), September 9, 1999, OL&RC:SB-1714, with enclosure, “Interim Guidance Pending Issuance of New U.S. Nuclear Regulatory Commission (NRC) Regulations (Revision 01).” ACC: MOL.19990910.0079. Entire N/AReference only 4.2 Interim guidance N/A N/A N/A N/A Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 Attachment I-6 March 2000 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev ANL-NBS-HS-000001/00 Change: Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 13. LaBolle, E.M.; Fogg; G.E.; and Tompson, A. F.B. 1996. “Random- Walk Simulation of Transport in Heterogeneous Porous Media: Local Mass- Conservation Problem and Implementation Methods.” Water Resources Research, 32 (4), 583–593. Washington, D.C.: American Geophysical Union. TIC: 245563. pp. 583– 593 N/A - Reference only 6.2.2 Equations 3 and 10 N/A N/A N/A N/A 14. Pruess, K. 1991. TOUGH2-A General Purpose Numerical Simulator for Multiphase Fluid and Heat Flow. Report LBL-29400. Berkeley, California: Lawrence Berkeley National Laboratory. ACC: NNA.19940202.0088. Entire N/A - Reference only 6.3 General software use N/A N/A N/A N/A Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 Attachment I-7 March 2000 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev ANL-NBS-HS-000001/00 Change: Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 15. Sudicky, E.A. and Frind, E.O. 1982. “Contaminant Transport in Fractured Porous Media: Analytical Solutions for a System of Parallel Fractures.” Water Resources Research, 18 (6), 1634–1642. Washington, D.C.: American Geophysical Union. TIC: 217475. PP. 1634– 1642 N/A - Reference only 6.4 Used analytical solution N/A N/A N/A N/A 16. Tang, D.H.; Frind, E.O; and Sudicky, E.A 1981. “Contaminant Transport in Fractured Porous Media: Analytical Solution for a Single Fracture.” Water Resources Research, 17 (3), 555–564. Washington, D.C.: American Geophysical Union. TIC: 225358. pp. 555– 564 N/AReference only 6.4.3 Analytical solution for matrix diffusion model. N/A N/A N/A N/A Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 Attachment I-8 March 2000 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev ANL-NBS-HS-000001/00 Change: Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 17. Wu, Y.S.; Ahlers, C.F.; Fraser, P.; Simmons, A.; and Pruess, K. 1996. Software Qualification of Selected TOUGH2 Modules. Report LBNL- 39490. Berkeley, California: Lawrence Berkeley National Laboratory. ACC: MOL.19970219.0104. Entire N/A – Reference only 6.3 General software use N/A N/A N/A N/A 18. Wu, Y. S. and K. Pruess. 1998. “A 3-D Hydrodynamic Dispersion Model for Modeling Tracer Transport in Geothermal Reservoirs.” Proceedings, Twentythird Workshop, Geothermal Reservoir Engineering, Stanford, California, January 26– 28, 1998, 139–146. Stanford, California: Stanford University. TIC: 245292. pp.139–146 N/A - Reference only 6.3 Projected area weighting method N/A N/A N/A N/A Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 Attachment I-9 March 2000 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev ANL-NBS-HS-000001/00 Change: Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 19. Wu, Y.S.; Haukwa, C. and Bodvarsson, G.S. 1999. “A Site-Scale Model for Fluid and Heat Flow in the Unsaturated Zone of Yucca Mountain, Nevada.” Journal Of Contaminant Hydrology, 38 (1–3), 185–215. Amsterdam, The Netherlands: Elsevier Science Publishers. TIC: 244160. pp. 187– 188 190–193 N/A – Reference only 5.1 5.3 Rationale for dual continua approach Discussion of approximations in numerical methods N/A N/A N/A N/A 20. Software Code: TOUGH2 V.1.4, STN: 10007-1.4-01. Entire N/AQualified/ Verified/ Confirmed 6 General software use N/A N/A N/A N/A 21. Software Code: T2R3D V.1.4, STN: 10006-1.4- 00. Entire N/AQualified/ Verified/ Confirmed 6 General software use N/A N/A N/A N/A 22. Software Code: DCPT V.1.0, STN: 10078-1.0- 00 Entire TBV-3156 6 General software use 1 . N/A N/A Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 Attachment I-10 March 2000 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev ANL-NBS-HS-000001/00 Change: Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 23. Software Code: FEHM V.2.0, STN: 10031- 2.00-00 Entire N/AQualified/ Verified/ Confirmed 6 General software use N/A N/A N/A N/A 24. Macro/Routine: MAKEPTRK V.1.0. ACC: MOL.19990915.0361. Entire N/AQualified/ Verified/ Confirmed 3 Routine to create transport parameter file for FEHM N/A N/A N/A N/A 25. Macro/Routine: PROCESS V.1.0. ACC: MOL.19990915.0360. Entire N/AQualified/ Verified/ Confirmed 3 Routine to post-process results of FEHM particle-tracking simulations N/A N/A N/A N/A 26. Macro/Routine: T2FEHM2 V.2.0. ACC: MOL.19990915.0359 Entire N/AQualified/ Verified/ Confirmed 3 Routine to create FEHM-readable files from TOUGH2. N/A N/A N/A N/A 27. Macro/Routine: PrepareKDfile V1.0. ACC: MOL.20000127.0120. Entire N/AQualified/ Verified/ Confirmed 3 Routine to create a DCPTreadable file from a TOUGH mesh file and a T2R3D input file N/A N/A N/A N/A 28. Macro/Routine: ExtractFlow V1.0. ACC: MOL.20000127.0121. Entire N/AQualified/ Verified/ Confirmed 3 Routine to create a DCPTreadable file from the TOUGH output file. N/A N/A N/A N/A Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 Attachment I-11 March 2000 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev ANL-NBS-HS-000001/00 Change: Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 29. Macro/Routine: ExBT V1.0. ACC: MOL.20000127.0122. Entire N/AQualified/ Verified/ Confirmed 3 Routine to extract breakthrough curve from T2R3D output files. N/A N/A N/A N/A 30. Macro/Routine: StatSpatial V1.0. ACC: MOL.20000202.0193 Entire N/AQualified/ Verified/ Confirmed 3 Routine to calculate the distribution of particles along y=10.25 based on the output file of DCPT. N/A N/A N/A N/A AP-3.15Q.1 Rev. 06/30/1999 Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 II-1 March 2000 ATTACHMENT II—INPUT & OUTPUT FILES FOR DCPT AND FEHM 1. FILES FOR DCPT All the files listed here will be submitted with this AMR. The typical steps used in simulation with DCPT (the technetium case as an example) are shown below: Step 1: Prepare the input files and copy “UZ99.in” to “PTInput.txt” Step 2: Execute ParticleTrack.exe Step 3: Use standard spreadsheet software (Corel Quattro Pro 7.0), which is not subject to QARD, to calculate statistics of the exit time of particles contained in the file “UZ99_out.txt”, e.g., cumulative frequency scaled by the total number of particles. Table II-1. Files Involved in Section 6.4.1 Filename Description “FM1DR.in” Control file. List of all input/output files and parameters used by DCPT "FM1D_m.tec" List of the flow field and other transport parameters of the matrix for each cell "FM1D_f.tec" List of the flow field and other transport parameters of the fracture for each cell “FM1D.txt” Mesh file (cells, columns, and segments) “FM1Dini.txt” List of initial distribution of particles “FM1DOutR.txt” Output file, list of the final status of particles “Fm1D.wb3” A “Corel Quattro Pro 7” file which contains all postprocess results and comparisons with the analytical solutions U0155 Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking ANL-NBS-HS-000001 REV00 II-2 March 2000 Table II-2. Files Involved in Section 6.4.2 Filename Description “Analy3D.in” Control file. List of all input/output files and parameters used by DCPT “Analy3D_m.tec” List of the flow field and other transport parameters of the matrix for each cell “Analy3D_f.tec” List of the flow field and other transport parameters of the fracture for each cell “Analy3D.txt” Mesh file (cells, columns, and segments) “Ana3Dtextini.txt” List of initial distribution of particles “Analy3DOut.txt” Temporary output file, list of the final status of particles “ana3D2M.out” Output file, distribution of particles along the specific line in space (y=0) “Analy3d.wb3” A “Corel Quattro Pro 7” file which contains all postprocess results and comparisons with the analytical solutions Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 II-3 March 2000 Table II-3. Files Involved in Section 6.4.3 (Comparison of DCPT Part Only) Filename Description “UZ97_1D.in” Control file. List of all input/output files and parameters used by DCPT for the case without sorption (Technetium) “UZ97_1Dm.TEC” List of the flow field and other transport parameters of the matrix for each cell “Uz97_1df.tec” List of the flow field and other transport parameters of the fracture for each cell “UZ97_1DR.mesh” Mesh file (cells, columns, and segments) “UZ97_1DPT.ini” List of initial distribution of particles “UZ97_1Dout.txt” Temporary output file, list of the final status of particles. The results are loaded into the spreadsheet file before another run of DCPT “UZ97_1DFMD.dat” List of the characteristic distances of the fracture systems in each cell “UZ97_1D.flow” List water flow rates (via both fracture and matrix) per connections of neighboring cells (part of TOUGH2 output) “UZ97_1Dcon.dat” Configuration of cell connections in the grid “UZ97_1D.kd” List of values of Kd and bulk density of related cells “UZ97_1DR.in” Control file. List of all input/output files and parameters used by DCPT for the case with sorption (Neptunium) “UZ971D.wb3” A “Corel Quattro Pro 7” file which contains all postprocess results and comparisons with the numerical solutions for the case without sorption (Technetium) “Uz971DR.wb3” A “Corel Quattro Pro 7” file which contains all postprocess results and comparisons with the numerical solutions for the case with sorption (Neptunium) U0155 Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking ANL-NBS-HS-000001 REV00 II-4 March 2000 Table II-4. Files Involved in Section 6.4.4 Filename Description “UZ99.in” Control file. List of all input/output files and parameters used by DCPT for the case without sorption (Technetium) “UZ99_m.tec” List of the flow field and other transport parameters of the matrix for each cell “UZ99_f.tec” List of the flow field and other transport parameters of the fracture for each cell “UZ99mesh.txt” Mesh file (cells, columns, and segments) “UZ99PTini.txt” List of initial distribution of particles “UZ99_out.txt” Output file, list of the final status of particles “UZ99.flow” List water flow rates (via both fracture and matrix) per connections of neighboring cells (part of TOUGH2 output) “UZ99mesh.con” Configuration of cell connections in the grid “UZ99.kd” List of values of Kd and bulk density of related cells “UZ99DR.in” Control file. List of all input/output files and parameters used by DCPT for the case with sorption (Neptunium) “UZ99.wb3” A “Corel Quattro Pro 7” file which contains all postprocess results and comparisons with the numerical solutions for the case without sorption (Technetium) Title: Analysis Comparing Advective-Dispersive Transport Solution to Particle Tracking U0155 ANL-NBS-HS-000001 REV00 II-5 March 2000 2. FILES FOR FEHM The software and files used in this analysis have been submitted to the Technical Data Management System (TDMS) as part of the records submittal of this analysis (DTN: SN9908T0581699.001). A complete explanation of the files is contained in README files in each directory. For the runs specific to the FEHM particle-tracking comparison, the files are contained in the tar file AMR_U0155_Ho.tar. The tar file may be zipped and contains a .gz suffix. Any decompression software (e.g., WinZip) should be able to decompress the files and un-tar (extract) the subdirectories. In Unix, type "gunzip AMR_U0155_Ho.tar" to unzip the file. Then, type "tar xvf AMR_U0155_Ho" to extract the subdirectories. The following provides a description of the files and how they are used in the development and implementation of the FEHM particle tracking simulations. The 1-D TOUGH2 flow field is described by three files: sd9_e9.dt1, sd9_e9.ot1, and sd9_mesh. The rock properties and hydrologic properties are contained in sd9_e9.dt1 along with the infiltration source. The grid information is in sd9_mesh, and the output from the simulation is contained in sd9_e9.ot1. These files are used by T2FEHM2 to create FEHM-readable files that contain the same information. A complete description of the FEHM files created by T2FEHM2 can be found in Attachment III. The actual files are included and documented in the subdirectory 't2fehm2_files' in DTN: SN9908T0581699.001. T2FEHM2 is run only once since only one flow field is used in the com-parison study. The resulting files have the prefix 'fmsd9_e9.' An additional file not created by T2FEHM2 is required for the FEHM particle-tracking simulations. The 'ptrk' macro file contains transport parameter information for different materials in the model and is created by MAKEPTRK (see Attachment III). This pre-processor uses the sd9_e9.dt1 and sd9_mesh files as input. In addition, it requires user-specified information on transport properties such as sorption coefficients, diffusion coefficients, and dispersivities. ATTACHMENT III SOFTWARE ROUTINES Software Routine Name: MAKEPTRK Version: 1.0 Development Software: FORTRAN 77, Sun OS 5.7 Description: This is a software routine that creates the ‘ptrk’ macro in FEHM that describes the transport parameters for the particle tracking simulation. The ‘ptrk’ macro file contains transport parameter information for different materials in the model and is created by MAKEPTRK. This pre-processor uses the TOUGH2 ROCKS property file and mesh file as input to identify the different materials and the elements (nodes) that belong to those materials. In addition, it requires user-specified information on transport properties such as sorption coefficients, diffusion coefficients, and dispersivities. Below is a sample user-specified input file (“Np_diff.inp”) for the Neptunium particle tracking simulation with diffusion and dispersion (an explanation of each line entry and the actual input parameter name from the source file is given following the dashed line): sd9_e9.dt1 sd9_mesh Np_diff.ptrk 1 4 1. 4. 0. 0. 20. 1.6e-10 1. 1 -------------------------------------- write(*,*)'What is the name of the file containing the TOUGH2' write(*,*)'ROCKS card?' read(*,*) rocks write(*,*)'What is the name of the file containing the TOUGH2' write(*,*)'ELEME and CONNE cards?' read(*,*) mesh write(*,*)'What would you like to name the output file?' read(*,*) out write(*,*)'What transport mechanisms apply for the matrix?' write(*,*)'1 - advection only (no dispersion or matrix diff)' write(*,*)'2 - advection and dispersion (no matrix diff)' write(*,*)'3 - advection and matrix diff (no dispersion)' write(*,*)'4 - advection, dispersion, and matrix diff' read(*,*) iflagm write(*,*)'What transport mechanisms apply for the fracture?' write(*,*)'1 - advection only (no dispersion or matrix diff)' write(*,*)'2 - advection and dispersion (no matrix diff)' write(*,*)'3 - advection and matrix diff (no dispersion)' write(*,*)'4 - advection, dispersion, and matrix diff' read(*,*) iflagf write(*,*)'What is the Kd (cc/g) for vitric units?' read(*,*)xkdv write(*,*)'What is the Kd (cc/g) for zeolitic units?' read(*,*)xkdz write(*,*)'What is the Kd (cc/g) for devitrified units?' read(*,*)xkdd write(*,*)'What is the matrix dispersivity (m)?' read(*,*) dispm write(*,*)'What is the fracture dispersivity (m)?' read(*,*) dispf write(*,*)'What is the molecular diffusion coefficient?' read(*,*) do write(*,*)'What is the retardation factor for fracture' write(*,*)'sorption? (1 = no fracture sorption)' read(*,*) rdfrac write(*,*)'Would you like to use the f/m reduction factor in' write(*,*)'calculating aperture parameter? (1=yes,0=no)' read(*,*) nfm Because this routine simply reads the input parameters and places them into a formatted output file, there is no limitation as to the range of input parameters that is used. The parameters can be visually inspected to ensure that the input values have been correctly transferred to the output file (see verification below). The output from MAKEPTRK is a file that contains transport parameter information in a format that is required by the FEHM ‘ptrk’ macro. The information is pasted into a ‘master.ptrk’ file and renamed. A sample of a resulting ‘ptrk’ file for Neptunium with diffusion, dispersion, and sorption (‘fmNp_diff.ptrk’) that is used by FEHM is provided below: ptrk /* Np simulation with diffusion and dispersion 100000 204853 /* 100,000 particles, random # seed 204853 */ 0 1.e20 0. 1.e20 /* time for starting, ending trans. simulation,and time for ending, starting flow simultaion */ 1 0 2 2 /* print out particle information and store it in *.fin */ 1 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.160E-09 1.0 0.890E-01 0.100E-03 # 12 tswM4 1 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.160E-09 1.0 0.115E+00 0.100E-03 # 13 tswM5 1 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.160E-09 1.0 0.920E-01 0.100E-03 # 14 tswM6 1 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.160E-09 1.0 0.200E-01 0.100E-03 # 15 tswM7 1 0.100E+01 0.000E+00 0.000E+00 0.000E+00 0.160E-09 1.0 0.265E+00 0.100E-03 # 16 ch1Mv 1 0.100E+01 0.000E+00 0.000E+00 0.000E+00 0.160E-09 1.0 0.321E+00 0.100E-03 # 17 ch2Mv 1 0.100E+01 0.000E+00 0.000E+00 0.000E+00 0.160E-09 1.0 0.321E+00 0.100E-03 # 18 ch3Mv 1 0.100E+01 0.000E+00 0.000E+00 0.000E+00 0.160E-09 1.0 0.321E+00 0.100E-03 # 19 ch4Mv 1 0.400E+01 0.000E+00 0.000E+00 0.000E+00 0.160E-09 1.0 0.193E+00 0.100E-03 # 20 ch1Mz 1 0.400E+01 0.000E+00 0.000E+00 0.000E+00 0.160E-09 1.0 0.240E+00 0.100E-03 # 21 ch2Mz 1 0.400E+01 0.000E+00 0.000E+00 0.000E+00 0.160E-09 1.0 0.240E+00 0.100E-03 # 22 ch3Mz 1 0.400E+01 0.000E+00 0.000E+00 0.000E+00 0.160E-09 1.0 0.169E+00 0.100E-03 # 23 ch4Mz 1 0.100E+01 0.000E+00 0.000E+00 0.000E+00 0.160E-09 1.0 0.274E+00 0.100E-03 # 24 pp3Mv 1 0.400E+01 0.000E+00 0.000E+00 0.000E+00 0.160E-09 1.0 0.197E+00 0.100E-03 # 25 pp2Mz 1 0.100E+01 0.000E+00 0.000E+00 0.000E+00 0.160E-09 1.0 0.274E+00 0.100E-03 # 26 bf3Mv 1 0.400E+01 0.000E+00 0.000E+00 0.000E+00 0.160E-09 1.0 0.197E+00 0.100E-03 # 27 bf2Mz 1 0.100E+01 0.000E+00 0.000E+00 0.000E+00 0.160E-09 1.0 0.274E+00 0.100E-03 # 28 tm3Mv 4 0.000E+00 0.200E+02 0.200E+02 0.200E+02 0.160E-09 1.0 0.890E-01 0.412E-02 # 74 tswF4 4 0.000E+00 0.200E+02 0.200E+02 0.200E+02 0.160E-09 1.0 0.115E+00 0.114E-01 # 75 tswF5 4 0.000E+00 0.200E+02 0.200E+02 0.200E+02 0.160E-09 1.0 0.920E-01 0.119E-01 # 76 tswF6 4 0.000E+00 0.200E+02 0.200E+02 0.200E+02 0.160E-09 1.0 0.200E-01 0.103E-01 # 77 tswF7 4 0.000E+00 0.200E+02 0.200E+02 0.200E+02 0.160E-09 1.0 0.265E+00 0.930E-02 # 78 ch1Fv 4 0.000E+00 0.200E+02 0.200E+02 0.200E+02 0.160E-09 1.0 0.321E+00 0.100E-03 # 79 ch2Fv 4 0.000E+00 0.200E+02 0.200E+02 0.200E+02 0.160E-09 1.0 0.321E+00 0.100E-03 # 80 ch3Fv 4 0.000E+00 0.200E+02 0.200E+02 0.200E+02 0.160E-09 1.0 0.321E+00 0.100E-03 # 81 ch4Fv 4 0.000E+00 0.200E+02 0.200E+02 0.200E+02 0.160E-09 1.0 0.193E+00 0.821E-04 # 82 ch1Fz 4 0.000E+00 0.200E+02 0.200E+02 0.200E+02 0.160E-09 1.0 0.240E+00 0.821E-04 # 83 ch2Fz 4 0.000E+00 0.200E+02 0.200E+02 0.200E+02 0.160E-09 1.0 0.240E+00 0.821E-04 # 84 ch3Fz 4 0.000E+00 0.200E+02 0.200E+02 0.200E+02 0.160E-09 1.0 0.169E+00 0.821E-04 # 85 ch4Fz 4 0.000E+00 0.200E+02 0.200E+02 0.200E+02 0.160E-09 1.0 0.274E+00 0.100E-03 # 86 pp3Fv 4 0.000E+00 0.200E+02 0.200E+02 0.200E+02 0.160E-09 1.0 0.197E+00 0.100E-03 # 87 pp2Fz 4 0.000E+00 0.200E+02 0.200E+02 0.200E+02 0.160E-09 1.0 0.274E+00 0.100E-03 # 88 bf3Fv 4 0.000E+00 0.200E+02 0.200E+02 0.200E+02 0.160E-09 1.0 0.197E+00 0.100E-03 # 89 bf2Fz 4 0.000E+00 0.200E+02 0.200E+02 0.200E+02 0.160E-09 1.0 0.274E+00 0.100E-03 # 90 tm3Fv 1 1.00 0.000 0.00 0.00 0.160E-09 1.0 0.100E+00 0.100E-03 # 58 chaMd 35 4 0.00 20.00 20.00 20.00 0.160E-09 1.0 0.100E+00 0.164E-03 #115 chaFd 36 1 0 0 1 -12 0 0 1 -13 0 0 2 -14 0 0 3 -15 0 0 4 -16 0 0 5 -17 0 0 6 -18 0 0 7 -19 0 0 8 -20 0 0 9 -21 0 0 10 -22 0 0 11 -23 0 0 12 -24 0 0 13 -25 0 0 14 -26 0 0 15 -27 0 0 16 -28 0 0 17 -74 0 0 18 -75 0 0 19 -76 0 0 20 -77 0 0 21 -78 0 0 22 -79 0 0 23 -80 0 0 24 -81 0 0 25 -82 0 0 26 -83 0 0 27 -84 0 0 28 -85 0 0 29 -86 0 0 30 -87 0 0 31 -88 0 0 32 -89 0 0 33 -90 0 0 34 -58 0 0 35 -115 0 0 36 -500 0 0 -1. 0. 1.E-5 /* release particles at zone 500 */ The first four lines (note that the third line is wrapped) contain information for FEHM and are not relevant to the routine MAKEPTRK. The next 36 lines contain transport properties for different geologic layers of the system. These lines were extracted from the output of the MAKEPTRK file and only those materials at or beneath the repository were retained (geologic layers above the repository are not needed for simulations of radionuclide transport between the repository and the water table). The verification section below discusses the transport parameters in more detail. Following the blank line, the next 37 lines assign zones of nodes to each of the geologic layers. The final line is also irrelevant to MAKEPTRK and specifies the release of radionuclides. Verification: The sample output file shown above can be verified by visual inspection. A sample of a spot check is performed as follows. For material #74 (tswF4), the first column contains a flag that denotes the transport mechanism for this material. As identified in the input file, the transport mechanism for this fracture material should be denoted as “1” (advection, diffusion, and dispersion). The second column is the sorption coefficient, and it is correctly listed as “0” (for fractures). The next three columns are the dispersivity values for the x-, y-, and z-directions, and they are correctly listed as 20 m. The next column is the diffusion coefficient, which is correctly listed for Neptunium in the input file as 1.6E-10. The next column is the fracture sorption parameter, which is correctly listed as “1.” The next column is the corresponding matrix porosity (not actually used in this version of FEHM), which can be verified as correct by looking at the TOUGH2 ROCKS card (in ‘sd9_e9.dt1’ in DTN: SN9908T0581699.001) for material tswM4. The last number to be verified in these rows of transport properties is the fracture aperture parameter that is used to simulate matrix diffusion. The aperture parameter is calculated as the fracture element volume divided by the fracture/matrix connection area for that fracture element. The fracture/matrix connection area can be found in the CONNE card (in ‘sd9_mesh’ in DTN: SN9908T0581699.001) for connections between fracture and matrix elements. In MAKEPTRK, the fracture/matrix connection area can also be calculated as the product of the connection area supplied in the CONNE card and the reduction factor, Xfm, found in the ROCKS card to accommodate reductions in fracture/matrix conductance due to sub-grid heterogeneities. The latter calculation was used for this analysis, but it was learned after these calculations were performed that the formulation in FEHM does not need a reduction in fracture/matrix area to be consistent with the prescribed flow fields (the fracture saturation, which represents this reduction factor, is already accommodated in the FEHM formulation for matrix diffusion). Future revisions of this analysis should revise the aperture parameter calculation to exclude the reduction factor, but the general trends and results are not expected to change significantly. The fracture volume for an element (‘FlE71’) belonging to the material ‘tswF4’ is given in the ELEME card (in ‘sd9_mesh’ in DTN: SN9908T0581699.001) as 355.8 m3. The fracture/matrix connection area for this element is given in CONNE as 1.079E06 m2. The reduction factor is given in the ROCKS card for ‘tswF4’ as 0.008. The aperture parameter (as used in this analysis) is therefore equal to (35.58 m3) ÷ (1.079E06 m2) ÷ (0.008) = 4.12E-3 m. This is exactly the value reported in the sample output file. Note that the aperture parameter is only relevant for fracture elements, so the values for the matrix elements are “dummy” parameters. This verification ensures that MAKEPTRK is performing correctly for the range of input parameters that is used in this analysis. Listing of Software Routine MAKEPTRK v. 1.0: c makeptrk_v1.f c________________________________________________________________________ c This program will create the transport models that are used in the c FEHM ptrk macro. The required input files are the TOUGH2 ROCKS card, c ELEME card, and CONNE card. This program will also ask the user for c parameters including fracture and matrix diffusion, dispersivity, and c Kd. The primary output is, for each ROCKS material, the Kd, c dispersivity, molecular diffusion, fracture sorption, matrix porosity, c and aperture parameter (for fracture->matrix diffusion). c c C.K.Ho c 3/12/99 c c Several modifications have been made: c 1) Kd's are not assigned to fracture materials c 2) Format for dispersivity value has been changed from f5.2 to e10.3 c 3) User is given an option to use fracture/matrix reduction factor in c calculating aperture parameter. c C.K.Ho c 4/20/99 c________________________________________________________________________ c23456789012345678901234567890123456789012345678901234567890123456789012 c implicit double precision (a-h,o-z) character*22 block,rocks,mesh,out character*5 matname(999),mat(99999),elemn(99999),elem1,elem2 real por(999),xfm(999),vf(99999),afm(99999),bf(999) write(*,*)'What is the name of the file containing the TOUGH2' write(*,*)'ROCKS card?' read(*,*) rocks write(*,*)'What is the name of the file containing the TOUGH2' write(*,*)'ELEME and CONNE cards?' read(*,*) mesh write(*,*)'What would you like to name the output file?' read(*,*) out write(*,*)'What transport mechanisms apply for the matrix?' write(*,*)'1 - advection only (no dispersion or matrix diff)' write(*,*)'2 - advection and dispersion (no matrix diff)' write(*,*)'3 - advection and matrix diff (no dispersion)' write(*,*)'4 - advection, dispersion, and matrix diff' read(*,*) iflagm write(*,*)'What transport mechanisms apply for the fracture?' write(*,*)'1 - advection only (no dispersion or matrix diff)' write(*,*)'2 - advection and dispersion (no matrix diff)' write(*,*)'3 - advection and matrix diff (no dispersion)' write(*,*)'4 - advection, dispersion, and matrix diff' read(*,*) iflagf write(*,*)'What is the Kd (cc/g) for vitric units?' read(*,*)xkdv write(*,*)'What is the Kd (cc/g) for zeolitic units?' read(*,*)xkdz write(*,*)'What is the Kd (cc/g) for devitrified units?' read(*,*)xkdd write(*,*)'What is the matrix dispersivity (m)?' read(*,*) dispm write(*,*)'What is the fracture dispersivity (m)?' read(*,*) dispf write(*,*)'What is the molecular diffusion coefficient?' read(*,*) do write(*,*)'What is the retardation factor for fracture' write(*,*)'sorption? (1 = no fracture sorption)' read(*,*) rdfrac write(*,*)'Would you like to use the f/m reduction factor in' write(*,*)'calculating aperture parameter? (1=yes,0=no)' read(*,*) nfm open(1,file=mesh,status='old') open(3,file=rocks,status='old') open(12,file=out,status='new') c...Data c...Assign a dummy aperture parameter for matrix materials. c...Matrix diffusion is not used for matrix materials. bfm=1.e-4 c...Read in fracture information from MESH n=1 read(1,1000) block 1000 format(a22) 99 read(1,65) elemn(n),mat(n),vf(n) 65 format(a5,10x,a5,e10.4) c...End of active elements is signified by boundary elements or a c...blank space if(elemn(n).eq.' ') go to 98 if(elemn(n)(1:2).eq.'TP'.or. & elemn(n)(1:2).eq.'BT') go to 99 N=N+1 c...Read next line as matrix (assumes alternating listing) read(1,*) GO TO 99 98 CONTINUE NMAX = n - 1 c...Check do i=1,nmax write(*,*) vf(i),' ',mat(i) end do c...End check c...NMAX is the total number of fracture elements read from MESH write(*,107) nmax 107 format('Have read in ',i8,' fracture elements from MESH...') c...nnodes is the total number of active nodes c...Read in connection information from MESH N=1 c...Read header line CONNE READ(1,1500) BLOCK 1500 FORMAT(A22,3X,25X,E10.4) c...Read elements 1 and 2 and the connection area for F/M pairs only 199 read(1,1502) elem1,elem2,afm(n) 1502 format(2a5,40x,e10.4) IF(elem1(1:5).EQ.' '.OR.elem1(1:3).EQ.'+++') GO TO 198 if(elem1(1:1).ne.'F'.or.elem2(1:1).ne.'M') go to 199 N=N+1 GO TO 199 198 CONTINUE NCMAX = N - 1 c...NCMAX is the total number of f/m connections read from MESH write(*,203) ncmax 203 format('Have read in ',i8,' f/m connections from MESH...') c...Check do i=1,ncmax write(*,207) afm(i) 207 format(e10.4) end do c...End check c...Read in ROCKS information from TOUGH2 input file 18 read(3,1000) block if(block(1:5).ne.'ROCKS') go to 18 i=1 nfmat=0 nmmat=0 408 read(3,410) matname(i),drok,por(i) 410 format(a5,5x,2e10.4) if(matname(i).eq.'REFCO') go to 408 if(matname(i).eq.' ') then c...ntotmat is the total number of materials in the ROCKS card ntotmat=i-1 go to 27 end if read(3,*) c...nfmat is the total number of fracture materials if(matname(i)(3:3).eq.'F'.or.matname(i)(4:4).eq.'F') then nfmat=nfmat+1 end if c23456789012345678901234567890123456789012345678901234567890123456789012 c...nmmat is the total number of matrix materials if(matname(i)(3:3).eq.'M'.or.matname(i)(4:4).eq.'M') nmmat=nmmat+1 read(3,33) xfm(i) if(xfm(i).eq.0.) xfm(i)=1. 33 format(60x,e10.4) read(3,*) i=i+1 go to 408 27 continue write(*,75) nmmat 75 format('Number of matrix materials in ROCKS = ',i5) write(*,77) nfmat 77 format('Number of fracture materials in ROCKS = ',i5) c...Check do i=1,ntotmat write(*,*) xfm(i) end do c...End check c...Determine matrix porosities corresponding to each fracture material c...Because the number of fracture and matrix materials are not equal, c...I am comparing the characters of the element names. I first c...determine where the 'F' is, and then I compare all other c...characters with the matrix material to get a match. write(*,*) ntotmat do i=1,ntotmat if(matname(i)(3:3).eq.'M'.or.matname(i)(4:4).eq.'M')goto83 if(matname(i).eq.'topbd'.or.matname(i).eq.'botbd')goto83 do j=1,ntotmat if(matname(i)(3:3).eq.'F') then if(matname(j)(3:3).eq.'M') then if(matname(j)(1:2).eq.matname(i)(1:2).and. & matname(j)(4:5).eq.matname(i)(4:5)) then por(i)=por(j) go to 83 end if end if elseif(matname(i)(4:4).eq.'F') then if(matname(j)(4:4).eq.'M') then if(matname(j)(1:3).eq.matname(i)(1:3).and. & matname(j)(5:5).eq.matname(i)(5:5)) then por(i)=por(j) go to 83 end if end if end if end do por(i)=0.1 c23456789012345678901234567890123456789012345678901234567890123456789012 write(*,113) matname(i) 113 format('Material ',a5,' does not have a matrix counterpart.'/ & 'It has been assigned a matrix porosity of 0.1') 83 end do c...Determine aperture parameter, bf do i=1,ntotmat c...If material is boundary then assign a dummy aperture parameter if(matname(i).eq.'topbd'.or.matname(i).eq.'botbd') then bf(i)=bfm goto87 end if if(matname(i)(3:3).eq.'F'.or.matname(i)(4:4).eq.'F') then do j=1,nmax if(mat(j).eq.matname(i).and.nfm.eq.1) then bf(i)=vf(j)/(xfm(i)*afm(j)) c...Check write(*,*) bf(i) go to 87 elseif(mat(j).eq.matname(i).and.nfm.eq.0) then bf(i)=vf(j)/afm(j) c...End check go to 87 end if end do c...If a material cannot be associated with an active fracture element, c...then assign the material a dummy aperture parameter. bf(i)=bfm end if 87 end do c...Write data to output file for PTRK macro do i=1,ntotmat c...Assign appropriate Kd if(matname(i)(5:5).eq.'v') then xkd=xkdv elseif(matname(i)(5:5).eq.'z') then xkd=xkdz else xkd=xkdd end if c23456789012345678901234567890123456789012345678901234567890123456789012 if(matname(i)(3:3).eq.'M'.or.matname(i)(4:4).eq.'M') then write(12,505)iflagm,xkd,dispm,dispm,dispm,do,1.,por(i), & bfm,i,matname(i) 505 format(i1,1x,4(e10.3,1x),e9.3,1x,f4.1,1x,2(e10.3,1x),'#', & i3,1x,a5) else write(12,505)iflagf,0.,dispf,dispf,dispf,do,rdfrac,por(i), & bf(i),i,matname(i) end if end do write(12,*) do i=1,ntotmat write(12,507) -i,0,0,i 507 format(i5,1x,i1,1x,i1,1x,i5) end do stop end Software Routine Name: PROCESS1 Version: 1.0 Development Software: Fujitsu FORTRAN 90 Description: This is a software routine that post-processes the results of the FEHM particle tracking to provide columns of time, mass flux (mol/year) and cumulative mass at the water table (mol). The post-processor PROCESS1 is executed with an input file “process.dat” that is modified to reflect the desired output name of the run. This processor takes the information from the particle tracking code and prints the information to an output file named by the user. A sample input file (‘process.dat’) for PROCESS1 is given below: ../fmsd9_e9.grid ../fmsd9_e9.fin fmNp_nodiff.output 0.5 100 100 1.000 4 0.01 0.05 0.1 0.5 The first two lines are the names of input files, and the third line is the desired output file name. The fourth line contains information about how the post-processor bins the particles for printing the time (years), mass flux (mol/years), and cumulative breakthrough (mol) at the water table. The fifth line indicates how many numbers are in the sixth line, and the sixth line contains values for the percent cumulative breakthrough at which times are desired to be printed to the screen. The output file contains three columns. The first column is the time in years. The second column is the mass flow (mol/year) recorded at the water table. The third column is the cumulative mass (moles) that has reached the water table at the specified time. A sample of the output file ‘fmNp_nodiff.output’ is extracted below: The results of the PROCESS1 can be plotted directly. 1.16189623 9.70895290e-02 9.99999975e-06 1.17165995 0.902203918 4.72000008e-03 1.17587113 1.42056239 9.42999963e-03 1.17873192 1.98488736 1.41399996e-02 1.18089843 2.28967118 1.88500006e-02 1.18277979 2.60619116 2.35600006e-02 1.18448448 3.12997580 2.82700006e-02 1.18586791 3.37901926 3.29799987e-02 1.18721294 3.84510279 3.76899987e-02 1.18845415 3.89304090 4.23999988e-02 1.18956292 4.66002941 4.71099988e-02 1.19055974 4.46301603 5.18199988e-02 1.19154954 5.20216608 5.65299988e-02 1.19247949 5.29256201 6.12399988e-02 . . . 366005.031 0.142051578 0.900529981 366005.062 1.48951869e+10 0.905250013 366005.062 0.142051578 0.909969985 366005.062 1.48951869e+10 0.914690018 366005.062 0.142051578 0.919409990 414162.375 3.62423052e-08 0.924130023 618860.563 3.03773398e-08 0.928849995 651621.313 1.48951869e+10 0.933570027 651621.313 7.10257888e-02 0.938290000 Verification: This section contains a verification test of the software routine PROCESS1. The test consists of one of the 1-D simulations used in the FEHM particle tracking analysis (Tc with diffusion 3.2e- 11 m2/s and dispersion=20 m in fractures. Only 20 particles are used in the test case so that the output in the ‘fmsd9_e9.fin’ file can be directly processed and compared to the results of PROCESS1. The last row of numbers in the ‘*.fin’ file contains the times at which each of the 20 particles left the system (exited at the water table). These times are plotted as a cumulative distribution function. The output from PROCESS1 is in ‘fmTc_diff_test.output’ (see /fehmruns/process/test_process in DTN: SN9908T0581699.001). The first column is the time in years. The second column is the interpolated mass flux (mol/yr), and the third column is the cumulative breakthrough (mol). Because only one mole was injected, the third column is the same as the cumulative percent breakthrough given in the CDF plotted directly from the ‘*.fin’ file. The following plots are reproduced from ‘process1.doc’ (see /fehmruns/process/test_process in DTN: SN9908T0581699.001), which show that the post-processor is producing results that are the same as the actual values in the output file. The post-processor provides interpolation, so that curve is smoothed in some regions. This verification indicates that PROCESS1 is performing correctly for the range of input parameters used in this analysis. 1 1 0 1 00 1 000 1 04 1 05 0 2 0 4 0 6 0 8 0 1 00 fm sd9_e9.fin_tes t_CDF.qpc T ime (ye a rs) C um u la tive Pe rcen t B re a kthrou gh o f T c Cumulative Percent Breakthrough Using CDF of Values in 'fmsd9_e9.fin' 20 particles of T c with D iffusion (3.2x10-11 m 2/s ) 0 2 0 4 0 6 0 8 0 1 00 1 1 0 1 00 1 000 1 0 4 1 0 5 process1_test.qpc C um u lative P e rcen t B re a kthro ugh o f T c T im e (ye a rs ) Cumulative Percent B reakthrough Using 'process1.f' 20 particles of T c w ith D iffusion (3.2x10-11 m 2/s ) Listing of Software Routine PROCESS1 v. 1.0: program process1 implicit none character(100) dummy_string, grid_file, fin_file, out_file character(4) gas_flag, ptrk_flag, dpdp_flag, dual_flag integer, allocatable :: ifinal(:), index(:) integer i, j, n0, npart, nseed, neq, ic, npartbin, nbins1, npbin integer nbins2, npstart, npartfin, nfraction_out, len_aaxy real(4) total_mass, flux, delta_time, sumtime, fraction1, fluxmax real(4) timemax real(8), allocatable :: rdum(:), a_axy(:) real(4), allocatable :: timep(:) real(4), allocatable :: fraction_out(:) integer, allocatable :: ifraction_out(:) c process.dat contains files names, histogram parameters open(1,file='process.dat') c Read name of grid file, .fin file, then open them read(1,'(a100)') grid_file read(1,'(a100)') fin_file read(1,'(a100)') out_file open(3,file= grid_file) open(4,file = fin_file) c Read number of nodes from grid file, then close read(3,'(a100)') dummy_string read(3,*) neq close(3) c open output file open(7,file= out_file) c Read 3 dummy lines, then get gas flag, ptrk flag read(4,'(a100)') dummy_string read(4,'(a100)') dummy_string read(4,'(a100)') dummy_string read(4,'(a4)') gas_flag read(4,'(a4)') ptrk_flag c Read dual and dpdp flags to tell if either option was used read(4,'(a100)') dummy_string read(4,'(a4)') dpdp_flag read(4,'(a4)') dual_flag c Set n0 based on ECM, DPDP, or DUAL if(dpdp_flag .eq. 'dpdp') then n0 = 2*neq elseif(dual_flag .eq. 'dual') then n0 = 3*neq else n0 = neq end if c Allocate space for the array to read state variables allocate(rdum(n0)) c read in state variables based on what type of gas option was used if(gas_flag .eq. 'ngas') then read(4,'(4g20.10)')(rdum(i),i=1,n0) read(4,'(4g20.10)')(rdum(i),i=1,n0) read(4,'(4g20.10)')(rdum(i),i=1,n0) read(4,'(4g20.10)')(rdum(i),i=1,n0) elseif(gas_flag(1:3) .eq. 'air') then read(4,'(4g20.10)')(rdum(i),i=1,n0) read(4,'(4g20.10)')(rdum(i),i=1,n0) else read(4,'(4g20.10)')(rdum(i),i=1,n0) read(4,'(4g20.10)')(rdum(i),i=1,n0) read(4,'(4g20.10)')(rdum(i),i=1,n0) end if c rdum no longer needed deallocate(rdum) c Determine if mass fluxes are in file, read if they are read(4,'(a100)') dummy_string if(dummy_string(1:4).eq.'mass') then read(4,*) len_aaxy allocate(a_axy(len_aaxy)) read(4,'(5g15.8)') (a_axy(i),i=1,len_aaxy) deallocate(a_axy) else backspace 4 end if c Read in number of particles, seed value read(4,*) npart, nseed c Allocate space for final node array and time array allocate(ifinal(npart)) allocate(index(npart)) allocate(timep(npart)) read(4,*)(ifinal(i),i=1,npart) c Skip through two other output arrays to get to the time array c if the user wrote these arrays out if(nseed .gt. 0) then read(4,*)(timep(i),i=1,npart) read(4,*)(timep(i),i=1,npart) end if c Read array of leaving times (or particle age if ifinal>0) read(4,*)(timep(i),i=1,npart) c Loop takes all particles that have left the system and shifts them c to the first positions in the arrays so that the sorting routine c will not have to deal with particles that are still in the system ic = 0 do i = 1, npart if(ifinal(i) .lt. 0) then ic = ic + 1 ifinal(ic) = ifinal(i) timep(ic) = timep(i) end if end do c The number of particles to bin are now only the ones that left c the system npartbin = ic c read in number of bins, total mass of radionuclides c fraction1 - the first set of bins applies to the first c fraction1*npartbin particles c nbins1 - number of bins in which to bin the first set of particles c nbins2 - number of bins for the remaining particles c read(1,*) fraction1, nbins1, nbins2, total_mass read(1,*) nfraction_out allocate(fraction_out(nfraction_out)) allocate(ifraction_out(nfraction_out)) read(1,*)(fraction_out(i),i=1,nfraction_out) do i = 1, nfraction_out ifraction_out(i)=fraction_out(i)*npart end do c Call routine to sort the particles and node array c lowest to highest c This routine is an indexing and ranking sort routine that returns c the index array, such that timep(index(i)), i = 1, npartbin c are in ascending order. It doesn't sort the timep array itself, c but supplies the index array so that the order can be obtained c by indirect indexing call sort_parti(npartbin, npart, index, timep) c Set max flux to small number initially fluxmax = 0. c Do average time and mass flux calculation for each bin of first c partition npartfin = fraction1*npartbin npbin = npartfin/nbins1 bin_loop1: do i = 1, npartfin, npbin sumtime = 0. if(i+npbin-1.gt.npartfin) exit bin_loop1 do j = i, i+npbin-1 sumtime = sumtime + timep(index(j)) end do sumtime = sumtime / npbin delta_time = timep(index(i+npbin-1))-timep(index(i)) if(delta_time.eq.0.) delta_time = 1.e-5 flux = npbin * total_mass / (npart*delta_time) if(flux.gt.fluxmax) then fluxmax = flux timemax = sumtime/31557600. end if write(7,*) sumtime/31557600., 31557600.*flux, 2 real(i)/real(npart) end do bin_loop1 c Do average time and mass flux calculation for each bin of second c partition npstart = i npartfin = npartbin npbin = (npartbin-i+1)/nbins2 bin_loop2: do i = npstart, npartfin, npbin sumtime = 0. if(i+npbin-1.gt.npartfin) exit bin_loop2 do j = i, i+npbin-1 sumtime = sumtime + timep(index(j)) end do sumtime = sumtime / npbin delta_time = timep(index(i+npbin-1))-timep(index(i)) if(delta_time.eq.0.) delta_time = 1.e-5 flux = npbin * total_mass / (npart*delta_time) if(flux.gt.fluxmax) then fluxmax = flux timemax = sumtime/31557600. end if write(7,*) sumtime/31557600., 31557600.*flux, 2 real(i)/real(npart) end do bin_loop2 write(6,*)(timep(index(ifraction_out(i)))/31557600., 2 i=1,nfraction_out), timemax, 31557600.*fluxmax end subroutine sort_parti(n, nsize, indx, time) c c Indexing and Ranking algorithm for sorting, from Numerical Recipes c Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T. c Vetterling, 1986, Numerical Recipes. The Art of Scientific c Computing, Cambridge University Press, Cambridge, pp. 232-234. c implicit none real(4) time(nsize), q integer indx(nsize) integer nsize integer n, j, l, ir, indxt, i do j = 1, n indx(j) = j end do l = n/2+1 ir=n 10 continue if(l.gt.1) then l=l-1 indxt=indx(l) q=time(indxt) else indxt=indx(ir) q=time(indxt) indx(ir)=indx(1) ir=ir-1 if(ir.eq.1) then indx(1)=indxt return end if end if i=l j=l+l 20 if(j.le.ir) then if(j.lt.ir) then if(time(indx(j)).lt.time(indx(j+1))) j=j+1 end if if(q.lt.time(indx(j))) then indx(i)=indx(j) i=j j=j+j else j=ir+1 end if goto 20 end if indx(i)=indxt goto 10 end Software Routine Name: T2FEHM2 Version: 2.0 Development Software: FORTRAN 77, Sun OS 5.7 Description: The software routine T2FEHM2 was written to reformat TOUGH2 files that contain information pertaining to unsaturated flow to FEHM-readable files that can be used for radionuclide particle tracking. This method maintains consistency with the hydrologic conditions (mass flow rates, liquid saturations, etc.) prescribed in the TOUGH2 flow fields. FEHM uses a cell-based particle tracking model that preserves the overall residence time through any portion of the model and probabilistically reproduces the migration of a solute through the domain. The requirement for the method is that the flow calculation be based on a control volume in which fluid flow rates into and out of each cell are computed. Since TOUGH2 is an integrated finite difference code, and FEHM employs a control volume finite element technique, the two codes are compatible for implementing the particle tracking technique. The required inputs for FEHM to use an externally-developed flow field are: (1) grid connectivity information and cell volumes; (2) properties and state variables (rock grain density, fluid saturation, and rock porosity at each grid point); (3) inter-nodal fluid mass flow rate for every connection in the numerical grid; and (4) fluid source and sink flow rates for each grid block. The post-processor, T2FEHM2, was written to generate these required data from existing TOUGH2 files. The remainder of this section describes the required inputs to T2FEHM2 and the corresponding output files. Required Input Files for T2FEHM2 When executed, T2FEHM2 will prompt the user for the names of three required files: (1) TOUGH2 input file; (2) TOUGH2 output file; and (3) TOUGH2 mesh file. T2FEHM2 will also prompt the user for the name of a fourth file containing the names of repository elements, but this file is optional. All input files can be found in the subdirectory ‘fromLBNL’ in DTN: SN9908T0581699.001. TOUGH2 Input File The TOUGH2 input file must contain the ROCKS and GENER cards. ROCKS contains material property information for fracture and matrix materials corresponding to a dual-permeability model. Fracture and matrix materials must have an ‘F’ or ‘M’, respectively, in the third or fourth character of the material name. Each material must have four lines associated with its entry. The GENER card should contain information on the infiltration source terms for prescribed elements. The generation rate is specified in units of kg/s. TOUGH2 Output File The TOUGH2 output file contains all simulated state variables (pressure, saturation) for each element and flux variables (mass flow rate) for each connection pair at user-specified print-out times. T2FEHM2 reads in these state and flux variables and puts them in a format that is compatible with FEHM. TOUGH2 Mesh File The TOUGH2 mesh file contains the ELEME and CONNE cards. ELEME contains the element names, material names, volumes, and coordinates of each element in the TOUGH2 model. The fracture and matrix elements should be listed alternately with a fracture element listed first. Also, all boundary elements must be listed at the end of the ELEME card. The material names associated with each element should be five-character names (not integers) that correspond identically to the name of one of the materials in the ROCKS card. The CONNE card contains all connection pairs and associated connection information for each element in the TOUGH2 model. T2FEHM2 stores these connection pairs to create connectivity arrays (ncon, istrw, nelmdg) for FEHM. File Containing Repository Elements A file containing the names of repository elements is optional. If present, T2FEHM2 will read the number of repository elements in the first line of the file. All repository element names will then be read from the file. These elements will be used to create special fracture and matrix zones in a FEHM file that will be used to define the location of radionuclide release for particle tracking. For the 1-D FEHM simulation used in this analysis, only one element is specified for the repository zone. 6.1.2 Output Files from T2FEHM2 After reading the required information from the input files, T2FEHM2 prints out nine (9) files that are used by FEHM. The user specifies a reference file name, and the code creates nine output files by appending the following nine suffixes to the reference file name: .dat .dpdp .files .grid .ini .rock .stor .zone .zone2 All T2FEHM2 output files can be found in the subdirectory ‘t2fehm2_files’ in DTN: SN9908T0581699.001. A tenth file, ‘file_name.check,’ is also printed but it is not used by FEHM. This file contains the node numbers and number of connections for each node. More detailed information on the contents of the FEHM macros can be found in Zyvoloski et al. (1997). The user should consult this information because a number of these macros have been created with T2FEHM2 using “dummy” variables that are either not needed by the particle tracking simulation (e.g., permeability, area coefficients, element specifications for nodes, etc.) or that can be modified by the user to suit the specific needs of the particle tracking simulation (e.g., date, time steps, print-out options, etc.). Most of these prescribed variables appear in the ‘*.dat’ file, so the user should become familiar with the macros listed in that file before using the default values prescribed in T2FEHM2. The prefix “fm” is placed in front of all T2FEHM2 files for identification purposes. The remainder of this section details the specific output files. To verify that T2FEHM2 is producing correct results, portions of the actual output files used in this analysis (‘fmsd9_e9*’) are included. The values are compared to those in the original TOUGH2 files by visual inspection to ensure correct results for the range of inputs used in this analysis. Output File ‘*.dat’ This file contains the required macros used by FEHM: ‘dpdp,’ ‘perm,’ ‘rlp,’ ‘rock,’ ‘flow,’ ‘time,’ ‘ctrl,’ ‘iter,’ ‘sol,’ ‘rflo,’ ‘air,’ ‘node,’ ‘zone,’ ‘ptrk.’ If the macros are not explicitly defined in this file, the names of macro files containing the actual information are listed here. Macros ‘perm’ and ‘rlp’ are not required by the particle tracking solution, so dummy values are inserted here. In addition, many of the values in the ‘*.dat’ file are prescribed within T2FEHM2 as default values, so the user should refer to Zyvoloski et al. (1997) to modify the values in the different macros to suit their needs. The output file ‘fmsd9_e9.dat’ is provided below : "fmsd9_e9.dat" 47 lines, 770 characters # input file for mean alpha, fitted fmx, present day q, ysw # AR 11/19/97 # Particle tracking for TOUGH2 flow field dpdp file fmsd9_e9.dpdp perm 1 0 0 0.100E-14 0.100E-14 0.100E-14 rlp 1 0. 0. 1. 1. 0. 1. 1 0 0 1 rock file fmsd9_e9.rock flow time 0.36525E+09 0.36525E+09 10 10 1997 10 ctrl -10 0.10E-03 40 1 0 0 1 0 1.00 3.00 1.00 5 0.20E+01 0.10E-09 0.10E+11 0 1 iter 0.10E-04 0.10E-04 0.10E-04 -0.10E-03 0.12E+01 0 0 0 0 0.14E+05 sol 1 -1 rflo air -1 20.0 0.1 node 1 1 zone file fmsd9_e9.zone2 ptrk file fmsd9_e9.ptrk stop Output File ‘*.dpdp’ This file contains a list of the zones corresponding to the fracture materials and lists the fracture porosities. It also contains dummy information regarding the length scale for matrix nodes that is not required for the TOUGH2-FEHM coupling. Here are the first few lines from ‘fmsd9_e9.dpdp’ that can be compared to the ROCKS card used by TOUGH2: dpdp 1 -63 0 0 0.2330E-03 -64 0 0 0.2990E-03 -65 0 0 0.7050E-04 -66 0 0 0.4840E-04 -67 0 0 0.4830E-04 -68 0 0 0.1300E-03 -69 0 0 0.6940E-04 -70 0 0 0.3860E-04 -71 0 0 0.8920E-04 -72 0 0 0.1290E-03 -73 0 0 0.1050E-03 -74 0 0 0.1240E-03 -75 0 0 0.3290E-03 -76 0 0 0.3990E-03 Output File ‘*.files’ This control file contains a list of files that FEHM reads for necessary information. Below is the ‘fmsd9_e9.files’ file: fmsd9_e9.dat fmsd9_e9.grid fmsd9_e9.zone fmsd9_e9.out fmsd9_e9.ini fmsd9_e9.fin fmsd9_e9.his fmsd9_e9.trc fmsd9_e9.con fmsd9_e9.stor fmsd9_e9.chk all 0 Output File ‘*.grid’ This file contains the ‘coor’ and ‘elem’ macros. The first line of the ‘coor’ macro gives the total number of fracture elements, followed by a list of all the nodes in the fracture domain and their respective x, y, and z coordinates. The ‘elem’ macro contains dummy information regarding the nodes associated with each element, but this is not required for the TOUGH2-FEHM coupling. Below are the first few lines of the ‘fmsd9_e9.grid’ file that can be compared to the values in ELEME used by TOUGH2: coor 25 1 171270.58 234054.36 1289.80 2 171270.58 234054.36 1285.89 3 171270.58 234054.36 1281.98 4 171270.58 234054.36 1275.11 5 171270.58 234054.36 1263.82 6 171270.58 234054.36 1255.06 7 171270.58 234054.36 1242.07 8 171270.58 234054.36 1224.77 9 171270.58 234054.36 1214.57 Output File ‘*.ini’ This file contains re-start information for FEHM. The liquid saturations of all fracture and matrix nodes are listed following eight header lines. The gas-phase pressures (MPa) are then listed for the fracture and matrix nodes. The fourth header line (‘air’) tells FEHM that the pressures are for the gas phase. Then, mass flux values (kg/s) are listed for each connection of each node, starting with node 1 (the ordering is the same as the ‘ncon’ array in ‘.stor’ without pointer information—see ‘*.stor’ below). The mass flux values include sources (infiltration) denoted as negative values and sinks (connection to water table) denoted as positive values for each node. Flow into a node is negative, and flow out of a node is positive. The mass flux values for the fracture domain are listed first followed by the mass flux values in the matrix domain. The mass flux between fracture and matrix elements are listed last. Flow from the fracture to the matrix is denoted as positive. The file ‘fmsd9_e9.ini’ is shown below: "fmsd9_e9.ini" 71 lines, 4506 characters # input file for mean alpha, fitted fmx, present day q, ysw # AR 11/19/97 This is a .ini file with saturations, pressures and mass flux values. 0. air ptrk nstr dpdp ndua 0.70676000E-01 0.76186000E-01 0.10618000 0.66245000E-01 0.61774000E-01 0.37347000E-01 0.23308000E-01 0.22572000 0.13596000 0.78086000E-01 0.69016000E-01 0.89885000E-01 0.10026000 0.10023000 0.10020000 0.10021000 0.12480000 0.12220000 0.17218000 0.26035000 0.30390000 0.26406000 0.78540000 0.32361000 0.28965000 0.95558000 0.95624000 0.99071000 0.61900000 0.57734000 0.54192000 0.43204000 0.41860000 0.76104000 0.56915000 0.80942000 0.91995000 0.85113000 0.85505000 0.85825000 0.85443000 0.96344000 0.59510000 0.74785000 0.98919000 0.99920000 0.99482000 0.94023000 0.96382000 0.99339000 0.91999800E-01 0.92000100E-01 0.92000000E-01 0.92000000E-01 0.92000000E-01 0.92000000E-01 0.92002000E-01 0.92000000E-01 0.92000000E-01 0.92000000E-01 0.92000000E-01 0.92000000E-01 0.91999600E-01 0.91999800E-01 0.92000100E-01 0.91999800E-01 0.91999500E-01 0.92000400E-01 0.91999600E-01 0.91999700E-01 0.92000300E-01 0.92000000E-01 0.92000000E-01 0.92000400E-01 0.91999900E-01 0.92000000E-01 0.92000000E-01 0.92000000E-01 0.91995000E-01 0.92000000E-01 0.92004000E-01 0.92004000E-01 0.92000000E-01 0.92000000E-01 0.92000000E-01 0.92000000E-01 0.92000000E-01 0.92000000E-01 0.92000000E-01 0.92000000E-01 0.92000000E-01 0.92000000E-01 0.91998000E-01 0.92000000E-01 0.92000000E-01 0.91999800E-01 0.92000000E-01 0.92000000E-01 0.91998000E-01 0.92005000E-01 mass flux values 171 #ntotmfv= 146, nnodes= 50, number of f-m connections= 25 -0.12390000E-02 0.12390000E-02-0.12390000E-02 0. 0.12390000E-02 -0.12390000E-02 0. 0.12387000E-02-0.12387000E-02 0. 0.20896000E-04-0.20896000E-04 0. 0.59733000E-05-0.59733000E-05 0. 0.61827000E-06-0.61827000E-06 0. 0.20613000E-08 -0.20613000E-08 0. 0.98473000E-05-0.98473000E-05 0. 0.12246000E-02-0.12246000E-02 0. 0.12229000E-02-0.12229000E-02 0. 0.12229000E-02-0.12229000E-02 0. 0.12224000E-02 -0.12224000E-02 0. 0.12211000E-02-0.12211000E-02 0. 0.12196000E-02-0.12196000E-02 0. 0.12175000E-02-0.12175000E-02 0. 0.12159000E-02-0.12159000E-02 0. 0.11714000E-02 -0.11714000E-02 0. 0.10451000E-02-0.10451000E-02 0. 0.83611000E-03-0.83611000E-03 0. 0.83487000E-03-0.83487000E-03 0. 0.86962000E-03-0.86962000E-03 0. 0.86860000E-03 -0.86860000E-03 0. 0.11925000E-02-0.11925000E-02 0. 0.11749000E-02-0.11749000E-02 0.11520000E-02 0. 0.60741000E-08 -0.60741000E-08 0. 0.31222000E-07-0.31222000E-07 0. 0.33953000E-06-0.33953000E-06 0. 0.12181000E-02-0.12181000E-02 0. 0.12330000E-02-0.12330000E-02 0. 0.12384000E-02 -0.12384000E-02 0. 0.12390000E-02-0.12390000E-02 0. 0.12292000E-02-0.12292000E-02 0. 0.14363000E-04-0.14363000E-04 0. 0.16063000E-04-0.16063000E-04 0. 0.16102000E-04 -0.16102000E-04 0. 0.16639000E-04-0.16639000E-04 0. 0.17937000E-04-0.17937000E-04 0. 0.19370000E-04-0.19370000E-04 0. 0.21459000E-04-0.21459000E-04 0. 0.23079000E-04 -0.23079000E-04 0. 0.67622000E-04-0.67622000E-04 0. 0.19386000E-03-0.19386000E-03 0. 0.40289000E-03-0.40289000E-03 0. 0.40413000E-03-0.40413000E-03 0. 0.36938000E-03 -0.36938000E-03 0. 0.37040000E-03-0.37040000E-03 0. 0.46502000E-04-0.46502000E-04 0. 0.64082000E-04-0.64082000E-04 0.86966000E-04 0.60741000E-08 0.25148000E-07 0.30831000E-06 0.12178000E-02 0.14922000E-04 0.53551000E-05 0.61621000E-06-0.98452000E-05-0.12148000E-02 0.16997000E-05 0.39210000E-07 0.53745000E-06 0.12978000E-05 0.14333000E-05 0.20891000E-05 0.16192000E-05 0.44544000E-04 0.12624000E-03 0.20903000E-03 0.12428000E-05-0.34750000E-04 0.10209000E-05-0.32390000E-03 0.17580000E-04 0.22884000E-04 This file can be verified by comparing the saturations and mass fluxes to the actual values reported in the TOUGH2 output file (‘sd9_e9.ot1’ in DTN: SN9908T0581699.001). The first fracture element listed (‘FaE71’) is used to spot check these values. In ‘sd9_e9.ot1’ the liquid saturation is reported to be 0.70676E-01, which corresponds exactly to the saturation reported for the first fracture element in ‘fmsd9_e9.ini’ above. The mass flow rate between ‘FaE71’ and the second element below it (‘FbE71’) is reported in the TOUGH2 output file in CONNE as 0.12390E-02 kg/s. In the ‘fmsd9_e9.ini’ file this value can be found under the heading ‘mass flux values’ beneath the first header line. This value is actually the second number listed. The first number, which is identical in value but negative, represents the generation of mass flow originating from infiltration in the upper boundary element. The mass flow between the fracture and matrix elements corresponding to ‘FaE71’ can also be verified. In the TOUGH2 output file, the mass flow between ‘FaE71’ and ‘MaE71’ is given in CONNE as -0.60741E-08 kg/s (which denotes flow from the fracture to the matrix. The corresponding value can be found in ‘fmsd9_e9.ini’ by noting that there are 25 active nodes in each continuum. Therefore, this value should be the 25th value from the last number in the file. A visual check in ‘fmsd9_e9.ini’ shown above confirms that this value is correctly listed. Output File ‘*.rock’ This file lists the zones of all fracture and matrix materials. For each zone, the rock grain density (kg/m3), specific heat (J/kg-K), matrix porosity, and intrinsic fracture porosity (1) are listed. The output file ‘fmsd9_e9.rock’ is shown below, and values can be confirmed with the values in the ROCK card of ‘sd9_e9.dt1.’ "fmsd9_e9.rock" 124 lines, 8481 characters rock -1 0 0 0.2480E+04 0.1000E+04 0.6600E-01 -2 0 0 0.2480E+04 0.1000E+04 0.6600E-01 -3 0 0 0.2480E+04 0.1000E+04 0.1400E+00 -4 0 0 0.2300E+04 0.1000E+04 0.3690E+00 -5 0 0 0.2300E+04 0.1000E+04 0.2340E+00 -6 0 0 0.2300E+04 0.1000E+04 0.3530E+00 -7 0 0 0.2300E+04 0.1000E+04 0.4690E+00 -8 0 0 0.2300E+04 0.1000E+04 0.4640E+00 -9 0 0 0.2480E+04 0.1000E+04 0.4200E-01 -10 0 0 0.2480E+04 0.1000E+04 0.1460E+00 -11 0 0 0.2480E+04 0.1000E+04 0.1350E+00 -12 0 0 0.2480E+04 0.1000E+04 0.8900E-01 -13 0 0 0.2480E+04 0.1000E+04 0.1150E+00 -14 0 0 0.2480E+04 0.1000E+04 0.9200E-01 -15 0 0 0.2480E+04 0.1000E+04 0.2000E-01 -16 0 0 0.2300E+04 0.1000E+04 0.2650E+00 -17 0 0 0.2300E+04 0.1000E+04 0.3210E+00 -18 0 0 0.2300E+04 0.1000E+04 0.3210E+00 -19 0 0 0.2300E+04 0.1000E+04 0.3210E+00 -20 0 0 0.2300E+04 0.1000E+04 0.1930E+00 -21 0 0 0.2300E+04 0.1000E+04 0.2400E+00 -22 0 0 0.2300E+04 0.1000E+04 0.2400E+00 -23 0 0 0.2300E+04 0.1000E+04 0.1690E+00 -24 0 0 0.2300E+04 0.1000E+04 0.2740E+00 -25 0 0 0.2300E+04 0.1000E+04 0.1970E+00 -26 0 0 0.2300E+04 0.1000E+04 0.2740E+00 -27 0 0 0.2300E+04 0.1000E+04 0.1970E+00 -28 0 0 0.2300E+04 0.1000E+04 0.2740E+00 -29 0 0 0.2390E+04 0.1000E+04 0.2000E+00 -30 0 0 0.2390E+04 0.1000E+04 0.2000E+00 -31 0 0 0.2390E+04 0.1000E+04 0.2000E+00 -32 0 0 0.2390E+04 0.1000E+04 0.2000E+00 -33 0 0 0.2390E+04 0.1000E+04 0.2000E+00 -34 0 0 0.2390E+04 0.1000E+04 0.2000E+00 -35 0 0 0.2390E+04 0.1000E+04 0.2000E+00 -36 0 0 0.2390E+04 0.1000E+04 0.2000E+00 -37 0 0 0.2390E+04 0.1000E+04 0.2000E+00 -38 0 0 0.2390E+04 0.1000E+04 0.2000E+00 -39 0 0 0.2390E+04 0.1000E+04 0.2000E+00 -40 0 0 0.2390E+04 0.1000E+04 0.2000E+00 -41 0 0 0.2390E+04 0.1000E+04 0.2000E+00 -42 0 0 0.2390E+04 0.1000E+04 0.2000E+00 -43 0 0 0.2390E+04 0.1000E+04 0.2000E+00 -44 0 0 0.2390E+04 0.1000E+04 0.2000E+00 -45 0 0 0.2300E+04 0.1000E+04 0.2650E+00 -46 0 0 0.2300E+04 0.1000E+04 0.3210E+00 -47 0 0 0.2300E+04 0.1000E+04 0.2740E+00 -48 0 0 0.2300E+04 0.1000E+04 0.2650E+00 -49 0 0 0.2300E+04 0.1000E+04 0.3210E+00 -50 0 0 0.2300E+04 0.1000E+04 0.2740E+00 -51 0 0 0.2480E+04 0.1000E+04 0.3600E-01 -52 0 0 0.2300E+04 0.1000E+04 0.2880E+00 -53 0 0 0.2300E+04 0.1000E+04 0.2880E+00 -54 0 0 0.2300E+04 0.1000E+04 0.3320E+00 -55 0 0 0.2300E+04 0.1000E+04 0.3320E+00 -56 0 0 0.2300E+04 0.1000E+04 0.3320E+00 -57 0 0 0.2300E+04 0.1000E+04 0.2660E+00 -58 0 0 0.2300E+04 0.1000E+04 0.1000E+00 -59 0 0 0.2390E+04 0.1000E+04 0.5000E-01 -60 0 0 0.2390E+04 0.1000E+04 0.5000E-01 -61 0 0 0.2390E+04 0.1000E+04 0.5000E-01 -62 0 0 0.2390E+04 0.1000E+04 0.1000E+00 -63 0 0 0.2480E+04 0.1000E+04 0.1000E+01 -64 0 0 0.2480E+04 0.1000E+04 0.1000E+01 -65 0 0 0.2480E+04 0.1000E+04 0.1000E+01 -66 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -67 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -68 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -69 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -70 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -71 0 0 0.2480E+04 0.1000E+04 0.1000E+01 -72 0 0 0.2480E+04 0.1000E+04 0.1000E+01 -73 0 0 0.2480E+04 0.1000E+04 0.1000E+01 -74 0 0 0.2480E+04 0.1000E+04 0.1000E+01 -75 0 0 0.2480E+04 0.1000E+04 0.1000E+01 -76 0 0 0.2480E+04 0.1000E+04 0.1000E+01 -77 0 0 0.2480E+04 0.1000E+04 0.1000E+01 -78 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -79 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -80 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -81 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -82 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -83 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -84 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -85 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -86 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -87 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -88 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -89 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -90 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -91 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -92 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -93 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -94 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -95 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -96 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -97 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -98 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -99 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -100 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -101 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -102 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -103 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -104 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -105 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -106 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -107 0 0 0.2480E+04 0.1000E+04 0.1000E+01 -108 0 0 0.2480E+04 0.1000E+04 0.1000E+01 -109 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -110 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -111 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -112 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -113 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -114 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -115 0 0 0.2300E+04 0.1000E+04 0.1000E+01 -116 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -117 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -118 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -119 0 0 0.2390E+04 0.1000E+04 0.1000E+01 -120 0 0 0.2300E+04 0.1000E+04 0.2000E+00 -121 0 0 0.2300E+04 0.1000E+04 0.2800E+00 stop Output File ‘*.stor’ The file contains connectivity arrays and control volumes for the grid. Following two header lines, four integers are listed: iwtotl: Total number of connections in a continuum (either fracture or matrix) for which inter-node fluxes and areas are assigned. This includes connections for a node to itself for sources and sinks. Equal to ncont-(neq+1). neq: Number of nodes in either the fracture or matrix continuum. ncont: Number of values in the ncon array (see below) sehtemp: Flag that is equal to 1 for particle tracking The following arrays are then read from .stor: sx1(i), i=1,neq: Primary (total) volume of each node in a continuum (includes fracture and matrix) ncon(i), i=1,ncont: Node connectivity array that contains the node numbers for each connection to a specified node in one continuum, starting with node 1. The node numbers in ncon associated with connections to a given node include the node of interest. All nodes connected to a given node are listed in ascending order. In the beginning of this array is pointer information with neq+1 entries. The entries identify the index of the array (i=1,ncont) that precedes the node denoted by the index of the pointer information. See the figure below for an example of a 9-node network. 1 2 3 4 5 6 7 8 9 index i ncon(i) 1 10 2 13 3 17 4 20 5 24 6 29 7 33 8 36 9 40 10 43 11 1 12 2 13 4 14 1 15 2 16 3 17 5 18 2 19 3 20 6 21 1 22 4 23 5 24 7 25 2 26 4 27 5 28 6 29 8 30 3 31 5 32 6 33 9 34 4 35 7 36 8 37 5 38 7 39 8 40 9 41 6 42 8 43 9 Pointer Information Node 1 Node 2 Node 3 Node 4 Node 5 Node 6 Node 7 Node 8 Node 9 9-Node Example of the ncon Array Used in FEHM. istrw(i), i=1,ncont: Not used in this application. The array is filled using following algorithm: do i = 1, ncont if(i.le.iwtotl) then istrw(i) = i else istrw(i) = 0 end if end do nelmdg(i), i=1,neq: Position (index) of node i in the ncon array: do i = 1, neq do j = ncon(i) + 1, ncon(i+1) if (ncon(j).eq.i) nelmdg(i) = j end do end do iwtotl numbers: Three groups of iwtotl numbers signifying the x, y, and z components of the nodes are divided by distance terms for all internode connections. Only place-holders are required: do i = 1,3 write(15,’(5(1pe16.8))’) (-1.0, j=1, iwtotl) end do The file ‘fmsd9_e9.stor’ is shown below: "fmsd9_e9.stor" 98 lines, 6309 characters # input file for mean alpha, fitted fmx, present day q, ysw # AR 11/19/97 This is a .stor file with dummy area coefficients 73 25 99 1 4.93991416E+04 4.94314381E+04 3.02836879E+04 6.70661157E+04 8.54244306E+04 3.21230769E+04 1.43904899E+05 1.97772021E+05 4.78251121E+04 2.86899225E+05 5.73809524E+05 2.86935484E+05 1.67386018E+05 1.91306991E+05 2.86899696E+05 2.15197568E+05 2.15187970E+05 2.67276423E+03 1.84173669E+04 1.26272727E+05 1.48363636E+05 1.48363636E+05 1.45614035E+05 1.50909091E+05 1.90909091E+05 26 28 31 34 37 40 43 46 49 52 55 58 61 64 67 70 73 76 79 82 85 88 91 94 97 99 1 2 1 2 3 2 3 4 3 4 5 4 5 6 5 6 7 6 7 8 7 8 9 8 9 10 9 10 11 10 11 12 11 12 13 12 13 14 13 14 15 14 15 16 15 16 17 16 17 18 17 18 19 18 19 20 19 20 21 20 21 22 21 22 23 22 23 24 23 24 25 24 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 27 30 33 36 39 42 45 48 51 54 57 60 63 66 69 72 75 78 81 84 87 90 93 96 99 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 -1.00000000E+00 Output File ‘*.zone’ This file contains definitions of zones that correspond to ROCKS materials in TOUGH2. The materials are listed sequentially in the same order as they appear in the ROCKS card. A comment (#) is added to identify the name of the material as it appears in ROCKS. The number of nodes within each zone is listed after the header ‘nnum’. Following that line, the nodes are listed in the order that they appear in the ELEME card in TOUGH2. Note that both the fracture and matrix materials are listed in this file. Additional comments are added after the ‘stop’ line of the file. The first few lines of ‘fmsd9_e9.zone’ are shown below: zone 1 #tcwM1 nnum 1 26 2 #tcwM2 nnum 1 27 3 #tcwM3 nnum 1 28 Output File ‘*.zone2’ This file is identical to the .zone file except that it contains two additional zones that define the repository nodes for the fractures and matrix. The repository elements are listed in another file that is specified by the user during one of the prompts by T2FEHM2. This external file should contain the total number of repository elements in the file followed by a line-by-line listing of all the repository elements. This zone (.zone2) is read at the end of the .dat file to identify nodes where particles will be released in the ptrk macro (note that ptrk is not created by this postprocessor). The nodes that are defined in .zone2 will retain the porosities and densities assigned to them previously in .rock and ‘.dpdp.’ For the 1-D FEHM simulations in this analysis, only one repository element (‘FlE71’) is identified. The last few lines from ‘fmsd9_e9.zone2’ are shown below: 500 #fracture repository nodes nnum 1 12 501 #matrix repository nodes nnum 1 37 stop #Total number of nodes = 50 #Total number of active boundary materials = 2 #Total number of active boundary nodes = 2 Verification: The output from T2FEHM2 has been verified by visual inspection in the previous section that detailed the output files. This ensures that T2FEHM2 is performing correctly for the range of inputs used in this analysis. All files relevant to T2FEHM2 can be found in DTN: SN9908T0581699.001. Listing of Software Routine T2FEHM2 v. 2 c t2fehm2_v2.f C**************************************************************** c This program creates column formatted files from TOUGH2.OUT c files of EOS3 simulations. c Files MESH, TOUGH2.INP, and TOUGH2.OUT must be present. c The format of the output files are amenable for an FEHM c restart. c C.K.Ho 5/27/97 c This version now re-formats TOUGH2.OUT files in either EOS3 or c EOS9 format. Multidimensional files can be post-processed. This c version assumes that the elements listed in ELEME alternate c between fractures and matrix, starting with a fracture element. c This can be generalized in the loop (do 3000...) by knowing how c how the fracture and matrix elements were listed and by arranging c the arrays accordingly. I started this by asking the user to c specify the ordering, but I didn't do much with it in this version. c So for now, the elements should be listed alternately starting with c a fracture element. Also, the matrix materials are assumed to be listed c first in the ROCKS card. c c C.K.Ho c 9/2/97-9/12/97,9/19/97 c This version (op1postv3.f) is tailored specifically for LBL site-scale c runs. The previous version (option1postv2.f) is still good for SNL c TOUGH2 simulations of flow fields. The major revisions include reading c information from external files (MESH, GENER). In MESH, the material c identifier is a 5-character name--not an integer, which was assumed in c the previous version. The coordinates will have to be c read from MESH. Changes will have to be made for recognizing c fracture or matrix materials to accomodate all the materials (there c are greater than 100 materials) in the site-scale model. The dimensions c will have to be greatly increased to accommodate the 80,000 element c site-scale model. c C.K.Ho c 10/23/97 c c This version (op1postv4.f) does not assume any ordering in the ROCKS c card. There can be different numbers of matrix and fracture c materials written to the FEHM zone macro. Also, this version can read c in a file containing repository element names to create a separate zone. c Another assumption is that the active elements are listed before any c boundary elements ('TP' or 'BT') in ELEME. c C.K.Ho c 11/5/97 c c A few things have been cleaned up and it appears to work for the LBNL c 3-D site scale model. The current version is 't2fehm2.f'. c C.K.Ho c 11/6/97 c c This version accommodates new output formatting used by LBNL. The c index field in the output has been changed from i6 to i12. Also, c the flux output has been shifted to the left a bit, and nlin3 is now equal c to 3 instead of 4 (this is the amount of header lines inserted in the flux c output periodically). c The liquid pressure now appears where the gas pressure used to appear in c the output file. To calculate the gas pressure: Pg=Pl-Pc c C.K.Ho c 3/9/99 c***************************************************************** c23456789012345678901234567890123456789012345678901234567890123456789012 C implicit double precision (a-h,o-z) DIMENSION X(99000),Y(99000),z(99000),SL(99000),vol(99000) dimension PG(99000) dimension gelem(99000),ifm(99000) dimension fluxl(990000),fmlfm(99000),ncord(99000) dimension icon2(990000),flol2(990000),istrw(990000) dimension drok(500),por(500),nelmdg(99000),ncon2(99000) double precision lblpor CHARACTER*22 BLOCK CHARACTER*5 ELEMN(99000),ELEM1(490000),ELEM2(490000),ELEMX character*5 genname,matname(500),matb,mat(99000) character*80 header character*40 filen,control,dat,grid,ini,stor,dpdp,rock,zone character*40 filein,fileout,meshfile,repfile,zone2,check character*1 char2 character*5 repname(1003) common/int/ ncon(99000),icon(99000,35) common/flux/ flol(99000,35) C write(*,*) 'This program will re-format TOUGH2 output files' write(*,*)'for FEHM restart files. The following files' write(*,*)'must be present: input, output, and MESH.' write(*,*)'The MESH file should contain 5-character material' write(*,*)'names.' write(*,*) write(*,*)'What is the name of the input file?' read(*,*) filein write(*,*)'What is the name of the output file?' read(*,*) fileout write(*,*)'What is the name of the MESH file?' read(*,*) meshfile write(*,4) 4 format('What type of run is this?'/,'1) SNL EOS3'/,'2) SNL EOS9?'/ & ,'3) LBNL EOS9'/,'4) LBNL EOS9 SR/LA') read(*,*) neos write(*,*)'What reference name would you like to use for the' write(*,*)'FEHM restart files? (no spaces in the name)' read(*,*) filen write(*,*)'In ELEME, how are the elements listed?' write(*,*)'(1) Alternatively with matrix first' write(*,*)'(2) Alternatively with fracture first' write(*,*)'(3) All matrix, then all fractures' write(*,*)'(4) All fractures, then all matrix' read(*,*) norder write(*,*)'For fracture-matrix connections, which element is' write(*,*)'listed first: (1) Fracture or (2) Matrix?' read(*,*) nfmc write(*,*)'What is the print-out time (sec) of interest?' read(*,*) tsec write(*,*)'The fracture volumes will be used as the primary' write(*,*)'control volume for each element. Have they been' write(*,*)'modified in TOUGH2.INP? (1=yes, 0=no)' read(*,*) nvol volscale=1. if(nvol.eq.1) then write(*,*)'What is the scaling factor to retrieve correct', & ' primary volumes from fracture volumes?' read(*,*) volscale end if write(*,7) 7 format('What is the geometry?'/'0) 3-D'/'1) X-Y Plane'/ & '2) X-Z Plane'/'3) Y-Z Plane') read(*,*) icnl write(*,*)'Is there a file with repository element names?' write(*,*)'1 = yes, 0 = no' read(*,*) nrepans if(nrepans.eq.1) then write(*,9) 9 format('What is the name of the file with repository elements?') read(*,*) repfile write(*,*)'Would you like to modify the 2nd character of the' write(*,*)'element name? 1=yes, 0=no' read(*,*) n2nd if(n2nd.eq.1) then write(*,*)'What character would you like to use?' read(*,'(a1)') char2 end if open(19,file=repfile,status='old') end if if(norder.eq.1.or.norder.eq.2) then nalt=2 else nalt=1 end if c...Define FEHM restart files based on reference name kend=index(filen,' ') control=filen(1:kend-1)//'.files' dat=filen(1:kend-1)//'.dat' grid=filen(1:kend-1)//'.grid' ini=filen(1:kend-1)//'.ini' stor=filen(1:kend-1)//'.stor' dpdp=filen(1:kend-1)//'.dpdp' rock=filen(1:kend-1)//'.rock' zone=filen(1:kend-1)//'.zone' zone2=filen(1:kend-1)//'.zone2' check=filen(1:kend-1)//'.check' if(neos.eq.1) then nlin1=5 nlin2=3 nlin3=3 elseif(neos.eq.2) then nlin1=6 nlin2=4 nlin3=4 elseif(neos.eq.3) then nlin1=6 nlin2=3 nlin3=4 elseif(neos.eq.4) then nlin1=6 nlin2=3 nlin3=3 end if write(*,*) 'Thank You! Please wait while I work...' open(1,file=meshfile,status='old') open(2,file=fileout,status='old') open(3,file=filein,status='old') open(11,file=control,status='unknown') open(12,file=dat,status='unknown') open(13,file=grid,status='unknown') open(14,file=ini,status='unknown') open(15,file=stor,status='unknown') open(16,file=dpdp,status='unknown') open(17,file=rock,status='unknown') open(18,file=zone,status='unknown') open(22,file=check,status='unknown') open(23,file=zone2,status='unknown') c....Data spht=1.e3 per1=1.e-15 per2=1.e-15 per3=1.e-15 day=365.25e6 tims=365.25e6 nstep=10 iprtout=10 iyear=1997 month=10 maxit=-10 epm=1.e-4 north=40 ja=1 jb=0 jc=0 igaus=1 as=1. grav=3. upwgt=1. iamm=5 aiaa=2. daymin=1.e-10 daymax=1.e10 lda=1 g1=1.e-5 g2=1.e-5 g3=1.e-5 tmch=-1.e-4 overf=1.2 irdof=0 islord=0 iback=0 icoupl=0 rnmax=14400. ntt=1 intg=-1 zero=1.d-10 ra=287. rv=461.52 C c...Read header from TOUGH2.INP read(3,'(a80)') header c_______________________________________________________________ c...Write information to .dat file write(12,510) header 510 format(a80/'# Particle tracking for TOUGH2 flow field') c...Write dpdp macro write(12,516) dpdp 516 format('dpdp'/'file'/a) c...Write perm macro write(12,518) per1,per2,per3 518 format('perm'/'1 0 0 ',3e10.3/) c...Write rlp macro write(12,520) 520 format('rlp'/'1 0. 0. 1. 1. 0. 1.'//'1 0 0 1'/) c...Write rock macro write(12,522) rock 522 format('rock'/'file'/a) c...Write flow macro write(12,524) 524 format('flow'/) c...Write time macro write(12,526) day,tims,nstep,iprtout,iyear,month 526 format('time'/2e13.5,4i8/) c...Write ctrl macro write(12,528) maxit,epm,north,ja,jb,jc,igaus,as,grav,upwgt, & iamm,aiaa,daymin,daymax,icnl,lda 528 format('ctrl'/i8,e10.2,i8/4i8/'0'/3f10.2/i8,3e10.2/2i8) c...Write iter macro write(12,530) g1,g2,g3,tmch,overf,irdof,islord,iback,icoupl, & rnmax 530 format('iter'/5e10.2/4i8,e10.2) c...Write sol macro write(12,532) ntt,intg 532 format('sol'/2i8) c...Write rflo macro write(12,534) 534 format('rflo'/'air'/'-1'/'20.0 0.1') c...Write node macro write(12,536) 536 format('node'/'1'/'1') c...Write zone macro that corresponds to the repository nodes write(12,515) zone2 515 format('zone'/'file'/a) c...Write ptrk macro write(12,538) filen(1:kend-1) 538 format('ptrk'/'file'/a,'.ptrk') c...Write stop write(12,540) 540 format('stop') c_______________________________________________________________ c...Write information to control file write(11,501) dat,grid,zone,filen(1:kend-1),ini,filen(1:kend-1) &,filen(1:kend-1),filen(1:kend-1),filen(1:kend-1),stor, &filen(1:kend-1) 501 format(a/a/a/a,'.out'/a/a,'.fin'/a,'.his'/a,'.trc'/a,'.con'// & a/a,'.chk'/'all'/'0') c...Read in repository element names if(nrepans.eq.1) then read(19,*) nrepelem numrep=nrepelem do i=1,nrepelem read(19,'(a5)') repname(i) repname(i)(1:1)='F' if(n2nd.eq.1) repname(i)(2:2)=char2 end do end if c...Read in grid information from MESH nbelm=0 nbmat=0 matb=' ' N=1 read(1,1000) block 1000 format(a22) 99 read(1,65) elemn(n),mat(n),vol(n),x(n),y(n),z(n) 65 format(a5,10x,a5,e10.4,20x,3e10.4) if(elemn(n).eq.' ') go to 98 if(elemn(n)(4:4).eq.'0') elemn(n)(4:4)=' ' c...Count number of boundary elements, nbelm, and number of boundary c...materials, nbmat. if(elemn(n)(1:2).eq.'TP'.or.elemn(n)(1:2).eq.'BT') then nbelm=nbelm+1 if(mat(n).ne.matb) then nbmat=nbmat+1 matb=mat(n) end if end if N=N+1 GO TO 99 98 CONTINUE NMAX = N - 1 c...NMAX is the total number of elements read from MESH write(*,107) nmax 107 format('Have read in ',i8,' elements from MESH...') c...nnodes is the total number of active nodes nnodes=nmax-nbelm c...Find maximum number of materials used in ROCKS (nmat) c nmat=0 c do i=1,nmax c nmat=max(mat(i),nmat) c end do c write(*,222) nmat c222 format('Maximum number of active materials = ',i8,'...') c...nfmat is the number of fracture materials c nfmat=(nmat-nbmat)/2 c...Read in connection information from MESH N=1 READ(1,1500) BLOCK 1500 FORMAT(A22,3X,25X,E10.4) 199 read(1,1502) elem1(n),elem2(n),ifm(n) c...ifm(n) is a flag in the 75th column of the CONNE card that Yu-Shu has c...specified as equal to '2' for fracture-matrix connections 1502 format(2a5,64x,i1) IF(elem1(n)(1:5).EQ.' '.OR.elem1(n)(1:3).EQ.'+++') GO TO 198 if(elem1(n)(4:4).eq.'0') elem1(n)(4:4)=' ' if(elem2(n)(4:4).eq.'0') elem2(n)(4:4)=' ' N=N+1 GO TO 199 198 CONTINUE NCMAX = N - 1 c...NCMAX is the total number of connections read from MESH write(*,203) ncmax 203 format('Have read in ',i8,' connections from MESH...') c...Read in ROCKS information from TOUGH2 input file 18 read(3,1000) block if(block(1:5).ne.'ROCKS') go to 18 i=1 nfmat=0 nmmat=0 408 read(3,410) matname(i),drok(i),por(i) 410 format(a5,5x,2e10.4) if(matname(i).eq.'REFCO') go to 408 if(matname(i).eq.' ') then c...ntotmat is the total number of materials in the ROCKS card c...nmat is the number of materials associated with non-boundary c...elements ntotmat=i-1 nmat=ntotmat-nbmat go to 27 end if c...LBNL uses columns 71-80 in the second line of each material card to c...identify the fracture porosity read(3,415) lblpor 415 format(70x,e10.4) c...nfmat is the total number of fracture materials if(matname(i)(3:3).eq.'F'.or.matname(i)(4:4).eq.'F') then nfmat=nfmat+1 if(neos.eq.3.or.neos.eq.4) por(i)=lblpor c...The perched water fractures do not have porosities listed in ROCKS. c...Yu-Shu said that they have the same porosity as the zeolitic fractures, c...which is 1.1e-5 (phone message 10/31/97). if(por(i).eq.0.) por(i)=1.1d-5 end if c...nmmat is the total number of matrix materials if(matname(i)(3:3).eq.'M'.or.matname(i)(4:4).eq.'M') nmmat=nmmat+1 read(3,*) read(3,*) i=i+1 go to 408 27 continue c...10/27/97 Ho c...Write grid macro file write(13,202) nnodes/2 202 format('coor'/i8) c...This assumes that all boundary elements ('TP' and 'BT') are listed c...after the active elements in ELEME do i=1,nnodes/2 write(13,204) i,x(i*nalt),y(i*nalt),z(i*nalt) 204 format(i8,3(3x,f10.2)) end do write(13,206) 206 format(/'elem'/'2 1'/'1 2 1'//'stop') c...Initialize generation array do i=1,nmax gelem(i)=0. end do c...Read in generation information from TOUGH2.INP i=1 33 read(3,1000,end=299) block if(block(1:5).ne.'GENER') go to 33 74 read(3,75) genname,g 75 format(a5,35x,e10.4) if(genname.eq.' ') go to 77 if(genname(4:4).eq.'0') genname(4:4)=' ' do ik=1,nmax if(genname.eq.elemn(ik)) then c...Assign a generation term for each element (flow into an element c...is defined as negative) c...The method used here is different than in v3. It eliminates a c...separate do-loop and the need for arrays igen and g. gelem(ik)=-g i=i+1 go to 74 end if end do write(*,*)'Could not find element name for generation' write(*,79) i,genname 79 format('element ',i8,': ',a5) stop 299 write(*,*)'***Warning*** No generation card in TOUGH2.INP' 77 ngentot=i-1 c...Write zone macro ntotin=0 write(18,'(a4)') 'zone' write(23,'(a4)') 'zone' do i=1,ntotmat write(18,512) i,matname(i) write(23,512) i,matname(i) 512 format(i4,5x,'#',a5) write(18,'(a4)') 'nnum' write(23,'(a4)') 'nnum' nin=1 do j=1,nmax c...Match nodes to respective materials. This assumes that the c...fractures and matrix elements are listed alternately in ELEME c...starting with the fractures first c...If element is a boundary element, go to next element if(elemn(j)(1:2).eq.'TP'.or.elemn(j)(1:2).eq.'BT') goto 517 if(mat(j).eq.matname(i)) then if(mat(j)(3:3).eq.'F'.or.mat(j)(4:4).eq.'F') then ncord(nin)=(j+1)/nalt nin=nin+1 go to 517 end if if(mat(j)(3:3).eq.'M'.or.mat(j)(4:4).eq.'M') then ncord(nin)=j/nalt+nnodes/2. nin=nin+1 end if end if 517 end do nin=nin-1 ntotin=ntotin+nin write(18,'(i10)') nin write(23,'(i10)') nin if(nin.gt.0) write(18,'(8i10)') (ncord(k),k=1,nin) if(nin.gt.0) write(23,'(8i10)') (ncord(k),k=1,nin) end do write(18,*) write(18,'(a4)') 'stop' c...Now write zones for nodes corresponding to repository elements nrp=1 do i=1,nmax do j=1,numrep if(elemn(i).eq.repname(j)) then ncord(nrp)=(i+1)/nalt nrp=nrp+1 go to 527 end if end do 527 end do nrp=nrp-1 write(23,*) '500 #fracture repository nodes' write(23,'(a4)') 'nnum' write(23,'(i10)') nrp if(nrp.gt.0) write(23,'(8i10)') (ncord(k),k=1,nrp) write(23,*) '501 #matrix repository nodes' write(23,'(a4)') 'nnum' write(23,'(i10)') nrp do i=1,nrp ncord(i)=ncord(i)+nnodes/2. end do if(nrp.gt.0) write(23,'(8i10)') (ncord(k),k=1,nrp) write(23,*) write(23,'(a4)') 'stop' c...Now write some additional information to the zone file write(18,*) write(23,*) write(18,514) ntotin,nbmat,nbelm write(23,514) ntotin,nbmat,nbelm 514 format(/'#Total number of nodes = ',i8/'#Total number of', & ' active boundary materials = ',i8/'#Total number of active', & ' boundary nodes = ',i8/) c...Write dpdp macro file write(16,550) 550 format('dpdp'/'1') c...Loop over the materials and print out fracture porosities do i=1,ntotmat if(matname(i)(3:3).eq.'F'.or.matname(i)(4:4).eq.'F') then write(16,552) -i,jb,jc,por(i) 552 format(3i8,5x,e10.4) end if end do write(16,554) ja,jb,jc 554 format(/,3i8,5x,'99.'//'stop') c...Write rock macro file write(17,556) 556 format('rock') do i=1,ntotmat porock=por(i) if(matname(i)(3:3).eq.'F'.or.matname(i)(4:4).eq.'F')porock=1. write(17,558) -i,jb,jc,drok(i),spht,porock 558 format(3i8,5x,e10.4,5x,e10.4,5x,e10.4) end do write(17,559) 559 format(/'stop') c...Search for "TOTAL TIME" in TOUGH2.OUT and then read in variables 89 READ(2,1000,END=90) BLOCK IF(BLOCK(1:12).NE.' TOTAL TIME') GO TO 89 READ(2,1001) TIME if(time.ne.tsec.and.tsec.gt.0) go to 89 1001 FORMAT(E13.5) do nl=1,nlin1 READ(2,1000) BLOCK end do C c23456789012345678901234567890123456789012345678901234567890123456789012 c...Read in state variables from TOUGH2.OUT 115 N1=1 N2=MIN(NMAX,45) DO 2000 I=N1,N2 if(neos.eq.1) then c... This is EOS3 format READ(2,1002) PG(I),SL(I) 1002 FORMAT(12x,e12.5,24x,7e12.5) elseif(neos.eq.2.or.neos.eq.3) then c... This is EOS9 format read(2,118) pg(i),sl(i) 118 format(12x,2e12.5) elseif (neos.eq.4) then c... This is EOS9 format with new index formatting of i12 read(2,119) pl,sl(i),pc pg(i)=pl-pc 119 format(18x,3e12.5) end if 2000 CONTINUE C 2100 CONTINUE c...Check to see if we've read in all the element variables IF(N2.EQ.NMAX) GO TO 91 N1=N2+1 N2=MIN(NMAX,N1+56) do nl=1,nlin2 READ(2,1000) BLOCK end do DO 2010 I=N1,N2 if(neos.eq.1) then c... This is EOS3 format READ(2,1002) PG(I),SL(I) elseif(neos.eq.2.or.neos.eq.3) then c... This is EOS9 format read(2,118) pg(i),sl(i) elseif (neos.eq.4) then c... This is EOS9 format with new index formatting of i12 read(2,119) pl,sl(i),pc pg(i)=pl-pc end if 2010 CONTINUE GO TO 2100 C 91 CONTINUE C c...Write saturations to .ini file (fractures saturations first followed c...by matrix saturations) write(14,302) header 302 format(a80/'This is a .ini file with saturations, pressures', & ' and mass flux values.'/'0.'/'air'/'ptrk'/'nstr'/ & 'dpdp'/'ndua') write(14,304) (sl(i),i=1,nnodes,2),(sl(i),i=2,nnodes,2) 304 format(4g16.8) c...Write pressures to .ini file in MPa (fractures first, then matrix) write(14,304) (pg(i)*1.d-6,i=1,nnodes,2), & (pg(i)*1.d-6,i=2,nnodes,2) write(*,*)'Have read in state variables from output file...' C c...Read in flux variables from TOUGH2.OUT 289 READ(2,1500,END=190) BLOCK if(neos.lt.4) then IF(BLOCK(11:22).NE.'ELEM1 ELEM2') GO TO 289 elseif (neos.eq.4) then IF(BLOCK(7:18).NE.'ELEM1 ELEM2') GO TO 289 end if READ(2,1500) BLOCK READ(2,1500) BLOCK C c...Read in mass flow liquid for each connection pair N1=1 N2=MIN(NCMAX,53) DO 1600 I=N1,N2 if(neos.eq.1) then READ(2,1003) fluxl(I) 1003 FORMAT(80x,4e13.5) elseif (neos.eq.2.or.neos.eq.3) then read(2,121) fluxl(i) 121 format(29x,e13.5) elseif (neos.eq.4) then read(2,122) fluxl(i) 122 format(31x,e13.5) end if 1600 CONTINUE C 2150 CONTINUE IF(N2.EQ.NCMAX) GO TO 191 N1=N2+1 N2=MIN(NCMAX,N1+56) do nl=1,nlin3 READ(2,1500) BLOCK end do DO 2020 I=N1,N2 if(neos.eq.1) then READ(2,1003) fluxl(I) elseif (neos.eq.2.or.neos.eq.3) then read(2,121) fluxl(i) elseif (neos.eq.4) then read(2,122) fluxl(i) end if 2020 CONTINUE GO TO 2150 C 191 CONTINUE C 190 CONTINUE c...Check write(*,*)'Have read in flux variables from output file...' c...Check c do i=1,ncmax c write(15,444) i,elem1(i),elem2(i),fluxl(i) c444 format(i8,2x,2(a5,2x),e10.4) c end do c stop c...End check C c...Loop over all elements to determine connections and fluxes for each c...element nmlfm=1 c...nmlfm is the total number of fracture-matrix connections DO 3000 I=1,NMAX if(mod(i,1000).eq.0) write(*,472) i 472 format('Still working... Element ',i8) c...fmlfm(i) is the flow (kg/s) between fracture and matrix fmlfm(i)=0.d0 c...jj is the number of connections for each element do jj=1,35 flol(i,jj)=0.d0 c...icon(i,jj) is the node number of the element for connection jj to element i icon(i,jj)=0 end do ELEMX=ELEMN(I) c...If element is a boundary element, go to next element if(elemx(1:2).eq.'TP'.or.elemx(1:2).eq.'BT') go to 3000 c...Write the element number and the number of connections for that element if(i.gt.1) write(22,*) i-1,ncon(i-1) c_________________________________________________________________________ c...For each element, loop over all connections to determine if c...the element is either the first or second element in each connection c...nc is the number of connections per element nc=1 DO 3001 J=1,NCMAX c...Say element is the first element in the connection if(elem1(j).eq.elemx) then nsign=-1 c...If connecting element is the top boundary, go to next connection if(elem2(j)(1:2).eq.'TP') go to 3001 c...If connecting element is the bottom boundary, treat the flow to the c...bottom boundary as a sink/source term and move on to the next connection if(elem2(j)(1:2).eq.'BT') then gelem(i)=fluxl(j)*nsign go to 3001 end if c...What is the second element in the connection? do ii=1,nmax if(elem2(j).eq.elemn(ii)) then k2nd=ii c...Determine if the connection is between a fracture and matrix element c...If it is a fracture-matrix connection (both elements have the same c...coordinates, or ifm=2), store this flux separately from fracture-fracture c...or matrix-matrix fluxes. dx=dabs(x(k2nd)-x(i)) dy=dabs(y(k2nd)-y(i)) dz=dabs(z(k2nd)-z(i)) if(dx.le.zero.and.dy.le.zero.and.dz.le.zero.or. & ifm(j).eq.2) then c...If the first element of f-m connection is a fracture, then process this if(nfmc.eq.1) then go to 3017 else go to 3001 end if end if icon(i,nc)=ii flol(i,nc)=fluxl(j)*nsign nc=nc+1 go to 3002 endif end do write(*,7001) elemx,j,elem2(j),elem2(j-1),elem2(j+1) 7001 format('***Could not find 2nd element in connection for', & ' first element ',a5,'***'/'Connection index = ',i8/ & 'Second element = ',a5/'j-1= ',a5/'j+1= ',a5) stop end if c...If no match in first element of connection, try second element if(elem2(j).eq.elemx) then nsign=1 c...If connecting element is the top boundary, go to next connection if(elem1(j)(1:2).eq.'TP') go to 3001 c...If connecting element is the bottom boundary, treat the flow to the c...bottom boundary as a sink/source term and move on to the next connection if(elem1(j)(1:2).eq.'BT') then gelem(i)=fluxl(j)*nsign go to 3001 end if c...What is the first element in the connection? do ii=1,nmax if(elem1(j).eq.elemn(ii)) then k2nd=ii c...Determine if the connection is between a fracture and matrix element c...If it is a fracture-matrix connection (both elements have the same c...coordinates), store this flux separately from fracture-fracture or c...matrix-matrix fluxes. dx=dabs(x(k2nd)-x(i)) dy=dabs(y(k2nd)-y(i)) dz=dabs(z(k2nd)-z(i)) if(dx.le.zero.and.dy.le.zero.and.dz.le.zero.or. & ifm(j).eq.2) then c...If the second element of f-m connection is a fracture, then process this if(nfmc.eq.2) then go to 3017 else go to 3001 end if end if icon(i,nc)=ii flol(i,nc)=fluxl(j)*nsign nc=nc+1 go to 3002 end if end do write(*,7000) elemx,j,elem1(j) 7000 format('***Could not find 1st element in connection for', & ' second element ',a5,'***'/'Connection index = ',i8/ & '1st element = ',a5) stop end if c...If neither element 1 or 2 for connection j is equal to elemx, then c...go on to the next connection goto 3001 3002 continue c______________________________________________________________________ c...go to next connection go to 3001 c______________________________________________________________________ c...Come here if this is a fracture-matrix connection AND the element c...being considered (elemx=elemn(i)) is a fracture c...Consider outflow to be positive and c...that the first element in the connection is a fracture 3017 continue fmlfm(nmlfm)=nsign*fluxl(j) nmlfm=nmlfm+1 c...Go to next connection c______________________________________________________________________ 3001 continue c...ncon(i) is the total number of connections for node i ncon(i)=nc-1 C c...Check c write(15,446) i,ncon(i),(icon(i,j),j=1,ncon(i)) c446 format(10(i8,2x)) c write(15,448) i,ncon(i),(flol(i,j),j=1,ncon(i)) c448 format(2(i8,2x),8(e10.4,2x)) c...End check c...Go to next element 3000 CONTINUE c...nmlfm is the total number of fracture-matrix connections nmlfm=nmlfm-1 c...Add connection for each element to itself using generation array c...nmfluxval is the total number of mass flux values c...Note: nodes 1-nnodes are still assumed to alternative between c...fractures and matrix. This will be adjusted later in the print-out c...to the FEHM files. nmfluxval=0 do i=1,nnodes ncon(i)=ncon(i)+1 icon(i,ncon(i))=i flol(i,ncon(i))=gelem(i) nmfluxval=nmfluxval+ncon(i) c...Check c write(15,448) i,ncon(i),flol(i,ncon(i)),nmfluxval c448 format(2(i8,2x),e10.4,2x,i8) c...End check c...nmfluxval is the total number of flux values for fracture and matrix c...elements excluding f-m fluxes end do c...Call sort subroutine to sort the necessary arrays in ascending order c...of elements for each connection pair of a given element call sort(nnodes) C c...Create 1-D arrays containing icon and flol information. The arrays c...will be icon2 and flol2. This assumes that the fractures and matrix c...elements alternate in ELEME and fractures are listed first. k=1 jj=1 ncont1=0 c...ncont1 is the total number of connections for each continuum c...do the fracture continuum first do i=1,nnodes,2 do j=1,ncon(i) c...The index k+nnodes/2+1 accounts for the leading pointer information icon2(k+nnodes/2+1)=(icon(i,j)+1)/2 flol2(k)=flol(i,j) k=k+1 end do ncont1=ncont1+ncon(i) c...ncon2(jj) is the number of connections for fracture node jj, where jj is c...now icremented 1,2,3...nnodes/2 ncon2(jj)=ncon(i) jj=jj+1 end do c...Now do the matrix continuum do i=2,nnodes,2 do j=1,ncon(i) flol2(k)=flol(i,j) k=k+1 end do end do c...ntotmfv is the total number of connections. This can be compared to c...nmfluxval as a cross-check to see if they're equal. ntotmfv=k-1 c...Write mass flux values to .ini file write(14,602) nmlfm+nmfluxval,ntotmfv,nnodes,nmlfm 602 format('mass flux values'/i8,5x,'#ntotmfv=',i8,', nnodes=',i8, & ', number of f-m connections= ',i8) write(14,604) (flol2(i),i=1,ntotmfv),(fmlfm(i),i=1,nmlfm) 604 format(5g15.8) c...Write .stor file write(15,702) header 702 format(a80/'This is a .stor file with dummy area coefficients') c...Add the pointer information (number of fracture nodes+1) to ncont1 neq=nnodes/2 ncont=ncont1+(neq+1) iwtotl=ncont-(neq+1) write(15,704) iwtotl,neq,ncont,1 704 format(4(i8,2x)) c...Write primary volume for each node to .stor c...If this is an LBNL run, then divide the fracture volumes by the c...fracture porosity, since the volumes in ELEME were multiplied by c...the fracture porosity. if(neos.eq.3.or.neos.eq.4) then do i=1,nnodes,2 do j=1,ntotmat if(mat(i).eq.matname(j)) then vol(i)=vol(i)/por(j) go to 833 end if end do 833 end do end if c...If the fracture volumes were globally modified, multiply the volume c...by a scaling factor, volscale, specified by the user to get the original c...volume back. write(15,706) (vol(i)*volscale,i=1,nnodes,2) 706 format(1p5e16.8) c...Compile and write ncon and pointer information c...Fill the icon2(i) array from i=1,neq+1 (recall that icon2(i) has c...already been filled from neq+2 to ncont1 (the total number of connections c...for the fracture continuum icon2(1)=neq+1 do i=2,neq+1 icon2(i)=icon2(i-1)+ncon2(i-1) end do write(15,708) (icon2(i),i=1,ncont) 708 format(5(i8,2x)) c...Compile and write istrw information to .stor file do i=1,ncont if(i.le.iwtotl) then istrw(i)=i else istrw(i)=0 end if end do write(15,708) (istrw(i),i=1,ncont) c...Compile and write nelmdg information to .stor file do i=1,neq do j=icon2(i)+1,icon2(i+1) if(icon2(j).eq.i) nelmdg(i)=j end do end do write(15,708) (nelmdg(i),i=1,neq) c...Write dummy area coefficients to .stor file do i=1,3 write(15,706) (-1.0,j=1,iwtotl) end do c______________________________________________________________________ write(*,1153) time 1153 format('Finished processing printout at ',e12.4,' sec') go to 722 C 90 CONTINUE write(*,*)'**Did not find desired print-out time in TOUGH2.OUT**' C 722 write(*,*) 'Done!!!' stop END subroutine sort(nnodes) c______________________________________________________________________ c This subroutine sorts variables using a multipass method. c C.K.Ho c 9/8/97 c______________________________________________________________________ implicit double precision (a-h,o-z) common/int/ ncon(99000),icon(99000,35) common/flux/ flol(99000,35) c...The objective here is to arrange the connections in ascending order c...of connecting node number. The associated flux should also be sorted. nsort=1 do i=1,nnodes 5 if(nsort.eq.1) then nsort=0 do j=1,ncon(i)-1 if(icon(i,j).gt.icon(i,j+1)) then itempicon=icon(i,j) icon(i,j)=icon(i,j+1) icon(i,j+1)=itempicon tempflol=flol(i,j) flol(i,j)=flol(i,j+1) flol(i,j+1)=tempflol nsort=1 end if end do go to 5 end if nsort=1 end do return end