Particle Tracking Model and Abstraction of Transport Processes Rev 00, ICN 00 ANL-NBS-HS-000026 April 2000 1. PURPOSE The purpose of the transport methodology and component analysis is to provide the numerical methods for simulating radionuclide transport and model setup for transport in the unsaturated zone (UZ) site-scale model. The particle-tracking method of simulating radionuclide transport is incorporated into the FEHM computer code and the resulting changes in the FEHM code are to be submitted to the software configuration management system. This Analysis and Model Report (AMR) outlines the assumptions, design, and testing of a model for calculating radionuclide transport in the unsaturated zone at Yucca Mountain. In addition, methods for determining colloid-facilitated transport parameters are outlined for use in the Total System Performance Assessment (TSPA) analyses. Concurrently, process-level flow model calculations are being carried out in a PMR for the unsaturated zone (CRWMS M&O 2000a). The computer code TOUGH2 (Pruess 1991) is being used to generate three-dimensional, dualpermeability flow fields, that are supplied to the Performance Assessment group for subsequent transport simulations. These flow fields are converted to input files compatible with the FEHM code, which for this application simulates radionuclide transport using the particle-tracking algorithm outlined in this AMR. Therefore, this AMR establishes the numerical method and demonstrates the use of the model, but the specific breakthrough curves presented do not necessarily represent the behavior of the Yucca Mountain unsaturated zone. The particle-tracking technique presented in this AMR, called the Residence Time Transfer Function (RTTF) particle-tracking technique, uses a cell-based approach that sends particles from node to node on a finite difference or finite element grid, after keeping each particle at the cell for a prescribed period of time. To incorporate transport mechanisms such as dispersion and matrix diffusion, the residence time of a particle at a cell is computed using a transfer function that ensures that the correct distribution of times at the cell is reproduced. This procedure is computationally very efficient, enabling large-scale transport simulations of several million particles to be completed rapidly on modern workstations. This requirement was needed for complex, three-dimensional simulation involving the simulation of multiple radionuclides. Furthermore, since the cell-based approach uses directly mass flow rate information generated from a numerical fluid flow solution, complex, unstructured computational grids and the dualpermeability flow model formulation pose no additional complications. For the present application, the technique is used for unsaturated, dual-permeability transport simulations, for which the method is well suited. The cell-based approach allows accurate simulation of dualpermeability systems in which there is a vast disparity in the travel times depending on whether the transport is in the fractures or the matrix. Furthermore, matrix diffusion and colloidfacilitated radionuclide transport can be simulated. Additionally, complex source terms and decay chain/ingrowth capabilities have been included in the model. Like all numerical methods, the particle-tracking technique has limitations that must be considered when deciding whether its use is appropriate for a given application. The key physical and chemical assumptions are advection-dominated transport and linear, equilibrium sorption. The possibility of grid orientation effects should also be considered. These effects can manifest themselves as an artificial lateral spreading of solute mass. Also, the accuracy of the ANL-NBS-HS-000026, Rev 00 12 April 2000 method for dual-permeability flow systems was investigated in detail and found to perform best when the flow regime undergoes abrupt transitions at unit interfaces, and in cases for relatively low diffusion. Given these results, this AMR demonstrates that the particle-tracking model can be used in three-dimensional radionuclide transport simulations of the Yucca Mountain unsaturated zone as long as the limits on the model are recognized and parameters are chosen accordingly. However, breakthrough curves presented in this report are not necessarily representative of the Yucca Mountain unsaturated zone. Also, given the accuracy of the model without diffusion, supporting particle-tracking model runs in the absence of diffusion should be performed to access the importance of matrix diffusion to the overall conclusions of the performance assessment. This analysis is governed by the following Office of Civilian Radioactive Waste Management (OCRWM) Work Direction and Planning Documents: “Improve Documentation and Verification of Radionuclide Decay Model (Rev. 01)” (CRWMS M&O 1999a); “Enhance Particle Tracking Features for TSPA (Rev. 01)” (CRWMS M&O 1999b); “Improve Matrix Diffusion Model (Rev. 01)” (CRWMS M&O 1999c). 2. QUALITY ASSURANCE The QA program applies to the development of this AMR. The Performance Assessment Operations (PAO) responsible manager has evaluated this activity in accordance with QAP-2-0, Conduct of Activities. The QAP-2-0 activity evaluation Conduct of Performance Assessment (CRWMS M&O 1999d), has determined that the preparation and review of this technical document of is subject to the Quality Assurance Requirements and Description (QARD) DOE/RW-0333P (DOE 2000) requirements. Preparation of this analysis/model did not require classification of items in accordance with QAP-2-3, Classification of Permanent Items. 3. COMPUTER SOFTWARE AND MODEL USAGE Table 1 lists the software and routines used in this study. The computer software code used as a starting point to perform/model the saturated-zone particle tracking in this AMR is FEHM, V2.0, CSCI/ STN# 10031-2.00-00 for a SUN Ultra Sparc, and is under Configuration Management (CM) control. A revised version of the code that includes the particle-tracking algorithm described in this AMR is currently being qualified under AP-SI.1Q, Software Management, as FEHM, Version 2.10 (Software Activity Number (SAN) LANL-1999-046, Software Tracking Number (STN): 10086-2.10-00). Version 2.10 is the version used for all calculations performed in this AMR; therefore, the FEHM software is currently unqualified and all results are to be considered TBV. The software is appropriate for the application and was used only within the range for which it was developed. The input and output files for these analyses are being submitted to the Project database for archival as DTN: LA0002BR12213S.002. ANL-NBS-HS-000026, Rev 00 13 April 2000 Table 1. Computer Software and Routines Software Name Version Software Tracking Number (STN) Computer Platform FEHM 2.10 10086-2.10-00 SUN with UNIX OS, FORTRAN Routines: Documentation BKPM 1.0 Attachment II SUN with UNIX OS, FORTRAN CHAIN 1.0 Attachment III PC with Windows NT, FORTRAN TEST 1.0 Attachment I SUN with UNIX OS, FORTRAN test_ebs_random 1.0 Attachment IV SUN with UNIX OS, FORTRAN In addition, the following commercially available software was used in this analysis and documentation: FORTNER PLOT SUN Workstations Version 1.3 This software was used for plotting graphs. Therefore, the software was used for presentation purposes only, and there were no additional routines or macros developed using this commercial software. 4. INPUTS 4.1 DATA AND PARAMETERS For the code development and testing portions of this AMR, generic values of hydrologic and transport parameters are used that are similar to those to be used in the unsaturated-zone flow and transport modeling activities. However, this code testing does not require referenceable parameter values to verify that the model is correctly implemented. For other site-specific models developed in this AMR, the following data sources were used: Table 2. Input Data Used DTN Data Description LB990501233129.001 Mean fracture aperture and spacing, variance in aperture GS980908312242.039 GS950608312231.008 GS960808312231.003 Moisture retention curves – used to calculate pore size distributions LA0003MCG12213.002 Cumulative Probabilities for Colloid Transport Between One Matrix and Another Calculated from Interpolation of Pore Volume Data from Yucca Mountain Hydrologic (Stratigraphic) Samples LA0002PR831231.003 Probabilities for constants and retardation factors from Cwells microsphere data. ANL-NBS-HS-000026, Rev 00 14 April 2000 4.2 CRITERIA This AMR complies with the DOE interim guidance (Dyer 1999). Subparts of the interim guidance that apply to this analysis or modeling activity are those pertaining to the characterization of the Yucca Mountain site (Subpart B, Section 15), the compilation of information regarding hydrology of the site in support of the License Application (Subpart B, Section 21(c)(1)(ii)), and the definition of hydrologic parameters and conceptual models used in performance assessment (Subpart E, Section 114(a)). 5. ASSUMPTIONS 5.1 GENERAL ASSUMPTIONS OF THE PARTICLE-TRACKING TECHNIQUE In this section, the assumptions and general approach taken to develop the UZ radionuclide transport model are outlined as the first step toward developing the computational and mathematical models needed in Performance Assessment (PA) calculations. Prediction of solute transport is a critical element of many groundwater flow studies, including contaminant transport and the movement of natural isotopes or dissolved ions in solution. Modeling efforts typically are motivated by the need to predict the movement of a pollutant or dissolved chemical in the subsurface to answer practical questions concerning the rate and direction of contaminant movement and the predicted concentration in solution. In a typical solute transport simulation, a dissolved chemical is introduced into a steady-state or time-varying flow field, and the fate of the chemical is tracked while undergoing physical and chemical processes such as advection, dispersion, chemical and biological reaction, or diffusion into deadend pore space. Often, a concentration front is established that must be tracked accurately. In addition, many field investigations employ natural or introduced tracers to study the flow and transport system. These studies also require models to simulate the movement of dissolved species. Traditional solutions to the advection-dispersion (CD) equation, such as those used in most finite-element or finite-difference flow and transport codes, are versatile and allow the simultaneous solution of multiple interacting species. One drawback of a finite-difference or finite-element solution to the CD equation is that significant numerical dispersion may arise in the portion of a computational domain occupied by a front of rapidly varying concentration. Reducing the numerical dispersion requires either increased grid resolution or higher-order approximation methods, both of which may lead to prohibitive computational costs. Numerical dispersion is identical in character to actual dispersion, so it is difficult to separate numerical from actual dispersion in complex transport simulations. Approaches to cope with this problem include front-tracking algorithms with multiple grids (e.g., Yeh 1990, Wolfsberg and Freyberg 1994), the method of characteristics (e.g., Chiang et al. 1989), hybrid Eulerian-Lagrangian solution techniques (Neuman 1984), and particle-tracking techniques (e.g., Tompson and Gelhar 1990). Front-tracking algorithms solve the CD equation in integrated form on a numerical grid while tailoring the mesh to increase the resolution of the calculation at fronts. In contrast, an Eulerian-Lagrangian technique casts the CD equation using ANL-NBS-HS-000026, Rev 00 15 April 2000 the total derivative, so that the advection portion of transport can be solved accurately using particle-tracking techniques or the method of characteristics, while the dispersion component of transport is solved on a finite-difference or finite-element grid using standard techniques. Particle-tracking transport models take a fundamentally different approach. The trajectory of individual molecules or packets of fluid containing molecules are tracked through the model domain. When the fluid path lines are the model result of interest (Pollack, 1988; Lu, 1994), a relatively small number of particles can be used to trace the streamlines. Particle tracking is also used to simulate solute transport, such as the migration of a contaminant plume (Akindunni et al. 1995) or the prediction of breakthrough curves in interwell tracer experiments (Johnson et al. 1994). For these applications, a relatively large number of such particles must be used to obtain accurate solutions to the transport problem. Particle tracking has also been used to solve the advective portion of complex reactive transport models that simulate chemical reactions among multiple species (Fabriol et al. 1993). In a typical particle-tracking algorithm, a particle is sent to a new position assuming that the magnitude and direction of the velocity vector are constant during a time step . If small enough time steps are taken, particle pathways can be tracked accurately. Dispersion is treated as a random process that diverts the particle a random distance from its dispersion-free, deterministic path. In these so-called “random walk” models (e.g., Kinzelbach 1988), dispersion is usually calculated stochastically subject to a Gaussian model to reproduce the specified dispersion coefficient. The technique has also been extended by employing non-Gaussian random walk functions to represent scale-dependent dispersion (Scheibe and Cole 1994). Linear equilibrium sorption can be handled through the use of a retardation factor to correct the magnitude of the particle velocity. A crucial component of most random-walk particle-tracking algorithms developed to date is the need to accurately estimate the velocity at every position in the model domain. In the context of a finite-difference or finite-element numerical code, this means that velocities at positions other than the node points of the fluid flow grid must be computed using an interpolation scheme. Many studies have proposed and studied the accuracy of different interpolation schemes, including methods developed for regular, two- or three-dimensional finite difference grids (Schafer-Perini and Wilson 1991, Zheng 1993), for two- and three-dimensional finite-element grids (Cordes and Kinzelbach 1992), and for codes that employ the boundary element method for computing fluid flow (Latinopoulos and Katsifarakis 1991). Special techniques have been developed to handle complexities such as point fluid sources and sinks and abrupt changes in the conductivity of the medium (Zheng 1994). Unfortunately, many of the velocity interpolation schemes used in conventional particle-tracking techniques are computationally intensive, thus limiting the number of particles that can practically be used. Another drawback to traditional particle-tracking approaches is that spatial and temporal discretization often results in numerical inaccuracy in the fluid flow solution upon which velocity determinations are based. Thus, precise and time-consuming velocity interpolation schemes may not be justified in finite-difference or finite-element models. Finally, and most important for the simulation of transport in the UZ at Yucca Mountain, dualpermeability models employ overlapping continua to represent fracture and matrix flow ANL-NBS-HS-000026, Rev 00 16 April 2000 (Zyvoloski et al. 1992, Zimmerman et al. 1993). To develop a streamline-based particle-tracking method for dual-permeability models, velocity interpolations on each continuum would have to be coupled to a transfer term that allows particles to move from one medium to the other. This additional complexity, along with the inherent approximations associated with the dualpermeability method itself, may make precise velocity interpolation calculations of limited validity. In this AMR, a new particle-tracking technique is developed for transient, multi-dimensional finite-difference or finite-element codes. The algorithm is designed for computing solute concentration fields quickly and easily with structured or unstructured numerical grids of arbitrary complexity. Both continuum and dual-permeability formulations can be simulated. This flexibility is accomplished by extending the cell-based strategy of Desbarats (1990) for mapping out the path of the particle. In this method, the calculation of an “exact” pathline is replaced with a cell-to-cell migration of the particle. The mass flux from cell to cell is used directly, and no velocity interpolations are required. Since numerical solutions for fluid flow are typically mass-conservative (though not necessarily accurate) the particle-tracking method automatically conserves mass. In subsequent sections, the mathematical basis for this algorithm is outlined, and theory is developed to incorporate the effects of sorption, dispersion, and matrix diffusion into this new particle-tracking framework. Then, the technique is verified by comparing the results to a variety of analytical solutions and by benchmarking the code against a finite-element solution of the advection-dispersion for a one-dimensional, dual-permeability flow field. Finally, specific capabilities are developed for the radionuclide transport simulations, including source term methods, decay chains/ingrowth options, and colloid-transport methods. 5.2 ASSUMPTIONS OF THE ACTIVE FRACTURE MODEL FOR DETERMINING FRACTURE SPACING AND APERTURE All assumptions stated in this section are used in Section 6.2.1. 1. Fracture frequency, aperture, and permeability are log-normally distributed. Basis: These properties are a quantities bounded at the low end by zero. Therefore, a lognormal distribution is a natural choice that can meet the measured means and standard deviations, and are constrained to be larger than zero. This assumption does not require verification. 2. The cubic law is a valid approximation for gas permeability in fractured rock at Yucca Mountain. Basis: The cubic law for fracture permeability is a basic relationship used to relate fracture properties and fracture permeability (Bear et al. 1993, p. 15). This is an adequate approximation, based on scientific judgement, for the purposes of establishing the fracture aperture variance. No verification is required. ANL-NBS-HS-000026, Rev 00 17 April 2000 3. Active fracture model appropriately accounts for reduced fracture/matrix interaction. Basis: The reduction in fracture/matrix contact area is a result of the active fracture unsaturated flow model. This reduction is justified on the basis of the desirability of maintaining consistency with the assumptions underlying the development of the flow fields developed for the UZ flow modeling effort. These assumptions are developed in Liu et al. (1998). No further confirmation of this assumption is required. 5.3 ASSUMPTIONS OF THE COLLOID TRANSPORT MODEL All assumptions stated in this section are used in Section 6.1.4. 1. Radionuclide mass sorbs reversibly to non-diffusing colloids with a constant equilibrium sorption parameter fluid coll c C C K / = , where coll C is the radionuclide concentration residing on the colloids (moles radionuclide on colloid per kg fluid), and fluid C is the corresponding concentration in the fluid phase (moles aqueous radionuclide per kg fluid). Basis: Most measurements of sorption onto bulk rock and colloids are interpreted using an equilibrium sorption model such as this one. For compatibility with the data collected on sorption, this assumption is adopted in the numerical model as well. If attachment to colloidal particles is thought to be truly irreversible, this model can be used with an extremely large value of Kc for that portion of the radionuclide inventory. This assumption does not require verification. 2. Colloids undergo reversible filtration in the porous medium, with a colloid retardation factor of coll R . Basis: To estimate retardation of colloids in the fracture continuum, field experiments at the Cwells complex near Yucca Mountain were examined, in which transport of microspheres was used as an analog for colloids. The microsphere breakthrough curves were fit to forward and reverse filtration rates (DTN: LA0002PR831231.003). These rate constants were then used to calculate a retardation factor for colloid transport through saturated fractured rock (Table 7 of this document, from CRWMS M&O 2000b). For compatibility with this analysis of field data, this assumption is adopted in the numerical model as well. This assumption does not require verification. 3. Colloids undergo dispersion with identical dispersivity as an aqueous solute. Basis: When dispersivity is used to model solute spreading in porous media, it is introduced to capture variability in the flow velocity that exists at smaller scales than are modeled in the numerical grid. To a first approximation, this variability will act similarly on aqueous and colloidal components. Therefore, the same dispersivity should be used for both. This assumption does not require verification. ANL-NBS-HS-000026, Rev 00 18 April 2000 6. ANALYSIS/MODEL 6.1 THE RTTF PARTICLE-TRACKING TECHNIQUE This section outlines the development of the general transport methods used for the RTTF particle-tracking model, and Section 6.2 discusses issues specific to the use of this model to simulate radionuclide transport for the Yucca Mountain UZ. 6.1.1 Basic Methods The particle-tracking method developed in the present study views the fluid flow computational domain as an interconnected network of fluid storage volumes. Particles travel only from cell to cell, requiring no greater resolution of the particle pathways. In this sense, the method is similar to the node-to-node routing method of Desbarats (1990, p. 156). This simple starting point has been extended to include many different transport submodels and complex flow domains. The description that follows is applicable for steady-state, single-porosity flow fields; the corrections to the method for treating transient flow systems and dual-porosity model formulations are discussed in subsequent sections. The two steps in the particle-tracking approach are: 1) determine the time a particle spends in a given cell; and 2) determine which cell the particle travels to next. These two steps are detailed below. The residence time for a particle in a cell is governed by a transfer function describing the probability of the particle spending a given length of time in the cell. Thus, this particle-tracking approach is called the “residence time/transfer function” (RTTF) method. The schematic plots shown in Figures 1 and 2 illustrate the basis of the RTTF approach. For a cumulative probability distribution function of particle residence times, the residence time of a particle in a cell is computed by generating a random number between 0 and 1 to determine the corresponding residence time from the distribution function. In this example, the advection-dispersion equation was used to generate the RTTF curve, but other transport mechanisms can be incorporated as well, as demonstrated below. Figure 1. Schematic of the Cell-Based Particle-tracking Technique ANL-NBS-HS-000026, Rev 00 19 April 2000 Figure 2. Schematic of the RTTF Technique for Determining Particle Residence Time in a Cell If a large number of particles pass through the cell, the cumulative residence time distribution (RTD) of particles in the cell will be reproduced. Particle-tracking models of single-fracture transport (Yamashita and Kimura, 1990) have employed this approach to simulate fracture transport of diffusion into the rock matrix. From the solution of the flow field in a numerical model, the mass of fluid in the computational cell and the mass flow rate to or from each adjacent cell is computed. In the simplest case, the residence time of a particle in a cell, part t , is given by . = = out f f part m M & t t (Eq. 1) where f M is the fluid mass in the cell and the summation term in the denominator refers to the outlet fluid mass flow rates from the cell to adjacent cells. In the absence of dispersion or other transport mechanisms, the transfer function describing the distribution of particle residence times is a Heaviside function (unit step function) that is unity at the fluid residence time f t , because for this simple case, all particles entering the cell will possess this residence time. Equilibrium, linear sorption is included by correcting the particle residence time by a retardation factor f R . Thus, f f part R t t = , and f R is given by f d b f S K R f . + =1 (Eq. 2) where ANL-NBS-HS-000026, Rev 00 20 April 2000 d K is the equilibrium sorption coefficient (mL fluid/g rock) b . is the bulk rock density (g rock/mL total) f is the porosity(mL pore space/mL total) f S is the saturation of the phase in which the particle is traveling (mL fluid/mL pore space). Once again, in the absence of other transport processes, the transfer function is a Heaviside function. Before discussing more complex transfer functions for the RTTF method, the method for determining which cell a particle travels to after completing its stay at a given cell is discussed. The approach that is consistent with the RTTF method is that the probability of traveling to a neighboring cell is proportional to the mass flow rate to that cell. Only outflows are included in this calculation; the probability of traveling to an adjacent node is 0 if fluid flows from that node to the current node. A uniform random number from zero to one is used to make the decision of which node to travel to. Thus, the particle-tracking algorithm is: 1) compute the residence time of a particle at a cell using the RTTF method; and 2) at the end of its stay, send the particle to an adjacent cell randomly, with the probability of traveling to a given cell proportional to the mass flow rate to that cell. 6.1.2 Dispersion Transport processes such as dispersion can be incorporated into the RTTF particle-tracking algorithm through the use of transfer functions. For dispersion, within a computational cell, the equation for one-dimensional, axial dispersion is applied. The transport equation and boundary conditions for the one-dimensional, advective-dispersion equation are: dz C dz C D dt C R eff f . - . = . . 2 2 (Eq. 3), 0 = C , 0 = t (Eq. 4), 0 C C = , 0 = z (Eq. 5), 0 = C , 8 . z (Eq. 6), where C is the concentration (moles/kg fluid) 0 C is the injection concentration (moles/kg fluid) . is the flow velocity (m/s) eff D is the effective dispersion coefficient (m2/s), given by v Deff a = , where a is the dispersivity of the medium (m). ANL-NBS-HS-000026, Rev 00 21 April 2000 Here the molecular diffusion coefficient is ignored, since in general it is much smaller than the flow dispersion component of eff D . A non-dimensional version of Equation 3 can be obtained by the following transformations: 0 /C C C = ) , L z z / = ) , and L R t f / . . = , where L is the flow path length. The solution to Equations 3-6 is obtained after manipulation of Freeze and Cherry (1979, p. 391, Equation 9.5), yielding: ( ) ( ) . .. . . .. . . .. . . .. . + + . .. . . .. . - = . . . . 2 1 2 1 2 1 Pe erfc e Pe erfc C Pe ) (Eq. 7) where Pe is the Peclet number (dimensionless), a . / / L D L Pe eff = = . The use of this solution in the RTTF particle-tracking method requires that the transport problem be advection-dominated, such that during the time spent in a computational cell, solute would not tend to spread a significant distance away from that cell. Then, the approximate use of a distribution of times within a single cell will be adequate. Quantitatively, the criterion for applicability is based on the grid Peclet number a / x Peg . = , where x . is the characteristic length scale of the computational cell. Note that in contrast to conventional numerical solutions of the advective-dispersion equation, coarse spatial discretization is helpful for satisfying this criterion. Of course, the mesh spacing must still be small enough to provide an accurate flow solution. Highly dispersive transport invalidates the assumptions of the RTTF particle-tracking technique. When dispersion coefficients are large, accurate solutions to the advective-dispersion equation are easily obtained by conventional finite-difference or finite-element techniques, so these techniques should be used instead under these circumstances. For multi-dimensional flow systems, the dispersion model developed for one-dimensional systems can be extended to include dispersion coefficient values aligned with the coordinate axes. For this case, the flow direction is determined by the vector drawn from the nodal position of the previous cell to the current cell, and the dispersivity for this flow direction is computed from the equation for an ellipsoid: 2 2 2 2 2 2 / / / z y x z y x L a a a a . + . + . = (Eq. 8) The RTTF particle-tracking technique cannot be simply formulated with a longitudinal and transverse dispersion coefficient model, with the tensor aligned with the flow direction, because the flow rates between cells are defined rather than the actual flow velocity at a position. For a dispersion model aligned with the flow direction, a random-walk particle-tracking method such as that of Tompson and Gelhar (1990), also implemented in the saturated zone particle-tracking algorithm of FEHM, or a conventional finite-element or finite-difference solution to the CD equation, such as the reactive transport solution module in FEHM, should be used instead. The numerical implementation of this technique requires the determination of the dimensionless time . in Equation 7 for a randomly determined value of the dimensionless concentration C ) . ANL-NBS-HS-000026, Rev 00 22 April 2000 This determination is accomplished numerically in the particle-tracking code by fitting Equation 7 at selected values of . between 1 and 1000 using a piecewise continuous series of straight lines spanning the entire range of values. Then, the value of . at an arbitrary value of Pe is computed by linear interpolation between values determined at the Peclet numbers that bracket the actual value. This technique, involving a simple search for the correct type curves, followed by the calculation of two values of . and an interpolation, is much more computationally efficient (about a factor of two in cpu time) and robust than an iterative approach to the exact solution using Newton’s method. Solutions of adequate accuracy (less than 1% error) for Peclet numbers between 1 and 1000 are easily obtained using this linear-interpolation method. 6.1.3 Matrix Diffusion Matrix diffusion has been recognized as an important transport mechanism in fractured porous media (e.g., Neretnieks 1980; Robinson 1994). For many hydrologic flow systems, fluid flow is dominated by fractures because of the orders-of-magnitude larger permeabilities in the fractures compared to the surrounding rock matrix. However, even when fluid in the matrix is completely stagnant, solute can migrate into the matrix via molecular diffusion, resulting in a physical retardation of solute compared to pure fracture transport. This effect has recently been demonstrated on the laboratory scale by Reimus (1995) and on the field scale both by Maloszewski and Zuber (1991) and in the saturated zone at Yucca Mountain by Reimus et al. (1999). To derive a transfer function for matrix diffusion, an idealized representation of the transport system must first be developed. In this particle-tracking algorithm, the model depicted in Figure 3 is used to provide a transfer function for the case of fracture flow and diffusion between equally spaced fractures. In this model, additional terms not defined earlier are as follows: t is the time (s) 2b is the fracture aperture (m) 2B is the mean fracture spacing (m) f is the porosity of the matrix Rf,m is the retardation factor in the matrix Rf,f is the retardation factor in the fracture D is the effective diffusion coefficient (m2/s) ANL-NBS-HS-000026, Rev 00 23 April 2000 f D Rf,m, Rf,f Figure 3. Schematic of the Matrix-Diffusion Submodel Transport in the fractures is governed by Equation 3 with an additional term describing diffusion from the fracture to the matrix: b q dz C dz C D dt C R eff f f - . - . = . . 2 2 , (Eq. 9) where the flux at the fracture/matrix interface is given by b x dx C D q = . - = f (Eq. 10) and transport within the matrix is described by the one-dimensional diffusion equation 2 2 , x C D t C R m f . . = . . (Eq. 11), The molecular diffusion coefficient is the product of the free diffusion coefficient of the solute in water and a tortuosity factor to account for the details of diffusion through a tortuous, fluid-filled pore network. In this model, D is treated as the fundamental transport parameter, recognizing that it is a property of both the solute and the medium. Although a particular boundary condition in between the fractures is depicted in Figure 3, there are a variety of analytical solutions to this transport problem depending on the nature of the boundary condition in the x-direction. An analytical solution is given by Tang et al. (1981) for the semi-infinite boundary condition 0 / = . . x C as B x . . For the case of plug flow (no ANL-NBS-HS-000026, Rev 00 24 April 2000 dispersion) in the fractures, the solution of Starr et al. (1985, Equation 8b) can be used. Replacing Starr’s term s v x / with f t , as is called for in the RTTF technique, and ignoring radioactive decay, their solution reduces to . .. . . .. . - = f f m f f R t b D R erfc C C t ft 2 , 0 (Eq. 12) for f f R t t > , and 0 / 0 = C C for f f R t t = . The semi-infinite boundary condition between fractures limits the validity of either of these solutions to situations in which the characteristic diffusion distance for the transport problem is small compared to the fracture spacing B . As long as the solute has insufficient time to diffuse to the centerline between fractures, the solution provided by Tang et al. (1981) can be used as the transfer function for the particle-tracking technique. Alternatively, the boundary condition depicted in Figure 3 can be used to provide a transfer function for the case of fracture flow and diffusion between equally spaced fractures. Under these conditions, the analytical solution provided by Sudicky and Frind (1982) can be used to obtain the transfer function. The derivation of a form of this solution suitable for incorporation into the particle-tracking methodology has been presented in CRWMS M&O (2000c) for the saturated zone (SZ) particle-tracking transport model. The same subroutines in FEHM are used for both the UZ and SZ transport models. Although in principle a solution such as Tang et al. (1981) that includes dispersion and matrix diffusion could be used directly for the transfer function, its complex form makes it very inconvenient for rapidly computing particle residence times. Instead, a two-step process is used wherein the residence time in the fracture is first computed using the transfer function for onedimensional dispersion (Equation 7) without sorption. This fracture residence time is then used in the plug-flow equation with matrix diffusion and sorption to compute the particle residence time. To use Equation 12 as a transfer function, a numerical algorithm was developed to determine the inverse of the error function, that is, the value of d x for a given value of d y , such that ( ) d d x erf y = (note that ( ) ( ) d d x erf x erfc - = 1 ). The numerical implementation of this method entails dividing the error function into piecewise continuous segments from which the value of d x is determined by interpolation. The use of the two-step approach is justified because of the principle of superposition, which allows the dispersive process in the fracture to be decoupled from diffusive transport in the matrix. Proof that this numerical technique is acceptable is presented in the Code Verification section (Section 6.3) of this document. For the finite fracture spacing model, there are two options for defining the fracture spacing. The first is a simple node-by-node assignment of the fracture spacing. The second, the so-called “active fracture model,” is more consistent with the UZ flow model and, thus, is outlined in detail in Section 6.2.1. ANL-NBS-HS-000026, Rev 00 25 April 2000 6.1.4 Colloid Transport For colloid-facilitated transport, the transport equations for matrix diffusion, with either the semi-infinite or finite fracture spacings, can be simply revised. Given the assumptions listed in Section 5.3, the expression for transport for contaminant on colloids is: dz C dz C D dt C R coll coll eff coll coll . - . = . . 2 2 (Eq. 13) where eff D is the same as for an aqueous solute (assumption 3, Section 5.3). Combining Equations 13 and 9, and making use of the relation fluid coll c C C K / = : ) 1 ( 1 2 2 c eff c coll c f K b q dz C dz C D dt C K R K R + - . - . = . . .. . . .. . + + . (Eq. 14) This derivation implies that the transport equation for matrix diffusion can be revised to include colloid facilitated transport by replacing the half-aperture b by ( )b K b c coll + = 1 (Eq. 15) And the retardation factor in the fracture by c coll c f coll f K R K R R + + = 1 , (Eq. 16) Alternatively, for transport in a porous continuum, the solution to Equation 3 is used, with the retardation factor given by Equation 16, and f R replaced by m f R , . These relationships are built into the FEHM particle-tracking code, so that the additional terms c K and coll R are provided as inputs. In addition to the transport of radionuclides bound to colloids, there are several mechanisms related to the migration of the colloids themselves that can be simulated in the model. Above, the reversible retardation factor for colloid migration coll R was introduced. Due to the colloid size and surface properties relative to the pore size and surface charge, colloids can also undergo size exclusion and/or filtration in porous media. In the particle-tracking module, models have been implemented for these processes. For advective flow from fracture to matrix in the dualpermeability model, a size-exclusion model is implemented whereby colloids can remain in the fracture in proportions greater than the relative flow rate entering the matrix. A size exclusion parameter 1 = coll f is defined such that the probability of particles entering the matrix due to advective transport is multiplied by this factor. Therefore, complete exclusion from the matrix is obtained by setting 0 = coll f , whereas aqueous solute behavior is retained by setting 1 = coll f . Filtration, resulting in complete immobilization of the particle, can also be simulated at specified ANL-NBS-HS-000026, Rev 00 26 April 2000 interfaces within either the fracture or matrix continua. To invoke this mechanism, a filtration factor filt f at an interface (the finite element connections between two specified zones in an FEHM simulation) is defined. If a particle is slated to pass from one zone to the other via advective transport, filt f is the probability the particle continues moving, ( filt f - 1 is the probability that it is irreversibly removed by filtration). When using the filtration option, a word of caution is warranted. Colloid simulations are typically used to provide a mechanism for radionuclides to travel in the water bound to colloids. Filtration renders the colloids immobile, which in reality, only renders the radionuclide immobile if it is irreversibly bound to the colloid. When the radionuclide is only weakly sorbed to the colloid, the filtration option will artificially remove radionuclide mass from the system, resulting in a non-conservative simulation. Therefore, the filtration option should only be invoked for irreversibly bound radionuclides or when simulating colloid tracer experiments. The reversible model, using coll R to delay the migration, should be used instead for colloid-facilitated transport of radionuclides. 6.1.5 Particle Sources and Sinks There are two methods for introducing particles into the flow system: (1) the particles are either injected with the source fluid entering the model domain or (2) released at a particular cell or set of cells. The first method is used to track source fluid as it passes through the system. The number of particles entering with the source fluid at each cell is proportional to the source flow rate at that cell, which is equivalent to injecting fluid with a constant solute concentration. For Method 2, an arbitrary number of particles are released at each specified cell, regardless of the source flow rate. In the present application, Method 2 is used to input particles, which are used to represent radionuclide mass into the system at the repository level. Within Method 2, there are various ways to input a time-varying source of particles. For standalone simulations, the particles are inserted at a constant rate for a specified duration. There is also an option, used when the FEHM code is interfaced with GoldSim, to input a time-varying and spatially varying source mass flux into the model. The details of the method for accepting complex sources of multiple radionuclides from the EBS model are discussed in a subsequent section. When fluid exits the model domain at a sink, the model treats this flow as another outlet flow from the cell. The decision of whether the particle leaves the system or travels to an adjacent cell is then made on a probabilistic basis, just as though the fluid sink were another connected cell. Thus, the complexities discussed by Zheng (1994) for handling a so-called weak sink are avoided in the RTTF particle-tracking model. 6.1.6 Transient Fluid Flow When the RTTF particle-tracking method is implemented for a time-varying fluid flow system, the approach is somewhat more complex but still tractable. Consider a numerical simulation in which a discrete time step is taken at time t and a new fluid flow field is computed. In this ANL-NBS-HS-000026, Rev 00 27 April 2000 model, the new fluid flow time t t tnew . + = is treated as an intermediate time at which the particle-tracking calculation must stop. The time is intermediate because if the flow field is at steady state, there is no reason to stop at any time before the end of the simulation except to record particle information for output or processing purposes. The fate of all particles is tracked from time t to time new t assuming that the flow field is constant over this time interval. When the simulation reaches new t , the position (cell number) of the particle is recorded, along with its fractional time remaining at the cell and the randomly generated y-coordinate of the transfer function used for that particle in the cell. When the new fluid flow solution is established, the remaining residence time for a particle is determined from the following steps: 1. Compute a new fluid residence time f t . 2. Using the y-coordinate of the transfer function previously computed and the new transfer function, calculate a new particle residence time. 3. Multiply this time by the fractional time remaining in the cell to obtain the remaining time in the cell. This method approximates the behavior in a transient system, while reducing to the behavior that would be obtained in an unchanging flow field had the calculation not been forced to stop at the intermediate time. Another transient effect that must be considered is that the sum of the outlet mass flow rates . out m& in Equation 1 does not necessarily equal the sum of the inlet mass flow rates. When there is net fluid flow into a cell, the particle-tracking algorithm uses the sum of the inlet flow rates in Equation 1, whereas Equation 1 itself is used when there is net outflow from a cell. 6.1.7 Dual-Permeability Formulation In the development of the dual-permeability flow model, the Yucca Mountain project recognized the need to explicitly account for fracture flow in some hydrologic units. A dual-permeability formulation allowed for a greater liklihood of fracture flow by relaxing the assumption inherent in the equivalent continuum model formulation that the fracture and matrix media must be in capillary pressure equilibrium at all locations. In the dual-permeability models, significant fracture flow occurs in units with low matrix permeability such as the Topopah Spring and zeolitic Calico Hills units. This advance allowed the project to reconcile the infiltration rates estimated in surface and soil moisture studies with the fact that in many units, matrix rocks were at less than complete saturation. In an equivalent continuum model formulation, this observation implied infiltration rates that were much lower than values estimated by more direct means. In the dual-permeability model, fractures are able to flow even though there is a high capillary pressure in the matrix blocks. Because of the success of dual-permeability models in their ability to simulate fracture flow, it would be tempting to use the same model formulation for radionuclide transport. Although this approach would have the advantage of simplicity and consistency with the flow model, there ANL-NBS-HS-000026, Rev 00 28 April 2000 may be significant limitations that would render a simple dual-permeability transport formulation less than adequate for radionuclide transport modeling. The restriction of a single matrix node for each fracture node means that sharp concentration gradients away from the fracture cannot be captured. For flow modeling, this may not be overly restrictive because of the generally more rapid migration of a pressure front into the matrix. Furthermore, the typical scenario of flow modeling is the simulation of fluid velocities, pressures, and saturations under the assumption of steady state flow. By contrast, transport of a contaminant plume is inherently a transient problem, and characteristic diffusional distances can be of the order of centimeters during unsaturated transport through a fractured medium. If the characteristic spacing between flowing fractures is much greater than centimeters, then the use of one grid point in the matrix is questionable for capturing the fracture/matrix diffusion term in the transport mass balance. Clearly, the biggest advantage of the dual-permeability model formulation for transport is the ability to capture, in a computationally efficient manner, the disparate flow velocities in the two media. However, because of the potential limitations of the dual-permeability formulation for transport, a hybrid model was developed using the cell-based particle-tracking algorithm and the RTTF approach. This method addresses the issue of potentially sharp gradients away from the fractures while relying on the dual-permeability flow formulation to adequately capture the distribution of advective velocities within the model domain. Extension of the advection term in the RTTF particle-tracking technique to handle dual-permeability systems is straightforward, because the fracture-matrix flow term is simply an additional inlet or outlet flow rate from the node. In the algorithm, this flow term is treated like flow to any other node, and the particle can shift from one continuum to the other. In addition, to simulate molecular diffusion as an additional fracture-matrix transport interaction term, the matrix diffusion option can be invoked for particles traveling in the fracture continuum. In this formulation, the simplification employed is that for the purposes of capturing the matrix diffusion component of the transport, particles delayed by matrix diffusion in the fracture continuum remain in the fracture continuum. Advective motion can drive particles from one continuum to the other, as often occurs in dualpermeability models with units of contrasting hydrologic properties. However, the matrix diffusion option assumes a model system in which the fracture and matrix flow systems are weakly coupled, such that matrix diffusion in a fracture can be simulating without explicitly accounting for the superimposed effect of advection. Of course, the method for modeling solute diffusion between the fractures and matrix is approximate, and needs to be tested against a model that more precisely captures the transport behavior. In Section 6.4, a series of dual-permeability benchmarking calculations are presented to explore the adequacy of the model for the purpose of simulating radionuclide migration in the Yucca Mountain unsaturated zone and to highlight potential strengths and limitations of the model. 6.1.8 Output Options There are several methods for obtaining output results for a particle-tracking model simulation. Generally, the results are reported in the form of a breakthrough curve, that is, a distribution of arrival times of particles at a fluid sink. These sinks can be lumped together, so that, for example, the breakthrough curve can be related to a mass flux of radionuclide at the water table in the UZ ANL-NBS-HS-000026, Rev 00 29 April 2000 transport model. The sink nodes can also be sub-divided into user-defined zones, and breakthrough curves can be obtained for each. In addition to breakthrough curves, in situ concentrations can be calculated and reported at a given node or nodes, or throughout the model domain. However, as is the case with all particletracking algorithms, concentration distributions can be subject to significant noise when used to compute concentration distributions. This problem can be especially acute for dual-permeability systems, for which transport velocities in the fractures can be so rapid that the probability of finding a statistically significant number of particles in a fracture grid node at a particular reporting period is very low. Alternatively, a provision in the code allows the cumulative number of particles passing through a cell to be reported. When a pulse input of particles is used, this type of concentration is equivalent to the response to a constant mass flux input of solute. Since all particles passing through a cell get counted, the problem just described is alleviated. 6.2 IMPLEMENTATION ISSUES FOR THE UZ TRANSPORT MODEL This section contains issues specifically related to the use of the particle-tracking model for the Yucca Mountain UZ. 6.2.1 Active Fracture Model The transport equation for a single fracture under unsaturated conditions can be obtained by a simple extension of Equation 9: b S q R dz C dz C D dt C R f fm A eff f f - . - . = . . 2 2 , (Eq. 17) Here Sf is the fracture saturation qfm is the fracture/matrix diffusive flux (moles/m2-s) RA is the fracture/matrix contact area reduction factor. The term Sfb may be thought of as the “water aperture” for the transport problem. In the active fracture model, only a portion of the fracture space participates in the flow and transport. The saturation of the active portion is the “active fracture saturation”. It is this saturation that applies to the transport equation because the inactive fractures can be thought of as unavailable, as if filled with minerals. Therefore, Sf is replaced by the active fracture saturation, Sa. The fracture/matrix area reduction factor is defined by Liu et al. (1998; Equation 13) to be: ae A S R = (Eq. 18) where Sae is the effective water saturation of the active fractures. The effective (or normalized) water saturation is given by Liu et al. (1998, Equation 4): ANL-NBS-HS-000026, Rev 00 30 April 2000 r r a ae S S S S - - = 1 (Eq. 19) where Sr is the residual water saturation. This area reduction factor was derived for unsaturated flow, for which it can be argued that the effective fracture/matrix flow process goes to zero at residual water saturation. For matrix diffusion and transport, however, the effective contact area should only go to zero when the physical saturation of the active fractures goes to zero. Therefore, a A S R = (Eq. 20) Substituting Sf = Sa and RA = Sa into the transport equation: b q dz C dz C D dt C R fm eff f f - . - . = . . 2 2 , (Eq. 21) which is the same equation as for transport in a saturated fracture. The adjustment of the fracture spacing to obtain the fracture spacing of active fractures is given by the following relationship (adapted from Liu et al., 1998, Equation 17): . - = S B B g ) (Eq. 22) where r r f S S S S - - = 1 ) (Eq. 23) g B is the geometric fracture spacing (including flowing and dry fractures) r S is the residual fluid saturation exponent . is a fitting parameter ( 1 0 < < . ). The model already has a local value of f S at every node, and thus requires g B , r S , and . as input parameters. In this section, the UZ parameter distributions to be used for the fracture aperture and spacing are derived for use in the finite fracture spacing matrix diffusion model. The inverse of the fracture aperture is used in the transport model as a measure of the amount of fracture surface area per unit volume of fracture pore space. Given that the fracture area per unit bulk volume, A, is available from the flow model, as listed in Table 3, and that the fracture pore volume per unit bulk volume (or fracture porosity, f f ) is also given in Table 3, it is a simple matter to compute the fracture aperture, b, as follows: ANL-NBS-HS-000026, Rev 00 31 April 2000 A b f f = (Eq. 24) The aperture derived from this formula will be used as the geometric mean of the distribution for aperture. The relationship for fracture spacing is simply the inverse of fracture frequency. Radionuclide transport is not expected to be very sensitive to fracture spacing. However, it is expected to be sensitive to fracture aperture. Therefore, a stochastic sampling for fracture aperture in the TSPA calculations is desired. Although the analysis above provides a simple means to evaluate the mean fracture aperture, the data for area and porosity do not lend themselves to the evaluation of a variance in fracture aperture for the purposes of sampling. However, using the data in Table 4, it is possible to establish a variance in fracture aperture by the following method. From the cubic law, the fracture aperture b 2 may be expressed as: 3 1 12 2 . .. . . .. . = f k b G (Eq. 25) where kG is the permeability (m2) of the fractured medium determined from estimates of gas permeability, and f is the fracture spatial frequency. Taking the logarithm of this expression gives: ) log( 3 1 ) log( 3 1 ) 12 log( 3 1 ) 2 log( f k b G - + = (Eq. 26) The mean and standard deviation for log(kG) and f are given in Table 4 for some of the model layers. Taking log(kG) and log(f) to be normally distributed with means mlogk and mlogf, respectively and variances, 2 log k s and 2 log f s , respectively, then the mean and variance for log(2b) are: f k b log log 2 log 3 1 3 1 ) 12 log( 3 1 µ µ µ - + = (Eq. 27) 2 log 2 log 2 2 log 9 1 9 1 f k b s s s + = (Eq. 28) ANL-NBS-HS-000026, Rev 00 32 April 2000 Table 3: Geometric Mean Fracture Aperture and Spacing Model Layer Fracture Frequency (1/m) Fracture Spacing (m) Fracture Porosity (dimension-less) Fracture/Matrix Connection Area (m2/m3) Mean Fracture Aperture (m) tcw11 0.92 1.09 2.8E-2 1.56 1.79E-02 tcw12 1.91 0.524 2.0E-2 13.39 1.49E-03 tcw13 2.79 0.358 1.5E-2 3.77 3.98E-03 ptn21 0.67 1.49 1.1E-2 1 1.10E-02 ptn22 0.46 2.17 1.2E-2 1.41 8.51E-03 ptn23 0.57 1.75 2.5E-3 1.75 1.43E-03 ptn24 0.46 2.17 1.2E-2 0.34 3.53E-02 ptn25 0.52 1.92 6.2E-3 1.09 5.69E-03 ptn26 0.97 1.03 3.6E-3 3.56 1.01E-03 tsw31 2.17 0.461 5.5E-3 3.86 1.42E-03 tsw32 1.12 0.893 9.5E-3 3.21 2.96E-03 tsw33 0.81 1.23 6.6E-3 4.44 1.49E-03 tsw34 4.32 0.231 1.0E-2 13.54 7.39E-04 tsw35 3.16 0.316 1.1E-2 9.68 1.14E-03 tsw36 4.02 0.249 1.5E-2 12.31 1.22E-03 tsw37 4.02 0.249 1.5E-2 12.31 1.22E-03 tsw38 4.36 0.229 1.2E-2 13.34 9.00E-04 tsw39 0.96 1.04 4.6E-3 2.95 1.56E-03 ch1z 0.04 25 1.7E-4 0.11 1.55E-03 ch1v 0.10 10 6.9E-4 0.3 2.30E-03 ch2v 0.14 7.14 8.9E-4 0.43 2.07E-03 ch3v 0.14 7.14 8.9E-4 0.43 2.07E-03 ch4v 0.14 7.14 8.9E-4 0.43 2.07E-03 ch5v 0.14 7.14 8.9E-4 0.43 2.07E-03 ch2z 0.14 7.14 4.3E-4 0.43 1.00E-03 ch3z 0.14 7.14 4.3E-4 0.43 1.00E-03 ch4z 0.14 7.14 4.3E-4 0.43 1.00E-03 ch5z 0.14 7.14 4.3E-4 0.43 1.00E-03 ch6 0.04 25 1.7E-4 0.11 1.55E-03 pp4 0.14 7.14 4.3E-4 0.43 1.00E-03 pp3 0.20 5 1.1E-3 0.61 1.80E-03 pp2 0.20 5 1.1E-3 0.61 1.80E-03 pp1 0.14 7.14 4.3E-4 0.43 1.00E-03 bf3 0.20 5 1.1E-3 0.61 1.80E-03 bf2 0.14 7.14 4.3E-4 0.43 1.00E-03 tr3 0.20 5 1.1E-3 0.61 1.80E-03 tr2 0.14 7.14 4.3E-4 0.43 1.00E-03 tcwf - faults 1.90 0.526 4.4E-2 1.3E+1 3.38E-03 ptnf - faults 0.54 1.85 1.6E-2 1.3E+0 1.23E-02 tswf - faults 1.70 0.588 3.6E-2 8.6E+0 4.19E-03 chnf – faults & below 0.13 7.69 1.6E-3 4.7E-1 3.40E-03 DTN: LB990501233129.001 ANL-NBS-HS-000026, Rev 00 33 April 2000 In fact, the only expression needed is the variance because the values in Table 4 will be used for the geometric mean fracture aperture. To use the relationship for 2 2 log b s , the relationship between the arithmetic mean and standard deviation for f, given in Table 4, to the variance for log (f) is desired. The relationship between 2 f s , 2 f µ , and 2 log f s is the following (Iman and Shortencarier 1984, p. 17): 1 ) exp( 2 log 2 2 2 - = f f f A s µ s (Eq. 29) where A= ln(10) is used to convert between logarithms of base 10 to base e. Solving for 2 log f s gives . .. . . .. . + . .. . . .. . = 2 2 2 2 log 1 ln ) 10 ln( 1 f f f µ s s (Eq. 30) The results are summarized in Table 4. Table 4. Aperture Data and Derived Geometric Standard Deviation of Fracture Aperture Model layer mlogk k log s f µ f s 2 log f s 2 log k s 2 2 log b s Geometric Std. Dev. of 2b tcw12 -11.279 0.778 1.91 2.09 0.14848 0.60528 0.08375 1.94715 tcw13 -11.344 1.147 2.79 1.43 0.04399 1.31560 0.15106 2.44722 ptn21 -11.491 0.885 0.67 0.92 0.19987 0.78322 0.10923 2.14044 ptn25 -12.784 0.101 0.52 0.6 0.15965 0.01020 0.01887 1.37207 tsw32 -12.146 0.658 1.12 1.09 0.12568 0.43296 0.06207 1.77477 tsw33 -12.112 0.612 0.81 1.03 0.18144 0.37454 0.06177 1.77235 tsw34 -12.474 0.546 4.32 3.42 0.09177 0.29811 0.04332 1.61486 DTN: LB990501233129.001 (aperture data) The average geometric standard deviation for these model layers is 1.9. This average geometric standard deviation is to be used for all the model layers beneath the potential repository to perform stochastic sampling of the fracture aperture in the transport calculations for TSPA. 6.2.2 Multiple Radionuclides With Decay/Ingrowth The FEHM (Zyvoloski et al. 1997) code allows particles to decay with or without the production of the daughter product. For multiple species with decay chains, the code uses the approach outlined below to determine the number of decayed particles, and the code performs the bookkeeping needed to keep track of the locations and numbers of each type of radionuclide. These multiple species can each have their own transport parameters such as sorption and diffusion coefficients. ANL-NBS-HS-000026, Rev 00 34 April 2000 For decay-ingrowth simulations with time-dependent release of tracer particles, the computational burden increases dramatically with the number of particles in the field. For example, the decay-ingrowth calculation for species i decays into species j at a decay rate . is: ( ) [ ] { } .= - - - = i N m m j t t N 1 exp 1 . (Eq. 31) where Nj is the number of particles of species j decayed from species i, Ni is the number of particles of species i, and tm is the time at which the m th particle is injected into the system. If 500,000 particles of species i are injected into the system, then at each time step, the number of mathematical operations for ingrowth calculations alone are around 2.5 million. For a simulation time period of 1 million years, the typical calculation requires about 100 time steps. Therefore, the total number of operations for ingrowth calculations will reach 0.25 billion. Therefore, for site-scale simulations, the use of Equation 31 would be extremely inefficient. To reduce the computational burden in simulations, the decay-ingrowth calculation in Equation 31 is approximated with an integral by assuming that particles are injected into the system uniformly in time domain. Multiplying both sides of Equation 31 by .t, the average injection time interval between particles, and approximate Equation 31 with an integral: ( ) ( ) ( ) [ ] t N j . . . . . . . - - - • + - ˜ 2 1 2 1 exp exp 1 .t .t . t t (Eq. 32) where 1 1 t t - = t and Ni t t - = 2 t , 1 t is the time at which particle injection starts, Ni t is the time of the Nth injected particle, and t is the time at which the decay-ingrowth calculation is carried out. The use of Equation 32 reduces the number of operations within one time step from millions of operations to just 10, which greatly increases the speed of simulations. Validity of this approach is demonstrated in Section 6.3. The decay-ingrowth model flow chart as implemented in FEHM is summarized below in Figure 4. Inside FEHM, mass is passed to the particle-tracking module either as direct input or by the Graphical Simulation Environment software package (GoldSim, Golder 2000). When the input is as mass, conversion factors (# of particles/mole) are introduced to convert mass into number of particles for use in particle-tracking simulations. The conversion factors are stored in memory and are passed to the daughter species when parent particles decay into daughter particles. The conversion factors guarantee that mass is conserved during the decay-ingrowth process. At the end of each simulation, FEHM uses the conversion factors to convert number of particles back into mass, then passes the results back to RIP or writes them to an output file. The accuracy of the integral approach depends on the number of particles and their release history. In general, the use of a greater number of particles increases the accuracy. When the same number of particles is used in simulations, the one with a constant release rate will cause less error. If the release rate changes with time, the release period is divided into segments so ANL-NBS-HS-000026, Rev 00 35 April 2000 that within each segment (corresponds to each time step) the release rate can be treated as a constant. Additional bookkeeping is required to handle the case of multiple species. For each species, particles released during a time step are stored in the corresponding particle array and are labeled as one segment of the entire array. The particle starting and ending release times, release time interval, and particle indexing information are stored in memory for use in decay-ingrowth calculations. Thus, a particle array consists of many segments with each segment corresponding to particles released at the corresponding time step. At each time step, particles are uniformly distributed in the field in two ways: 1. When the number of injected particles, np, is larger than the number of releasing nodes, n, in the field, each releasing node is assigned m=int(np/n) particles. Then, the code loops through all release nodes and at each node, evenly releases m particles into the field during the current time step. The remaining particles, if np is not a multiple of n, are injected as in 2. below. 2. When the number of particles, np, is smaller then the number of releasing nodes, n, a uniform random number generator is used to randomly select np nodes from a poll of n releasing nodes to release particles during the current time step. By doing so, particles are uniformly distributed over the repository region. In the FEHM decay-ingrowth model, a first in, first out approach is used to select which particles undergo decay. Alternatively, the approach could have been to first calculate the number of decay particles, nd, for each segment, then select decay candidates from a group of particles which are sequentially released during the corresponding particle injection period. Then, for the case in which several particles are injected simultaneously, to make sure that the selection of decay particles is uniform in space and representative, a uniform random number generator is used to select nd decay particles from the particle array. But this approach would require that at each time step, the selecting process be run to pick out decay particles, which would increase computational burden. Instead, an alternative approach was developed in which particles are first distributed uniformly in the field, after which they are sorted in ascending order in time domain based on their release time. Finally, the array index of particles having the same release time are randomly changed. By doing this, the sequence of particles with the same release time but at different release nodes are mixed in the array. Once the number of decay particles, nd, is calculated, the code transforms the first nd particles in the parent particle array into daughter species particles. This approach skips the process of selecting decay particles within each segment at each time step, and hence increases the efficiency of the decay model. ANL-NBS-HS-000026, Rev 00 36 April 2000 is a decayingrowth simulation? distribute particles uniformally in time and space over designated nodes no yes get mass from GoldSim or input data file no yes decay only decay-ingrowth decay only or decayingrowth? call decay-ingrowth subroutine carry out particle tracking return mass to GoldSim or print out results call decay subroutine convert mass into number of particles mass in # of particles? Figure 4. Flow Chart of FEHM Decay-Ingrowth Model For simulations involving decay without tracking of the daughter species, an indicator is set so that the code skips the decayed particle in the particle-tracking simulations. Figure 5 shows the implementation of decay simulations at each time step in FEHM. ANL-NBS-HS-000026, Rev 00 37 April 2000 no Does the current species decay? For each species no yes For each segment of the particle array, calculate number of particles decayed, nd Mark decayed particles in the segment and record the index of the last decayed particle. Done with decay calculations? Start particle tracking simulations yes Figure 5. Flow Chart of Decay Subroutine At each time step, for a decay species, the code loops through each segment and calculates the number of parent particles decayed, nd, within the current time step using the integration approach in Equation 32 as follows. Suppose the particle decay index which records the position of the last decayed particle during previous time step for the same segment is m. Then, at the current time step, the code marks nd particles sequentially in the particle array as decayed by pointing the decay index pointer at (m+nd) in the array of the parent species. The updated decay index pointers for each segment are stored in the memory for use in next time step. Once the decay calculations are performed, the particle-tracking module is started to simulate the movement of particles in the field. Using the decay index pointers stored in the memory, the code skips particles already decayed previously, thereby reducimg unnecessary operations and CPU time. Decay calculation for the decay-ingrowth model is the same as for the decay only model, except that after a particle decays, it is transformed into a daughter species particle. The flow chart for the decay-ingrowth simulations at each time step is shown in Figure 6. ANL-NBS-HS-000026, Rev 00 38 April 2000 Decay-ingrowth simulation? For each species no yes For each segment of the particle array, calculate number of particles decayed Mark decayed particles in the segment and record the index of the last decayed particle. no Done with the calculations? Start particle tracking simulations yes Convert the decayed particles into new particles of the daughter species. Figure 6. Flow Chart of the Decay-Ingrowth Subroutine A new particle generated due to decay-ingrowth takes over the location of the parent particle but possesses its own transport properties. This process requires that related particle information be passed from the parent particle to the daughter particle. This information includes: particle injection time, particle current residence time in the cell, and the corresponding particle conversion factor, etc. Newly generated daughter particles form a new segment in the daughter particle array. The corresponding starting and ending particle release times for the new segment is stored in the memory for use in sequential decay-ingrowth calculations. 6.2.3 Interface with TOUGH2 and GoldSim The UZ transport model can be run either as a standalone FEHM simulation or within a performance assessment simulation using the GoldSim application. Although FEHM can perform flow simulations as well, the UZ transport application utilizes only the particle-tracking portion of the code, and thus requires the flow field as input. Flow fields are being determined as part of the UZ PMR using the computer code TOUGH2 (Pruess 1991). The results are converted to FEHM input files using the computer code T2FEHM2 (CRWMS M&O 2000d). The information required to read the flow fields are the grid and connectivity information, the fluid saturation at each grid point, and the fluid mass flux between each node (and any fluid sources ANL-NBS-HS-000026, Rev 00 39 April 2000 and sinks). This conversion program produces input files that can be read by FEHM in lieu of FEHM computing the flow field. The interface between GoldSim and FEHM establishes a protocol for defining the radionuclide sources to the UZ transport model (provided by GoldSim), the definition of a particular flow field for FEHM to use, and exit mass fluxes of radionuclides from the UZ model (from FEHM to GoldSim based on the particle-tracking simulation). This protocol is quite flexible, allowing an arbitrary number of source regions, radionuclides, exit regions, and flow fields to be defined and passed between GoldSim and FEHM through the FEHM subroutine call statement. During a GoldSim simulation, FEHM cedes control of the time step to GoldSim. At each time, GoldSim provides a flag defining the flow field and the mass flux inputs of radionuclides. By changing the flow field during a simulation, the model can simulate the impact of a climate-induced change in the UZ system. When this occurs, FEHM reads in the new flow field and proceeds with the transport simulation. Since each flow field is a steady state flow field, the implicit assumption is a quasi-steady one, that is, the system establishes a new steady rapidly in comparison to transport velocities through the unsaturated zone. At the end of each time step, FEHM passes back to GoldSim the exiting mass flux values from the UZ model. 6.2.4 Engineered Barrier System (EBS) Random Release Model for Radionuclide Source Terms If waste packages containing high level radionuclide material in the repository eventually fail due to corrosion, the process will almost certainly be variable in both space and time. At early times, a few packages may fail, releasing radionuclides into the UZ. At later times, a greater number of packages may fail. In VA (Viability Assessment) calculations (DOE 1998, Figure 4- 10), releases from the EBS to UZ (Unsaturated Zone) were spread over the entire repository subregions. Such treatment of the EBS release could result in significant artificial dilution of the UZ transport source term in some circumstances. In reality, waste packages may not fail uniformly in space and time. Rather, a few waste packages may fail at early times, while others may fail gradually over longer time periods. An EBS random release model was developed in FEHM to allow the model to simulate early failed packages and time- and spatially variable radionuclide releases. To begin, a repository is defined consisting of N_large sub-regions. Each sub-region contains certain number of waste packages. Initially, M_fine packages fail at locations designated by package x-y locations (x,y). The M_fine failed packages release radionuclides at a mass flow rate of M_flux_i, where i is the ith failed package. As time proceeds, packages fail in the sub-regions of the repository. At each time step, there are a certain number of failed packages in each sub-region i. The mass flux released from those packages is denoted as N_large_flux for the ith sub-region. In this model, the release nodes in the numerical grid for the failed packages are randomly selected from the available repository nodes within that sub-region to mimic the failure process of the waste packages. The mass release of M_fine packages is separated from those of the other failed packages. ANL-NBS-HS-000026, Rev 00 40 April 2000 To simulate the impact of the EBS random release on system performance at the Yucca Mountain site, the FEHM EBS random release model was developed to perform the following tasks: • Locate the M_fine early failed package nodes in repository sub-regions based on given failed package coordinates. If no node matches a given coordinate, then select the nearest node to the given coordinate. • Randomly select the failed package nodes in the designated sub-region i. The existing FEHM subroutine getrip was modified to handle EBS random release. From FEHM particle-tracking subroutine part_track, subroutine getrip is called to determine the particle release locations. First, the subroutine obtains information passed by GoldSim in an input array called in[ ]. The structure of the in[ ] array is shown in Figure 7. time flow field index parameter index M_fine failed package (x,y)coordinates of early failed package (M_fine) N_large, # of repository sub-regions list # of failed packages in each sub-region (N_large) # of species mass input flag. 0: no mass input ; 1, there is mass input # of input buffers, M_fine+N_large # of output buffers mass release for each species during current time step. Values are passed for all species from the first releasing node to the M_fine th node, then, from the first sub-region to the N_large th sub-region. (M_fine+N_large)*number_of_species values. Figure 7. The Structure of the in[ ] Array Passed to FEHM from GoldSim The algorithm used in FEHM EBS random release model is summarized in Figure 8, the flow chart of the EBS random release model. Starting with the M_fine early failed packages, getrip extracts the (x,y) coordinates of the early failed packages and loops through each repository sub-region node to select the one that is closest to the given coordinates. To prevent a node being selected more than once for two or more given coordinates, getrip checks the selected nodes for overlapping. If overlapping is found, getrip prints out error messages to the error file fehmn.err, then stops the program. ANL-NBS-HS-000026, Rev 00 41 April 2000 no get the M_fine group coordinates and select nodes that have the shortest distance to the given coordinates randomly select N_newly_failed nodes in each subregion for radionuclide release in each sub-region, sum the mass flux of M_fine group with N_large group, then set number of particles released based on mass fraction of the M_fine group and N_large group, respectively is a node selected twice due to similar coordinate values yes no mass input = 0 From GoldSim? yes call subroutine setmptr to release particles from the selected nodes into the system Back to calling subroutine part_track stop, print out error message call getrip subroutine to get the M_fine and N_large group mass flux and node information passed by GoldSim within part_track subroutine Figure 8. Flow Chart of FEHM Random EBS Random-Release Model When the selection of M_fine nodes is complete, getrip extracts the number of failed packages in each sub-region for the N_large sub-regions. The number of failed packages at the current time step is compared with the values at the previous time step to determine the number of newly failed packages, N_newly_failed, within the current time step in each sub-region. Then, getrip randomly selects N_newly_failed nodes within the corresponding sub-region. The selected failed nodes are stored in the memory for use in releasing radionuclides and are removed from the base of available repository nodes. If the number of failed packages is larger than the number of nodes in a sub-region, then radionuclide will be released from all nodes within the sub-region. Once all nodes of failed packages are determined, getrip allocates the number of released particles proportionally to the mass flux values of each failed package. Then, subroutine setmptr is called to inject particles into the system for each species. To ensure that the code is functioning as designed, a verification example is included in a subsequent section. ANL-NBS-HS-000026, Rev 00 42 April 2000 6.2.5 Colloid Transport Parameters for the UZ Transport Model To simulate colloid size exclusion between the fracture and matrix and filtration processes at the interfaces between matrix unit, arguments are invoked related to the colloid size versus the pore size to determine conservative parameters. The pore size distribution for hydrogeologic units considered in the model can be estimated from moisture retention curves (DTN: GS980908312242.039, GS950608312231.008, GS960808312231.003). The pore sizes were calculated from Marshall et al. (1996, Equation 8.1a), which is . . .cos 2 - = r (Eq. 33) where r is the radius of pores, . is surface tension of the water, . is contact angle (assumed to be zero), and . is matric potential. This equation is valid for the range of pore sizes between 100 nm and 1000 µm. Many of the pore sizes in units such as the Topopah Springs welded (TSw4) are smaller. A correction for pore sizes less than 100 nm can be calculated if the relative vapor pressure of the soil water is known. This data was not available, so the correction was not made, and a slight error will exist in the pore sizes calculated less than 100 nm. To determine the percent of pores within each size range, the change in moisture content between measurements was used. The percentage of colloids that can enter a matrix unit from the fractures can be determined based on the percentage of the pores that were greater than the colloid size of interest. For example, for the TSw4, 71% of the pore sizes are estimated to be smaller than 100 nm (DTN: LA0003MCG12213.002). Therefore, for the size exclusion between the fracture and matrix only 29% of 100 nm colloids are allowed to enter the matrix (Table 5). The remainder of the colloids must therefore remain in the fracture continuum. The values used for the size exclusion are based on an expected colloid size of 100 nm, but values based on other colloid sizes could also be used. Colloid filtration at interfaces between matrix units is based on the colloid size distribution and the cumulative probability of colloids entering a unit based on the pore size distributions (Table 6). The implementation of distributions allows for a wide range of colloid sizes to be considered in a single simulation. For filtration at the interface between two units, the properties of the unit that flow is entering is used. The data in Table 6 show that the only interfaces where colloid filtration is likely to occur is between units like the vitric and zeolitic Calico Hills where there is a large contrast in pore size distribution. When a colloid can not enter a matrix unit at an interface, the colloid is removed from the simulation and is considered permanently filtered. As a result, the use of filtration at interfaces is not recommended for natural colloids associated with actinides because when the colloids are removed they no longer participate in reactions, such as desorption. These estimates are only based on the pore size distribution, discounting any effect of partial saturation. In the Tsw4 and the Ch1z, the average saturation is close to one, so any variations in water content would not significantly affect the values provided. For the CH1v, the average saturation is much lower. The use of the pore size distribution under fully saturated conditions is ANL-NBS-HS-000026, Rev 00 43 April 2000 justified by recognizing that the saturation is probably higher near the fracture matrix interface where the size exclusion rule is applied and lower further away from the fracture. Therefore, the matrix would probably be sufficiently saturated near the fracture that a reduction in the percent of colloids that can enter the matrix due to unsaturated conditions was not required. Table 5. Fraction of 100-nm Colloids That Can Enter the Matrix from the Fracture Based on Size Units TMN (TSW4) 0.29 TLL (TSW5) 0.39 TM2 (TSW6) 0.35 TMN1 (TSW7) 0.07 PV3 (TSW8) 0.10 PV2 (TSW9) 0.61 BT1a (CH1) 0.15 CHV 0.61 CHZ 0.27 BT (CH6) 0.08 PP4 0.02 PP3 0.79 PP2 0.35 PP1 0.43 BF3 0.26 BF2 0.04 DTN: LA0003MCG12213.002 To summarize the size exclusion and filtration information, our recommendation is to use colloid-size distribution information to set the colloid transport parameters based on the information in Tables 5 and 6. To estimate retardation of colloids in the fracture continuum, field experiments at the C-wells complex near Yucca Mountain were examined, in which transport of microspheres was used as an analog for colloids. The tracer experiments were conducted in both the conductive Bullfrog unit and less conductive Prow Pass unit. The microsphere breakthrough curves were then fit to forward and reverse filtration rates (DTN: LA0002PR831231.003). These were then used to calculate a retardation factor for colloid transport through saturated fractured rock (Table 7). Data on colloid filtration through partially saturated fractures is not available. Therefore, it is recommended that either: 1) the distribution of retardation factors for the saturated zone be used to represent colloid retardation through unsaturated fractures, or 2) conservatively assume no colloid retardation in fractures. ANL-NBS-HS-000026, Rev 00 44 April 2000 Table 6. Cumulative Probabilities for Colloid Transport between One Matrix Unit and Another Colloid Size (nm) Units 2000 1000 450 200 100 50 6 TMN (TSW4) 1.00 0.92 0.87 0.81 0.71 0.55 0.31 TLL (TSW5) 1.00 0.80 0.79 0.70 0.61 0.51 0.19 TM2 (TSW6) 1.00 0.94 0.90 0.82 0.65 0.51 0.21 TMN1 (TSW7) 1.00 0.99 0.99 0.99 0.93 0.68 0.36 PV3 (TSW8) 1.00 0.98 0.96 0.94 0.90 0.89 0.68 PV2 (TSW9) 1.00 0.72 0.57 0.47 0.39 0.35 0.22 BT1a (CH1) 1.00 0.91 0.89 0.87 0.85 0.83 0.53 CHV 1.00 0.58 0.49 0.43 0.39 0.36 0.07 CHZ 1.00 0.79 0.76 0.73 0.68 0.56 0.30 BT (CH6) 1.00 0.95 0.94 0.92 0.92 0.85 0.40 PP4 1.00 0.99 0.99 0.98 0.98 0.96 0.32 PP3 1.00 0.49 0.34 0.26 0.21 0.16 0.07 PP2 1.00 0.91 0.86 0.81 0.65 0.53 0.22 PP1 1.00 0.79 0.68 0.63 0.57 0.48 0.21 BF3 1.00 0.97 0.94 0.83 0.74 0.66 0.14 BF2 1.00 0.98 0.97 0.96 0.96 0.83 0.25 DTN: LA0003MCG12213.002 Table 7. Cumulative Probability Density Functions for the Retardation of Colloids Through Saturated Fractured Rock R (retardation factor) Probability 1.06 0.0105 1.1 0.039 6 0.08125 100 0.2605 280 0.7605 800 1 DTN: LA0002PR831231.003. ANL-NBS-HS-000026, Rev 00 45 April 2000 6.3 CODE VERIFICATION In this section, a series of comparisons are presented of the particle-tracking model against analytical solutions or other numerical model results. In all cases, the acceptance criteria for these simulations is based on a visual comparison. The comparison is performed to ensure that the particle-tracking model is reproducing the analytical or other solution result to within the resolution of the plot. The development of more quantitative criteria is difficult because comparison of particle-tracking results with other solutions is subject to considerable random error associated with the number of particles injected. Therefore, scientific judgement in the visual comparison of results is a more appropriate means for evaluating the particle-tracking model. 6.3.1 One-Dimensional Transport with Dispersion and Sorption The first set of simulations are one-dimensional flow problems for which the dispersion and sorption can be tested against analytical solutions. The flow system is a simple one-dimensional grid in the direction of flow. A steady-state flow field is established for a ten-node, onedimensional grid, after which the RTTF particle-tracking technique is used to simulate the transport. In the first suite of one-dimensional tests, a short pulse of particles is injected into the system to test the dispersion transfer function. The responses of these pulse injections of solutes are compared to the time-derivative of Equation 7, derived by Nauman and Buffham (1983, p. 107, Equation 3.47): ( ) p. . . . Pe Pe C .. . .. . - - = 4 1 exp 2 1 2 (Eq. 34) The cumulative RTD is also compiled from the particle statistics at the outlet node and compared to Equation 7. Figure 9 shows that the cumulative RTD solution is reproduced almost exactly for 10 > Pe , corresponding to a grid Peclet number of 1, as long as a sufficient number of particles is used in the simulation. For this transport system, the 20 = Pe results become affected by statistical fluctuations when fewer than 10,000 particles are used, whereas for 100,000 particles, the breakthrough curve is represented quite accurately, as shown in Figure 10. This effect is problem dependent; a sufficient number of particles must be determined empirically for each application of the particle-tracking technique by comparing the results with a simulation using more particles. ANL-NBS-HS-000026, Rev 00 46 April 2000 DTN: LA0002BR12213S.002 NOTE: The family of breakthrough curves of mean travel time of 0.00015 days are for Rf = 1. Figure 9. Cumulative Breakthrough Curves Produced by the RTTF Particle-tracking Technique for One- Dimensional Dispersion (Filled Circles) Compared to Analytical Solution (Solid Curves) Figure 9 also shows that for low enough values of Pe , such as the 10 = Pe case in this example, the RTTF particle-tracking simulation begins to deviate noticeably from the analytical solution. This case corresponds to a grid Peclet number of 1, which is a practical cut-off below which one cannot obtain accurate results using this technique. The solution can be made more accurate using a coarser finite element grid, keeping in mind that the restrictions for obtaining an accurate fluid flow solution still exist. The practical implication of this result is that the RTTF particletracking technique is most useful for advection-dominated problems (large values of Pe ). This will generally be the case for large-scale two- and three-dimensional problems, which necessarily result in large spatial discretization due to computational limitations. ANL-NBS-HS-000026, Rev 00 47 April 2000 DTN: LA0002BR12213S.002 NOTE: The RTD is solute concentration normalized so that the area under the curve is unity. All simulations used Pe=20, Rf=1. Figure 10. Breakthrough Curves Calculated Using Different Numbers of Particles Finally, the agreement of the sorption curve in Figure 9 with the analytical solution (Equation 7) shows that linear, equilibrium sorption is also properly accounted for by the RTTF particletracking model. 6.3.2 One-Dimensional Transport with Matrix Diffusion Figure 11 shows a series of simulations of the RTTF particle-tracking algorithm, along with the solution of Tang et al. (1981, Equation 35) for matrix diffusion with and without dispersion in the fractures. The agreement with the analytical solution is very close for both small or moderate amounts of diffusion (curves a and b), for the case of diffusion and sorption on the fracture surface and in the matrix (curve c), and for dispersion, diffusion, and sorption in both media (curve d). This agreement shows that the two-step RTTF approach outlined in the Mathematical Development section provides an excellent approximation of the combined diffusive-dispersive transport system. When the finite-fracture-spacing option in the FEHM particle-tracking model is used, the model (Figure 12) captures the extremes of behavior ranging from no diffusion, to complete diffusion between fractures, to an intermediate case (shown to agree closely with the infinite spacing curve in the figure). ANL-NBS-HS-000026, Rev 00 48 April 2000 DTN: LA0002BR12213S.002 NOTE: Key to breakthrough curves: a) low matrix diffusion; b) moderate matrix diffusion; c) matrix diffusion and sorption on fracture and in matrix; d) dispersion, matrix diffusion, and sorption on the fracture and matrix. Figure 11. Breakthrough Curves Computed Using the RTTF Particle-tracking Technique for the Infinite Spacing Matrix Diffusion Model, Compared to the Analytical Solution of Tang et al. (1981) (Solid Curves) DTN: LA0002BR12213S.002 NOTE: Key to breakthrough curves: Dashed curves - extreme values of no diffusion, and complete diffusion between fractures. Solid curve - infinite fracture spacing assumption. Dotted curve - finite spacing assumption. Figure 12. Breakthrough Curve Computed Using the RTTF Particle-tracking Technique for the Finite Spacing Matrix Diffusion Model In this simulation, a one-dimensional, saturated flow model was used to establish that over a wide range of transport parameters, the code behaves as anticipated. The two dashed curves reproduce the plug-flow breakthrough curves for transport only in the fractures (no diffusion), and complete solute diffusion to the centerline between fractures, resulting in a transport time governed by the matrix porosity. These breakthrough times agree closely with the expected ANL-NBS-HS-000026, Rev 00 49 April 2000 times for fracture-dominated and matrix-dominated transport. At intermediate diffusion coefficients, a comparison of the infinite fracture spacing solution and the finite spacing solution are shown. For the particular parameters shown, the curves agree with each other, indicating that the characteristic diffusion distance is less than fracture spacing for this case. The significant result for this validation exercise is that the two transfer function models agree closely, despite the completely different formulations and governing equations. This is strong evidence that the finite spacing solution is properly implemented in the particle-tracking model. 6.3.3 Transient Wetting of a 1-D Column In this simulation, a transient, one-dimensional, unsaturated flow system is simulated. Using a relative permeability model proposed by van Genuchten (1980), a model is constructed of the movement of a water front under the influence of gravity and associated tracer particles. The entire flow path is initially set at a saturation below the residual liquid saturation of the medium, and a constant infiltration rate is applied at the inlet of the column. Two particle-tracking simulations are performed: one in which the particles are injected with the infiltrating fluid, and one in which the particles start at a location halfway down the column (cell 54 out of a total of 107). Although many particles are injected, all transport with identical travel times, since this is an advection-only transport case. This validation problem tests the ability of the code to handle the situation of a transient flow field. The exercise is thus important to justify the use of the particle-tracking method for examining radionuclide migration in transient flow fields. When the particles are injected with the infiltrating fluid, the particles should move with the saturation front. Figure 13 shows that the expected behavior is approximated (the filled circles in the plot), with the particles lagging the front edge of the wetting front by two cells. This slight lag of the particle movement is caused by a slightly dispersed wetting front that does not maintain a sharp interface. For the situation in which the particles are introduced halfway down the length of the column (the X’s), the particles remain stationary until the saturation front reaches the particles, after which time the particles move with the front as before. The figure demonstrates that for this transient flow case, the RTTF particle-tracking approach provides an adequate solute transport simulation. ANL-NBS-HS-000026, Rev 00 50 April 2000 DTN: LA0002BR12213S.002 NOTE: Key to plot. Curves - saturation profiles at regular time intervals. Filled circles - locations of particles started at the entrance. X’s - particles started part way down the column. Figure 13. Movement of Particles with the Saturation Front for the Transient Wetting Case 6.3.4 Decay Calculations Using the Integral Approach A C++ software routine TEST was written to test the integration approach against the discrete formula in Equation 31. The software routine was checked following AP-SI.1Q procedures for correctness. The results showed that routine TEST performed as designed. The source code and the checking process are documented in Attachment I. Theoretically, the accuracy of the integration approach depends on two factors, the number of particles injected into the system, and the particle injection distribution. In general, large number of particles and uniform injection of particles in time domain increases the accuracy of the integration approach. In a typical performance assessment simulation, on the order of 1,000,000 particles might be released for each species over the repository area (represented by about 500 nodes) during a time period of 10,000 to 1,000,000 years, and the time step used in FEHM simulations is of order 1,000 years. To test the accuracy of the integration approach, the software routine TEST was run by injecting 1,000 particles into the repository during a period of 1,000 years. The results were used to compare the integration approach against the discrete approach. Radionuclides 239Pu (half life of 2.406E4 years) and 237Np (half life of 2.14E6 years) (CRC 1991) were chosen for these tests. Figures 14 and 16 compare the numbers of particles decayed using the two approaches for 237Np and 239Pu, respectively. The results show good agreement between the integration approach and the direct discrete calculation. Figures 15 and 17 are the percentage errors of integration approach corresponding to discrete calculation results at different times. The large relative errors at early times are due to the fact that few particles have decayed. At greater times, the relative percentage error decreases. The relative error can also be interpreted as a time delay of the integration approach corresponding to the discrete calculations in the time domain. The figures also show that increasing the number of particles released during the same time period increases the accuracy of the integration approach. ANL-NBS-HS-000026, Rev 00 51 April 2000 Therefore, in general, the results from the integral approach are very close to those based on the discrete calculations in Equation 31. The relative error for the majority of the runs was less than 1%. However, the computational burden was greatly reduced from millions of operations to just 10 operations within each time step, and thus is suitable for use in the particle-tracking simulations. Precautions should be taken to ensure that a sufficient number of particles are injected to avoid error. This is best done by increasing the number of particles until the results no longer show significant change from a simulation with fewer particles. DTN: LA0002BR12213S.002 Figure 14. Number of 237Np Particles Decayed Based on the Integration Approach and the Discrete Calculations Time (year) 103 104 105 106 Np237, number of particles decayed 100 101 102 103 104 discrete_1,000 particles integration_1,000 particles discrete_10,000 particles integration_10,000 particles ANL-NBS-HS-000026, Rev 00 52 April 2000 Time (year) 103 104 105 106 Np237, % error corresponding to discrete calculations -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 integration_1,000 particles integration_10,000 partcles DTN: LA0002BR12213S.002 Figure 15. 237Np Integration Approach Percentage Errors Corresponding to Discrete Calculations ANL-NBS-HS-000026, Rev 00 53 April 2000 DTN: LA0002BR12213S.002 Figure 16. Number of 239Pu Particles Decayed Based on the Integration Approach and the Discrete Calculations Time (year) 103 104 105 106 Pu239, number of particles decayed 101 102 103 104 discrete_1,000 particles integration_1,000 particles discrete_10,000 particles integration_10,000 particles ANL-NBS-HS-000026, Rev 00 54 April 2000 Time (year) 103 104 105 106 Pu239, % error corresponding to discrete calculations -40 -30 -20 -10 0 10 integration_1,000 particles integration_10,000 particles DTN: LA0002BR12213S.002 Figure 17. 239Pu Integration Approach Percentage Errors Corresponding to Discrete Calculations 6.3.5 Decay/Ingrowth Validation Simulations In this test suite, simulation results from FEHM particle-tracking model with decay and ingrowth verified against the CHAIN semi-analytical solutions for a 4 species chain decay-ingrowth model (van Genuchten 1985). The use of this routine, the qualification procedures, source code, and test results are listed in Attachment III. The method of comparison for these runs a visual comparison of the plotted results. These runs are designed to test the accuracy of the particle-tracking model, and should be interpreted as representative results of the Yucca Mountain system. The Software routine BKPM was used to post process FEHM simulation results to extract breakthrough information for comparison with CHAIN analytical results. Routine BKPM was qualified according to AP-SI.1Q. The use of the BKPM routine, the qualification procedures, source code, and test results are listed in Attachment II. For all comparisons of CHAIN and FEHM for decay-ingrowth simulations, a flow system was developed with the following attributes: • Saturated steady state flow in a 1-D system • Porosity of 0.3 • Average pore-water velocity was 1.05192x10-4 m/year ANL-NBS-HS-000026, Rev 00 55 April 2000 • Solute with constant injection concentration of 1.0 mol/l released from 0 to 5,000 years at origin x=0. Breakthrough data were collected at x=1 m down stream. • FEHM grid resolution: 0.005 m • Longitudinal dispersivity of 0.005 m. The first comparison is for a conservative tracer with no decay and no sorption. The purpose of this initial run is to lay a foundation for FEHM decay-ingrowth model verifications by ensuring that the two codes yield similar transport behavior in the absence of decay and ingrowth. The extracted FEHM particle breakthrough curve is plotted together with the breakthrough curve calculated by CHAIN in Figure 18. The close agreement of the FEHM and CHAIN breakthrough curves proves that the two codes compute the transport of a conservative solute correctly. Time (year) 0 5000 10000 15000 20000 Relative concentration 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 CHAIN FEHM DTN: LA0002BR12213S.002 Figure 18. Comparison of Base Case Transport Results: CHAIN and FEHM The second test uses the same flow system, but employs a 4 member decay-ingrowth chain. The sequential decay-ingrowth chain was formed by Ra Th U Pu 226 230 234 238 . . . . The decay rates for the 4 species are 0.0079, 0.0000028, 0.0000087, and 0.00043 yr-1 for 238Pu, 234U, 230Th, and 226Ra, respectively (van Genuchten 1985, p. 144, appendix 5). The grid Peclet number of the system was 1.0. Retardation factors were set to 1.0 for 238Pu and 234U, but to avoid some ANL-NBS-HS-000026, Rev 00 56 April 2000 parameters being divided by 0 in software routine CHAIN, the retardation factors for 230Th and 226Ra were set to 1.001 in the CHAIN input file. The initial concentrations for all species were 0. From time 0 to 5,000 years, only 238Pu was released into the system at the origin, with a constant concentration of 1.0. In FEHM simulations, the 5,000 year release period was divided into 50 segments so that each segment corresponding to 100 years. Within each segment, 10,000 particles were injected into the system uniformly over the 100 year period. There were no source terms for 234U, 230Th, and 226Ra. The CHAIN breakthrough data and the post processed FEHM breakthrough data were plotted together in Figure 19. Due to the short half life of 238Pu and the relatively long half lives of 234U and 230Th, the concentrations of 238Pu and 226Ra were extremely low and were omitted in this figure. The figure exhibits good agreement between FEHM and CHAIN breakthrough curves. DTN: LA0002BR12213S.002 Figure 19. Comparison of CHAIN and FEHM Transport Results for a 4-Member Decay-Ingrowth Chain, Ra Th U Pu 226 230 234 238 . . . . Peclet Number = 1.0 In the third test case, a pseudo sequential decay chain with species 4 3 2 1 . . . , with half lives for species 1 through species 4 of 10,000, 3,000, 10,000, and 4,000 years, respectively. The transport process was dominated by advection and dispersion only with a grid Peclet number of 1.0. The particle injection pattern used in this simulation was the same as for case 1 and case 2. The purpose of this problem is to select half-lives that allow all species to exhibit a visible breakthrough at the exit, thus allowing a more complete comparison to be made. The CHAIN Time (year) 0 5000 10000 15000 20000 Relative concentration 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 CHAIN species 2 CHAIN species 3 FEHM species 2 FEHM species 3 ANL-NBS-HS-000026, Rev 00 57 April 2000 and FEHM breakthrough curves are compared in Figure 20. All species agree fairly closely with one another, with only small differences for species 3 and 4, especially near the peak. The differences are probably due to either numerical truncation errors in CHAIN for approximation errors in FEHM. Nevertheless, the visual comparison yields no systematic differences, so the agreement is satisfactory. Time (year) 0 5000 10000 15000 20000 25000 Relative concentration 0.0 0.1 0.2 0.3 0.4 0.5 CHAIN species 1 CHAIN species 2 CHAIN species 3 CHAIN species 4 FEHM species 1 FEHM species 2 FEHM species 3 FEHM species 4 DTN: LA0002BR12213S.002 Figure 20. Comparison of CHAIN and FEHM Transport Results for a Case of 4-Member Sequential Decay Chain without Retardation. Peclet Number = 1.0 The fourth test case compares FEHM and CHAIN for a single solute with a retardation factor of 1.9. All other parameters are similar to the previous test cases. The pupose of this run is to ensure that the codes are handling sorbing solutes in the same manner as before, including decay and sorption. The comparison shown in Figure 21 shows that this is the case. ANL-NBS-HS-000026, Rev 00 58 April 2000 Time (year) 0 5000 10000 15000 20000 25000 30000 35000 40000 45000 Relative concentration 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 CHAIN FEHM DTN: LA0002BR12213S.002 Figure 21. Comparison of CHAIN and FEHM Transport Results for the Base Case with a Retardation Factor of 1.9 and a Peclet Number of 1.0 The fifth and final test case of this suite includes sorption and decay. The retardation factors for species 1 through species 4 were 1.0, 1.0, 1.9, and 1.001, respectively. The boundary and initial conditions and half-lives were the same as those of case 3. The FEHM particle injection pattern was the same as for case 3. The FEHM and CHAIN breakthrough curves were plotted in Figure 22. In general, good agreements were observed between FEHM and CHAIN curves. Taken together, this suite of verification simulations demonstrates that the particle-tracking model accurately handles decay chains and, thus, is suitable for use in simulating multiple radionuclides in the UZ transport model. ANL-NBS-HS-000026, Rev 00 59 April 2000 Time (year) 0 5000 10000 15000 20000 25000 Relative concentration 0.0 0.1 0.2 0.3 0.4 0.5 CHAIN species 1 CHAIN species 2 CHAIN species 3 CHAIN species 4 FEHM species 1 FEHM species 2 FEHM species 3 FEHM species 4 DTN: LA0002BR12213S.002 Figure 22. Comparison of CHAIN and FEHM Transport Results for a Case with a 4-Member Decay- Ingrowth Chain and a Retardation Factor of 1.9 for Species 3. Peclet Number = 1.0 6.3.6 EBS Random Release Model Verification The FEHM EBS random release model was tested in standalone mode to verify that the model worked as designed. The model was tested with a “pseudo problem” to make sure the model selects the right nodes. The test results from the test problem were checked manually for correctness. The EBS random release model was tested for selecting early failed package nodes based on given coordinates, and for randomly selecting failed package nodes based on given number of failed packages within a sub-region. In the discussion below, the term “M_fine” is used to describe the category of packages that fail initially, and “N_large” is used to describe the regions in which the remaining package failures occur. They are used as descriptive terms, and do not refer to the number of packages of either type that fail. To perform the test, a routine called test_ebs_random.f Version 1.0, was written in Fortran 90. The core of the test routine is the EBS random release model, with code added to generate random node locations within a given field, loading nodes into an array containing node coordinates, and randomly assigning each node an index. When started, the user inputs the total number of nodes, number of sub-regions, random seed, number of M_fine individually failed packages, and number of species, and the on/off switch for releasing particles at the current time ANL-NBS-HS-000026, Rev 00 60 April 2000 step. Then, the test routine generates randomly located nodes in each sub-region, loading those nodes into an array, and randomly assigns locations for the M_fine early failed packages. The test routine next reads in the controlling parameter for number of simulation time steps and the number of failed N_large packages in each sub-region at each time step. Once the nodes are generated and the parameters read in, test_ebs_random.f starts the core part of the calculation, the EBS random release model, to select nodes for the M_fine early failed packages and randomly select release nodes for the N_large failed packages in each sub-region. To test the release model, a system was constructed consisting of two rectangular sub-regions defined by range of the (x,y) coordinates. Each sub-region has 20 nodes. Within each sub-region there are three early failed packages defined by location coordinates. A total of four time steps are simulated. At each time step, the number of N_large failed packages increase with time. The EBS random release model was used to select release nodes at each time step. The input parameters used for this test problem are listed in Table 8. Table 8. List of Test Problem Input Parameters: EBS Random Release Model Verification Parameters Sub-region 1 Sub-region 2 Rectangular Sub-region Range (x_min,x_max)=(0,100) (y_min,y_max)=(0,100) (x_min,x_max)=(120,250) (y_min,y_max)=(120,250) Number of Nodes 20 20 M_fine Early Failed Packages 3 3 Time Step 1: N_large Failed Packages 10 10 Time Step 2: N_large Failed Packages 15 15 Time Step 3: N_large Failed Packages 17 17 Time Step 4: N_large Failed Packages 20 20 The M_fine early failed package nodes and the randomly located N_large failed package nodes at different time steps are plotted in Figures 23 to 26. At time step 1, there were 10 failed N_large packages in each sub-region containing 20 nodes. The EBS random release model should find nodes corresponding to the M_fine early failed packages and randomly select N_large number. The results from this test indicate that showed that the algorithm successfully generates node locations and loads the generated nodes into the node coordinate array. The node index of the generated nodes are not in sequential order but rather randomly assigned to mimic the fact that node index for repository nodes may not be in sequential order. The test code also successfully generates 3 early failed package locations in each sub-region, as expected. The selected M_fine early failed package nodes and the randomly located N_large failed package nodes at different time steps are plotted in Figures 23 to 26. Figure 23 shows the locations of nodes and the M_fine early failed packages. ANL-NBS-HS-000026, Rev 00 61 April 2000 x (m) -50 0 50 100 150 200 250 300 y (m) -50 0 50 100 150 200 250 300 sub-region boundary sub-region nodes M_fine failed packages sub-region 1 sub-region 2 DTN: LA0002BR12213S.002 Figure 23. Distribution of Nodes and M_fine Failed Packages in the Test Field The EBS random release model was run to select nodes that are closest to the given locations of the early failed packages. At time step 1, there were 10 failed N_large packages in each subregion containing 20 nodes. The EBS random release model should find nodes corresponding to the M_fine early failed packages and randomly select N_large number of nodes for the N_large failed packages in each sub-region. Figure 24 shows that the EBS random release model successfully selects nodes for the three early failed M_fine packages and randomly selects 10 nodes from the remaining 17 nodes for the N_large failed packages in each sub-region. ANL-NBS-HS-000026, Rev 00 62 April 2000 x (m) -50 0 50 100 150 200 250 300 y (m) -50 0 50 100 150 200 250 300 sub-region boundary sub-region nodes M_fine failed packages N_large failed packages sub-region 1 sub-region 2 DTN: LA0002BR12213S.002 Figure 24. Distribution of Selected M_fine and N_large Failed Packages at Time Step 1 At time step 2, there are a total of 15 failed N_large packages in each sub-region. Thus, the number of newly failed N_large packages in each sub-region during this time step is five. The three early failed package nodes and the 10 previously selected nodes were not altered once they were selected. The EBS random release model then randomly selected five new nodes from the seven nodes remaining. Figure 25 shows the selected nodes failed during this time step. ANL-NBS-HS-000026, Rev 00 63 April 2000 x (m) -50 0 50 100 150 200 250 300 y (m) -50 0 50 100 150 200 250 300 sub-region boundary sub-region nodes M_fine failed packages N_large failed packages sub-region 1 sub-region 2 DTN: LA0002BR12213S.002 Figure 25. Distribution of Selected M_fine and N_large Failed Packages at Time Step 2 At time step 3, the total number of failed N_large packages increased from 15 for the previous time step to 17 during the current time step. The number of newly failed N_large packages during this time step is two. The EBS random release model was used to select two nodes in each sub-region. Figure 26 shows the selected nodes during this time step. ANL-NBS-HS-000026, Rev 00 64 April 2000 x (m) -50 0 50 100 150 200 250 300 y (m) -50 0 50 100 150 200 250 300 sub-region boundary sub-region nodes M_fine failed packages N_large failed packages sub-region 1 sub-region 2 DTN: LA0002BR12213S.002 Figure 26. Distribution of Selected M_fine and N_large Failed Packages at Time Step 3 At time step 4, although the total number of failed N_large failed packages was 20, there are no available nodes for selection at this time step. All nodes in each sub-region failed. The distribution of failed nodes at this time step is therefore same as that at time step 3. In summary, the results of this manual test indicates that the code is performing as expected and designed. 6.4 DUAL-PERMEABILITY BENCHMARKING SIMULATIONS In Section 6.3 test cases were presented comparing the particle-tracking transport model against analytical solutions in a saturated, single-continuum system. Here a series of one-dimensional, dual-permeability calculations are carried out under unsaturated conditions. In all cases, the breakthrough curve at the water table is used to judge the adequacy of the model. Since the intended application for the UZ transport model is to simulate breakthrough curves at the water table, these tests exercise the code options in a manner equivalent to the way the code will be used in the TSPA calculations. Analytical solutions for dual-permeability systems are difficult to obtain, even in one dimension, for any but the simplest cases. Therefore, the particle-tracking model for dual permeability, unsaturated flow and transport is compared to a discrete fracture model that uses FEHM’s finite element transport module. In addition, simulation results are presented for the one-dimensional, ANL-NBS-HS-000026, Rev 00 65 April 2000 dual-permeability model, also an option in FEHM, with transport solved using a finite element solution to the transport equation. These comparisons serve to illustrate the accuracy and, in some cases, assess the limitations of the particle-tracking model compared to other transport techniques. In this suite of benchmarking runs, a one-dimensional system is selected because in the discrete fracture model, a one-dimensional, vertical fracture is easily discretized to capture the gradients in fluid pressure and solute concentration. The schematic in Figure 3 represents the model system geometry being used in these runs. The level of discretization used in this discrete fracture model is impossible to achieve in a site-scale model, but serves as an idealized benchmark for judging the suitability of a simplified particle-tracking model or other transport model. If the benchmarking is successful, the complexities in the three-dimensional flow system can be captured with a dual-permeability flow field, while the particle-tracking model can be used to capture the detailed physics of transport with matrix diffusion and sorption at the smaller scale. Our primary criterion for judging the adequacy of the particle-tracking model is that it adequately capture the behavior simulated in the discrete fracture model over a wide range of transport model parameters, for a wide range of flow system behavior. If the particle-tracking model can pass this test, its advantage of computational efficiency can be exploited in site-scale transport simulations. The purpose of also presenting one-dimensional, dual-permeability transport model results using a finite element solution is to demonstrate the relative ability of the more traditional model formulation to mimic the behavior of the discrete fracture model. Practically speaking, the dual-permeability transport simulator is one of the few alternative formulations available if a three-dimensional flow model is used directly for transport simulations. The benchmarking system consists of a two-dimensional discrete fracture model grid that extends from an elevation of 630 m to 1150 m (meant to represent roughly the top of the TSw). The grid has a vertical spacing of 2.5 m, and a variable horizontal spacing that goes from less than 1 mm to about 0.9 m at the opposite side of the grid from the fracture. The model domain to the center line between fractures (B in Figure 3) is 5 m. This system is designed to simulate half of an open fracture and the adjacent matrix block out to the middle position between fractures (i.e. fracture spacing is 10 m). Four flow situations are considered, each modeled by both the discrete-fracture and the dualpermeability models: Case 1: Permeable fracture, impermeable matrix. Case 2: Permeable fracture and matrix, with parallel flow in the fractures and matrix, but little or no fluid flow between the fracture and matrix from the radionuclide release point to the water table. Case 3: Permeable fracture, varying fluid fluxes in the fractures and matrix in three hydrologic units. Transitions from fracture flow to matrix flow (and back) occur abruptly very close to the unit interfaces for the parameter values chosen. ANL-NBS-HS-000026, Rev 00 66 April 2000 Case 4: Permeable fracture and matrix, with a flow restriction at the fracture/matrix interface such that fluid fluxes in the fractures gradually decrease with depth over the entire length of the model domain. All cases use 5 mm/y as an average flux along the 5 m top surface. In the discrete fracture model, fluid is injected only in the fracture and allowed to imbibe into the matrix. At the radionuclide release level, water travels in parallel in the fracture and matrix (except for Case 4). The discrete fracture model domain was extended about 100 m below the assumed water table elevation of 730 m so as to not introduce any boundary-condition related artifacts at the 730 elevation. In the transport simulations, a row of nodes was selected 730 m, and a solute sink was placed there. Solute mass was introduced only at the fracture node at 1087.5 m, by fixing the concentration at a value of unity. The solute mass arrival rate at the water table plane is recorded in the solute mass balance output in the *.out file of FEHM. This time-varying curve was normalized by the injection solute mass rate. At long times, the input and output mass rates equal each other, and the breakthrough curve approaches unity. For the one-dimensional, dual-permeability simulations, the same 5 mm/y rate was used, but in Cases 2 and 3 the flow distribution (fracture and matrix) was adjusted at the top of the model to ensure that the flow fraction at the release level mimicked the discrete fracture model. Particletracking model results in the form of arrival times for each particle at the water table are recorded in the *.fin file for post-processing to obtain the cumulative breakthrough curve. The one-dimensional, dual-permeability model simulations performed with the finite element transport simulation module are performed using the same flow fields as the particle-tracking models, and the post-processing is executed similarly to the discrete fracture model runs. The nodiffusion cases used 1.E-30 for the diffusion coefficient, while the diffusion cases assume 3.2E- 11 m2/s, and in some cases 1.0E-12 m2/s. A longitudinal dispersivity of 1 m was assumed throughout the model domain in the vertical direction. 6.4.1 Case 1 Results The first case isolates flow to only the fracture by applying an extremely low matrix permeability. In this way, the transport solution can be checked against the analytical solution for matrix diffusion. The analytical solution of Tang et al. (1981) is used for comparison of both the particle-tracking and discrete-fracture models. Figure 27 shows the breakthrough curves for the particle-tracking and discrete-fracture models compared to the analytical solution. With no diffusion into the matrix, breakthrough times are on the order of one year. The particle-tracking and analytical solutions agree very closely, showing that the particle-tracking model with dispersion is implemented properly. The deviation of the discrete fracture model from the analytical result is caused by numerical dispersion, but the breakthrough time agrees with the analytical solution. ANL-NBS-HS-000026, Rev 00 67 April 2000 DTN: LA0002BR12213S.002 Figure 27. Breakthrough Curves for Case 1 Comparing the Discrete Fracture Model with the Particle-tracking Method and an Analytical Solution When diffusion into the stagnant fluid in the matrix is included, breakthrough times ranging from about 100 to 10,000 years are obtained, and again the particle-tracking and analytical solutions agree closely. The particle-tracking model follows the discrete fracture model closely for most of the breakthrough curve, deviating at long times. Since the infinite fracture spacing solution was used in this comparison, this departure at long times is expected because at late times, diffusion to the centerline between fractures will result in somewhat shorter travel times than if mass can diffuse greater distances away from the fracture than the half-spacing B. One factor to note about the particle-tracking solution is that the aperture required for the model to agree with the analytical solution under unsaturated conditions is the total aperture times the fracture saturation. This factor must be considered when using the particle-tracking model to perform unsaturated transport simulations. This test case, simple enough to invoke an analytical solution, ensures that the numerical parameters associated with the discrete fracture model are suitable for performing comparisons to the particle-tracking model. 6.4.2 Case 2 Results The second case considers parallel flow in the matrix and fracture. Water is injected at the top of the model only in the fracture node. It redistributes at steady state such that there is parallel flow in the fracture and matrix. This redistribution at the top occurs at a higher elevation than the repository horizon of 1087 m, so that the flow field relevant to transport is parallel flow in the fracture and matrix, with virtually no water flux between fracture and matrix along the transport ANL-NBS-HS-000026, Rev 00 68 April 2000 pathway. Thus, the difference between Case 1 and Case 2 is that for Case 2, there is somewhat less flux in the fracture and some flux in the matrix, rather than stagnant fluid. The discrete fracture and particle-tracking results without diffusion agree closely (Figure 28), as expected since this simulation isolates transport to the fractures, making it very similar to the nodiffusion results of Case 1. For a diffusion coefficient of 3.2E-11 m2/s, three results are presented in Figure 28: the discrete fracture model, particle tracking with a infinite fracture spacing assumption, and particle tracking with the finite spacing assumption. The improved agreement of the finite spacing result (compared to the infinite spacing case) indicates that the matrix diffusion model with finite spacing captures the transport behavior more adequately than the infinite spacing case. The particle-tracking model agrees closely with the discrete fracture benchmark at early times, and departs from it somewhat at later times. This discrepancy is caused by the inherent assumption in the transfer function for matrix diffusion that the matrix fluid is stagnant. When enough solute diffuses into the matrix to result in travel times that approach the advective travel time in the matrix, the assumptions in the particle-tracking model are no longer strictly valid. Under those conditions, in the discrete fracture flow system, advection in the matrix carries solute to the outlet somewhat more rapidly than the particle-tracking model. DTN: LA0002BR12213S.002 Figure 28. Breakthrough Curves for Case 2 Comparing the Discrete Fracture Model with the Particle-tracking Method For comparison, the same simulation results are presented in Figure 29, but transport results using the finite element, dual-permeability transport module are also included. All of the models agree adequately without diffusion, but the one-dimensional, dual-permeability solution predicts much earlier breakthrough for this case. At early times, the dual-permeability solution fails to match the discrete fracture model because the fracture/matrix diffusion term is underestimated, and mass remains in the fracture to an unrealistic degree. Conversely, this model does a better ANL-NBS-HS-000026, Rev 00 69 April 2000 job than the particle-tracking model at late times, fracture and matrix concentrations have equilibrated for this particular transport system. Nevertheless, the particle-tracking model adequately reproduces the discrete fracture model results over the entire range of travel times. When the diffusion coefficient is even smaller (Figure 30, diffusion coefficient of 1.0E-12 m2/s), the finite element, dual-permeability transport model completely fails to capture the behavior of the discrete fracture model, whereas the particle-tracking model continues to provide a good representation of the behavior of the system. DTN: LA0002BR12213S.002 Figure 29. Breakthrough Curves for Case 2 Comparing the Discrete Fracture Model with the Particle-tracking Method and the Finite Element, Dual-Permeability Solution To summarize this test case, systems with small to moderate amounts of diffusion are well captured (using the comparison to the discrete fracture model as a benchmark) by the particletracking algorithm, providing a superior transport model to a dual-permeability transport solution. The method avoids the artificial early-time breakthrough of mass that the dualpermeability transport solution yields due to the inability of the latter to capture sharp gradients in concentration away from the fractures at early times. At late times, the particle-tracking method deviates somewhat from the discrete fracture model. In estimates of unsaturated zone system performance, this deviation is probably more tolerable than the severe overestimate of early breakthrough of mass in the dual-permeability transport solution, but of course this factor must be judged in light of the goals of the model analysis being performed. ANL-NBS-HS-000026, Rev 00 70 April 2000 DTN: LA0002BR12213S.002 Figure 30. Breakthrough Curves for Case 2 (Low Diffusion) Comparing the Discrete Fracture Model with the Particle-tracking Method and the Finite Element, Dual-Permeability Solution 6.4.3 Case 3 Results This test case simulates the transition from predominantly fracture flow in the upper unit (given properties similar to the TSw) to predominantly matrix flow from an elevation of 898 m to 833.5 m (CHnv-like properties) to significantly more fracture flow from 833.5 m to the water table (CHnz-like properties). This case is realistic for testing the dual-permeability particle-tracking algorithm for a system containing several hydrogeologic units of contrasting properties, such as the Yucca Mountain unsaturated zone. To facilitate the model comparison, the discrete fracture model assumes identical properties within the entire length of the fracture, and does not attempt to account for contrasting fracture volume fraction and fracture spacing from unit to unit. Likewise, the dual-permeability model uses constant fracture properties throughout the model. Figure 31 compares the simulation results for the three model options assuming no matrix diffusion. The particle-tracking and the dual-permeability finite element solution agree with one another closely, and compare fairly well with the discrete fracture model solution. The discrepancy with the discrete fracture simulation can be attributed to the inability to exactly match the advective velocities and detailed flow field of the discrete fracture model near the unit interfaces where transition from fracture to matrix (or matrix to fracture) flow occurs. This is a natural limitation of this benchmarking exercise for such a flow field, and cannot be easily overcome. Nevertheless, despite this disagreement, the characteristic behavior of the system is reproduced with both the dual-permeability transport model and the particle-tracking model. A minimum travel time through the system of the order of 1000 years is obtained due to matrix flow through the middle layer. ANL-NBS-HS-000026, Rev 00 71 April 2000 DTN: LA0002BR12213S.002 Figure 31. Breakthrough Curves for Case 3 (No Diffusion) Comparing the Discrete Fracture Model with the Particle-tracking Method and the Finite Element, Dual-Permeability Solution Figure 32 is the same comparison of the three models for a matrix diffusion coefficient of 3.2E- 11 m2/s. The particle-tracking model exhibits a somewhat later arrival time distribution than the discrete fracture model, whereas the finite element dual-permeability transport solution exhibits an earlier arrival. The dual-permeability transport solution probably better represents the overall breakthrough curve than does the particle-tracking result. However, the problem of rapid transport of some mass through the fractured units has not been eliminated (refer to Figure 29) but merely masked by the long travel time through the matrix-flow dominated unit. On the other hand, the particle-tracking solution tends to predict somewhat longer residence times than the discrete fracture model when diffusion is introduced. As in Case 2, this is the result of the assumption of diffusion into stagnant liquid as formulated using the transfer function approach. When diffusion is extensive enough to result in travel time delays that approach the travel time of fluid through the matrix block, the particle-tracking method as currently formulated will overestimate travel times. This conclusion is in contrast to the excellent agreement of the particle-tracking model to the discrete-fracture model for small to moderate amounts of diffusion. ANL-NBS-HS-000026, Rev 00 72 April 2000 DTN: LA0002BR12213S.002 Figure 32. Breakthrough Curves for Case 3 (High Diffusion) Comparing the Discrete Fracture Model with the Particle-tracking Method and the Finite Element, Dual-Permeability Solution 6.4.4 Case 4 Results This test case is designed to examine a situation in which fluid is imbibed from the fracture to the matrix over the entire length of the fracture, including the portion of the model domain where the transport is occurring. The first step in developing this test case is to obtain discrete fracture and dual-permeability flow fields that match one another. Because our purpose is to compare transport models for comparable flow fields, the exact approach to obtain flow fields that agree is not important, as long as the same fluid flux values are achieved in the fracture and matrix at all locations in the model. The basic model consists of a single matrix unit with uniform properties, as in Case 2. To achieve the flow field for the discrete fracture model, a flow resistance was applied between the fracture nodes and the adjacent matrix, so that instead of redistributing within the top 10 to 20 meters of the model, flow imbibes from the fracture into the matrix along the entire length of the model. The dual-permeability flow field used parameters similar to those in Case 2 except for a larger matrix permeability and the injection of all inlet fluid into the fracture node instead of distributing it between the fracture and matrix. By simultaneously adjusting the matrix permeability in the dual-permeability model and the resistance factor in the discrete fracture model, comparable flow fields were obtained. Figure 33 shows the relative flow rate in the fracture portion of each model versus depth, illustrating that comparable flow fields were obtained for performing the transport runs. ANL-NBS-HS-000026, Rev 00 73 April 2000 DTN: LA0002BR12213S.002 Figure 33. Fracture Flux Fraction Versus Elevation for Case 4 Flow Fields: Discrete Fracture Model and One-Dimensional, Dual-Permeability Model Figure 34 shows the results of the three transport models for the case of no diffusion. In all cases, early arrival times of a portion of the injected solute are observed due to the persistent fracture flow along the entire length of the system, and the absence of diffusion to bring mass into the matrix. The height of the plateau in the breakthrough curve represents the ratio of the flow rate of fluid reaching the water table in the fracture to the fracture flow rate at the radionuclide release location. The slight mismatch of the particle-tracking and dual-permeability transport models with the discrete fracture benchmark result is due to the difficulty of achieving an exact representation of the discrete fracture flow field with the one-dimensional, dual-permeability model. Nevertheless, the models are quite consistent with one another, illustrating that the particle-tracking model reproduces the idealized behavior of the system for this flow system for advective transport without diffusion. ANL-NBS-HS-000026, Rev 00 74 April 2000 DTN: LA0002BR12213S.002 Figure 34. Breakthrough Curves for Case 4 (No Diffusion) Comparing the Discrete Fracture Model with the Particle-tracking Method and the Finite Element, Dual-Permeability Solution Figure 35 shows the results of the model for the low diffusion case (diffusion coefficient of 1.E- 12 m2/s), and Figure 36 compares the models for the higher diffusion coefficient of 3.2E-11 m2/s. These figures suggest that the particle-tracking model continues to provide an adequate representation of the transport system for lower diffusion, but departs from the discrete-fracture model when higher diffusion is used. In contrast, the dual-permeability transport solution fails to capture the early-time behavior of the breakthrough curve, yielding an artificial early arrival of some solute mass. The extent of this problem is lessened for the higher diffusion coefficient case. Finally, Figure 37 shows the same model parameters as Figure 35, but includes sorption on the matrix rock (Kd = 5 mL/g). The early arriving mass continues to be present in the dualpermeability transport solution. Note the scale change in Figure 37 when evaluating the model behavior compared to that presented in Figure 35. The particle-tracking model also predicts an earlier first arrival of solute than the discrete fracture model result, but to a much lesser degree than the dual-permeability transport solution. The particle-tracking model could be improved in the future by adding a term to account for diffusive transport that carries mass into the matrix domain. For now, the transport result with diffusion and matrix sorption can be considered conservative for portions of a flow domain in which significant fracture-to-matrix flow occurs along the entire flow path. In dual-permeability models in which the transition occurs abruptly at the interface, the error will be isolated to nodes over which the flow transition occurs. ANL-NBS-HS-000026, Rev 00 75 April 2000 DTN: LA0002BR12213S.002 Figure 35. Breakthrough Curves for Case 4 (Low Diffusion) Comparing the Discrete Fracture Model with the Particle-tracking Method and the Finite Element, Dual-Permeability Solution DTN: LA0002BR12213S.002 Figure 36. Breakthrough Curves for Case 4 (High Diffusion) Comparing the Discrete Fracture Model with the Particle-tracking Method and the Finite Element, Dual-Permeability Solution ANL-NBS-HS-000026, Rev 00 76 April 2000 Note: 1 cc = 1mL DTN: LA0002BR12213S.002 Figure 37. Breakthrough Curves for Case 4 (Low Diffusion, Matrix Sorption) Comparing the Discrete Fracture Model with the Particle-tracking Method and the Finite Element, Dual-Permeability Solution 7. CONCLUSIONS The particle-tracking algorithm developed for this AMR incorporates the transport processes determined to be relevant in the site characterization program, including advection, dispersion, sorption, and matrix diffusion. In addition, new model development was required to allow for finite spacing between fractures in the matrix-diffusion model, multiple-species transport with decay/ingrowth, and the integration with the TOUGH2 and GoldSim applications. These capabilities were incorporated into the current version of FEHM, and validation testing demonstrated that the various processes are captured in the code. Therefore, this version of the code can be used to perform the UZ transport calculations for TSPA-SR. However, several caveats discussed below must be considered when executing the model runs to predict UZ radionuclide transport. The particle-tracking method developed in this AMR, called the Residence Time Transfer Function (RTTF) particle-tracking technique, employs a new, efficient particle-tracking algorithm that is suitable for performing large-scale transport simulations. Like most particletracking methods, this algorithm provides advection-dominated solutions that are free of numerical dispersion. Thus, the technique can track a sharp front in solute concentration without the usual numerical dispersion or concentration-profile oscillations encountered in conventional finite-difference or finite-element solutions of the advection-dispersion equation. The development that distinguishes the RTTF particle-tracking method from existing groundwater transport particle-tracking techniques is the conceptualization of transport as the ANL-NBS-HS-000026, Rev 00 77 April 2000 movement of particles from cell to cell. This approach eliminates the need to resolve the velocity vectors by interpolation at all particle positions. Instead, the mean residence time and the probabilities of travel to adjacent cells are computed once for each cell (until the flow solution changes), regardless of how many particles travel through the cell. Another advantage of the cellto- cell approach is that the information needed to compute particle residence times and pathways is readily available in finite-difference or finite-element solutions of the flow problem. The fluid mass and intercell mass flows rates are a direct result of the fluid flow solution, and thus the method can be implemented without regard to the nature of the numerical grid (structured versus unstructured grids, element shapes, etc.). The particle-tracking technique is also simple to implement for dual-permeability flow models, since the fracture-matrix interaction term is treated as another flow rate to or from a given node. As particles are tracked through the model domain, no particle transport time stepping is necessary; all particles are tracked from the starting time to the ending time of a segment of the simulation. The only reasons to terminate and restart a simulation are either to update a flow field in a transient flow simulation or to pause to write the particle positions to output files. In the UZ transport application, the GoldSim application performs the time step control, and FEHM is called at each time step to compute the transport. The method is computationally efficient: simulations of several million particles are practical on conventional workstations without employing parallel processing hardware. Like most particle-tracking techniques, the algorithm would parallelize naturally, as the fate of each particle can be computed independently of the movement of other particles. The ability to track the movement of a large number of particles allows large-scale transport simulations to be carried out. The concept of particle residence times in each cell has been extended to relax the assumption of pure plug flow without dispersion. Transfer functions have been developed to simulate dispersion, equilibrium sorption, matrix diffusion, and colloid-facilitated transport. For accuracy, dispersion coefficients must be small enough that the grid Peclet number remains greater than one throughout the grid. This limitation is not viewed as overly restrictive, because systems with large dispersion coefficients can be simulated accurately on relatively coarse grids with conventional finite-difference or finite-element solutions to the advection-dispersion equation. Thus, this particle-tracking method nicely complements the solute transport capabilities present in existing flow and transport codes such as FEHM. The matrix diffusion transfer function fills a void in our ability to simulate tracer experiments in fractured media using equivalent-continuum approaches. Transport time scales in tracer experiments are such that characteristic diffusional distances into the rock matrix are of the order of millimeters. The transfer function is a logical alternative to the explicit simulation of large concentration gradients in the matrix rock in the vicinity of fractures. Linear, equilibrium sorption is also incorporated into the matrix diffusion transfer function. This feature allows a variety of pollutant transport problems to be solved efficiently. When the assumptions inherent in the linear sorption isotherm are inappropriate, more complex chemical transport models should be employed. The ability to implement the RTTF particle-tracking technique efficiently in a dual-permeability model framework is another advantage. Numerically, the communication between fractures and matrix is treated as one additional connection with a flow rate known from the fluid flow calculation. Rapid transit times through fractures can be duplicated easily using the method, whereas a typical equivalent-continuum representation will overestimate the travel time by using ANL-NBS-HS-000026, Rev 00 78 April 2000 the matrix porosity in the computation. Finite-element solutions to the advection-dispersion equation in dual permeability, though possible, are problematic because the short travel times in fractures impose a very stringent constraint on the time step. Furthermore, our benchmarking exercise showed that significant error is introduced for diffusing species due to the lack of grid resolution in the matrix. In principle, a greater number of grid blocks could be placed in the matrix, but this solution is only practical for relatively small systems. It is not a reasonable alternative for site-scale models consisting of a large number of nodes. This limitation does not exist in the RTTF particle-tracking technique. Like all numerical methods, the residence time/transfer function particle-tracking technique has limitations that must be considered when deciding whether its use is appropriate for a given application. The assumptions of advection-dominated transport and linear, equilibrium sorption have just been discussed. The possibility of grid orientation effects must also be considered. These effects can manifest themselves as an artificial lateral spreading of solute mass. However, for the UZ transport model application, the breakthrough of radionuclide mass at the water table is typically desired, rather than in situ concentration. A slight artificial spreading transverse to the flow direction would not be expected to change the travel time of an individual particle significantly. As unstructured numerical grid generation techniques become more common, the simplicity of implementing the RTTF particle-tracking technique for these grids should be a great benefit. In the UZ, for example, unstructured grids are made to follow the complex stratigraphy present at the site, so that units with contrasting hydrologic properties can be represented. It may be that for these grids, orientation errors are small because the grid is aligned with the hydrostratigraphic units and thus are more likely to be aligned with the flow field. Thus, one of the RTTF particletracking technique’s possible limitations should be minimized. Finally, to judge the suitability of the particle-tracking model for unsaturated, dual-permeability transport, a benchmarking exercise was carried out comparing the particle-tracking and finite element, dual-permeability transport models to the presumably more representative discrete fracture model. The results of this suite of test cases highlighted the strengths and weaknesses of the particle-tracking model. The following conclusions were drawn from this analysis: • Particle-tracking solutions that closely track the discrete fracture model result were obtained under all flow scenarios when matrix diffusion is not included. • For low diffusion, the particle-tracking model reproduced the results of the discrete fracture model adequately for a variety of different flow scenarios. In contrast the dual-permeability transport solution yielded an artificial early arrival of solute for flow systems with significant fracture flow along the entire length of the fracture. For situations in which matrix-dominated flow occurred during a portion of the flow path, the early arrival is masked by the longer transport times within that unit, but the problem still exists. • For higher diffusion, the particle-tracking model deviates from the discrete fracture model. This is not surprising, given the assumption used during the development of the method, whereby the fracture transport with matrix diffusion is conceptualized as essentially decoupled from the matrix flow and transport. Whether the model accurately represents ANL-NBS-HS-000026, Rev 00 79 April 2000 reality depends on the actual behavior of the unsaturated zone, particularly with respect to the degree to which diffusion is capable to homogenize concentrations between the flowing fractures and the surrounding matrix. • Not surprisingly, the dual-permeability transport model performs well in the high-diffusion regime. This result is due to the model formulation, consisting of a single matrix node for each fracture node. As in all such dual-porosity and dual-permeability models, accuracy in the low-diffusion regime is sacrificed in favor of a model that captures the long-time behavior of the system. Given these results, it has been demonstrated that the particle-tracking model can be used in three-dimensional radionuclide transport simulations of the Yucca Mountain unsaturated zone as long as the limits on the model are recognized and parameters are chosen accordingly. Also, given the accuracy of the model without diffusion, supporting particle-tracking model runs in the absence of diffusion should be performed to access the importance of matrix diffusion to the overall conclusions of the performance assessment. In addition, an effort should be undertaken to relax the limiting assumptions that reduce the accuracy of the particle-tracking model when large amounts of diffusion are present. An appropriate gage of the extent of diffusion relative to the block size is a characteristic diffusion time to the center of a block. When the diffusion time is on the same order of magnitude or smaller than the transit time through the fracture, the fracture/matrix system in a dual-permeability model could be treated as an intimately coupled entity for the purposes of computing the particle transit time. This enhancement would do a better job at capturing the behavior in the presence of extensive diffusion without sacrificing the model’s ability to handle more moderate diffusion scenarios. This document may be affected by technical product input information that requires confirmation. Any changes to the document that may occur as a result of completing the confirmation activities will be reflected in subsequent revisions. The status of input information quality may be confirmed by review of the Document Input Reference System database. 8. INPUTS AND REFERENCES 8.1 DOCUMENTS CITED Akindunni, F.F.; Gillham, R.W.; Conant, B.; and Franz, T. 1995. "Modeling of Contaminant Movement Near Pumping Wells: Saturated-Unsaturated Flow with Particle Tracking." Ground Water, 33, (2), 264-274. Worthington, Ohio: National Water Well Association. TIC: 226455. Bear, J.; Tsang, C.F.; and de Marsily, G., eds. 1993. Flow and Contaminant Transport in Fractured Rock. San Diego, California: Academic Press. TIC: 235461. Chiang, C.Y.; Wheeler, M.F.; and Bedient, P.D. 1989. "A Modified Method of Characteristics Technique and Mixed Finite Elements Method for Simulation of Groundwater Solute Transport." Water Resources Research, 25, (7), 1541-1549. Washington, D.C.: American Geophysical Union. TIC: 226452. ANL-NBS-HS-000026, Rev 00 80 April 2000 Cordes, C. and Kinzelbach, W. 1992. "Continuous Groundwater Velocity Fields and Path Lines in Linear, Bilinear, and Trilinear Finite Elements." Water Resources Research, 28, (11), 2903- 2911. Washington, D.C.: American Geophysical Union. TIC: 226449. CRC 1991. CRC Handbook of Chemistry and Physics, Lide, D.R., ed., 1991-1992. 72nd Edition. Boca Raton, Florida: CRC Press. TIC: 3595. CRWMS M&O 2000a. Unsaturated Zone Flow and Transport Model and Abstractions for TSPA-SR Process Model Report. TDR-NBS-HS-000002 REV 00. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.20000320.0400. CRWMS M&O 2000b. Input Transmittal from LANL for Abstraction of Colloid-Facilitated Pu Transport Modeling for TSPA (ANL-NBS-HS-000031). 00145T. Las Vegas, Nevada: CRWMS M&O. URN-0145. CRWMS M&O 2000c. Input Transmittal from LANL for Saturated Zone Transport Methodology and Transport Component Integration (MDL-NBS-HS-000010). 00146T. Las Vegas, Nevada: CRWMS M&O. URN-0146. CRWMS M&O 2000d. Abstraction of Flow Fields for RIP (ID:U0125). ANL-NBS-HS-000023 REV 00. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.20000127.0089. Desbarats, A.J. 1990. "Macrodispersion in Sand-Shale Sequences." Water Resources Research, 26, (1), 153-163. Washington, D.C.: American Geophysical Union. TIC: 236576. DOE (U.S. Department of Energy) 1998. Overview - Viability Assessment of a Repository at Yucca Mountain. DOE/RW-0508. Washington, D.C.: U.S. Department of Energy, Office of Civilian Radioactive Waste Management. ACC: MOL.19981007.0027. Fabriol, R.; Sauty, J.P.; and Ouzounian, G. 1993. "Coupling Geochemistry with a Particle Tracking Transport Model." Journal of Contamination Hydrology, 13, 117-129. New York, New York: Elsevier Science. TIC: 222395. Freeze, R.A. and Cherry, J.A. 1979. Groundwater. Englewood Cliffs, New Jersey: Prentice- Hall. TIC: 217571. Golder Associates 2000. User's Guide GoldSim Graphical Simulation Environment. Version 6.02. Draft #3. Redmond, Washington: Golder Associates Inc. TIC: 247324. Iman, R.L. and Shortencarier, M.J. 1984. A FORTRAN 77 Program and User's Guide for the Generation of Latin Hypercube and Random Samples for Use with Computer Models. NUREG/CR-3624. Albuquerque, New Mexico: Sandia National Laboratories. TIC: 224592. Johnson, J.A.; Ravi, V.; and Rumery, J.K. 1994. "Estimation of Solute Concentrations Using the Pathline Counting Method." Ground Water, 32, (5), 719-726. Worthington, Ohio: National Water Well Association. TIC: 226453. ANL-NBS-HS-000026, Rev 00 81 April 2000 Kinzelbach, W 1988. The Random Walk Method in Pollutant Transport Simulation. Groundwater Flow and Quality Modelling. p. 227-245. Norwell, Massachusetts: D. Reidel Publishing Company. TIC: 224064. Latinopoulos, P. and Katsifarakis, K. 1991. "A Boundary Element and Particle Tracking Model for Advective Transport in Zoned Aquifers." Journal of Hydrology, 124, 159-176. Amsterdam, The Netherlands: Elsevier Science. TIC: 226466. Liu, H.H.; Doughty, C.; and Bodvarsson, G.S. 1998. "An Active Fracture Model for Unsaturated Flow and Transport in Fractured Rocks." Water Resources Research, 34, (10), 2633-2646. Washington, D.C.: American Geophysical Union. TIC: 243012. Lu, N. 1994. "A Semianalytical Method of Path Line Computation for Transient Finite- Difference Groundwater Flow Models." Water Resources Research, 30, (8), 2449-2459. Washington, D.C.: American Geophysical Union. TIC: 211997. Maloszewski, P. and Zuber, A. 1991. "Influence of Matrix Diffusion and Exchange Reactions on Radiocarbon Ages in Fissured Carbonate Aquifers." Water Resources Research, 27, (8), 1937-1945. Washington, D.C.: American Geophysical Union. On Order Library Tracking Number-L1369 Marshall, T.J.; Holmes, J.W.; and Rose, C.W. 1996. Soil Physics, Third Edition. New York, New York: Cambridge University Press. TIC: 246638. Nauman, E.B. and Buffham, B.A. 1983. Mixing in Continuous Flow Systems. New York, New York: John Wiley & Sons. TIC: 247226. Neretnieks, I. 1980. "Diffusion in the Rock Matrix: An Important Factor in Radionuclide Retardation." Journal of Geophysical Research, 85, (B8), 4379-4397. Washington, D.C.: American Geophysical Union. TIC: 221345. Neuman, S.P. 1984. "Adaptive Eulerian-Lagrangian Finite Element Method for Advection- Dispersion." International Journal for Numerical Methods in Engineering, 20, 321-337. New York, New York: John Wiley and Sons. TIC: 226440. Pollack, D.W. 1988. "Semianalytical Computation of Path Lines for Finite-Difference Models." Ground Water, 26, (6), 743-750. Worthington, Ohio: National Water Well Association. TIC: 226464. Pruess, K. 1991. TOUGH2-A General-Purpose Numerical Simulator for Multiphase Fluid and Heat Flow. LBL-29400. Berkeley, California: Lawrence Berkeley Laboratory. ACC: NNA.19940202.0088. ANL-NBS-HS-000026, Rev 00 82 April 2000 Reimus, P.W. 1995. The Use of Synthetic Colloids in Tracer Transport Experiments in Saturated Rock Fractures. Ph.D. Dissertation. LA-13004-T. Los Alamos, New Mexico: Los Alamos National Laboratory. TIC: 240694. Reimus, P.W.; Adams, A.; Haga, M.J.; Humphrey, A.; Callahan, T.; Anghel, I.; and Counce, D. 1999. Results and Interpretation of Hydraulic and Tracer Testing in the Prow Pass Tuff at the C-Holes. Milestone SP32E7M4. Los Alamos, New Mexico: Los Alamos National Laboratory. TIC: 246377. Robinson, B.A. 1994. "A Strategy for Validating a Conceptual Model for Radionuclide Migration in the Saturated Zone Beneath Yucca Mountain." Radioactive Waste Management and Environmental Restoration, 19, (1-3), 73-96. Yverdon, Switzerland: Harwood Academic Publishers. TIC: 222513. Schafer-Perini, A.L. and Wilson, J.L. 1991. "Efficient and Accurate Front Tracking for Two- Dimensional Groundwater Flow Models." Water Resources Research, 27, (7), 1471-1485. Washington, D.C.: American Geophysical Union. TIC: 226448. Scheibe, T.D. and Cole, C.R. 1994. "Non-Gaussian Particle Tracking: Application to Scaling of Transport Processes in Heterogeneous Porous Media." Water Resources Research, 30, (7), 2027-2039. Washington, D.C.: American Geophysical Union. TIC: 226446. Starr, R.C.; Gillham, R.W.; and Sudicky, E.A. 1985. "Experimental Investigation of Solute Transport in Stratified Porous Media, 2. The Reactive Case." Water Resources Research, 21, (7), 1043-1050. Washington, D.C.: American Geophysical Union. TIC: 222358. 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. Tompson, A.F.B. and Gelhar, L.W. 1990. "Numerical Simulation of Solute Transport in Three- Dimensional, Randomly Heterogeneous Porous Media." Water Resources Research, 26, (10), 2541-2562. Washington, D.C.: American Geophysical Union. TIC: 224902. van Genuchten, M.T. 1980. "A Closed-Form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils." Soil Science Society American Journal, 44, (5), 892-898. Madison, Wisconsin: Soil Science Society of America. TIC: 217327. van Genuchten, M.T. 1985. "Convective-Dispersive Transport of Solutes Involved in Sequential First-Order Decay Reactions." Computers in Geosciences, 11, (2), 129-147. New York, New York: Pergamon Press. TIC: 245188. ANL-NBS-HS-000026, Rev 00 83 April 2000 Wolfsberg, A.V. and Freyberg, D.L. 1994. "Efficient Simulation of Single Species and Multispecies Transport in Groundwater with Local Adaptive Grid Refinement." Water Resources Research, 30, (11), 2979-2991. Washington, D.C.: American Geophysical Union. TIC: 226450. Yamashita, R. and Kimura, H. 1990. "Particle-Tracking Technique for Nuclide Decay Chain Transport in Fractured Porous Media." Journal of Nuclear Science and Technology, 27, (11), 1041-1049. Tokyo, Japan: Nihon Genshryoku Gakkai. TIC: 236573. Yeh, G.T. 1990. "A Lagrangian-Eulerian Method with Zoomable Hidden Fine-Mesh Approach to Solving Advective-Dispersion Equations." Water Resources Research, 26, (6), 1133-1144. Washington, D.C.: American Geophysical Union. TIC: 224067. Zheng, C. 1993. "Extension of the Method of Characteristics for Simulation of Solute Transport in Three Dimensions." Ground Water, 31, (3), 456-465. Worthington, Ohio: Water Well Publishing Company. TIC: 226465. Zheng, C. 1994. "Analysis of Particle Tracking Errors Associated with Spatial Discretization." Ground Water, 32, (5), 821-828. Dublin, Ohio: Water Well Publishing Company. TIC: 226463. Zimmerman, R.W.; Chen, G.; Hadgu, T.; and Bodvarsson, G.S. 1993. "A Numerical Dual- Porosity Model with Semianalytical Treatment of Fracture/Matrix Flow." Water Resources Research, 29, (7), 2127-2137. Washington, D.C.: American Geophysical Union. TIC: 236571. Zyvoloski, G.; Dash, Z.; and Kelkar, S. 1992. FEHMN 1.0: Finite Element Heat and Mass Transfer Code. LA-12062-MS, Rev. 1. Los Alamos, New Mexico: Los Alamos National Laboratory. ACC: NNA.19910625.0038. Zyvoloski, G.A.; Robinson, B.A.; Dash, Z.V.; and Trease, L.L. 1997. Summary of Models and Methods for the FEHM Application - A Finite Element Heat- and Mass-Transfer Code. LA- 13307-MS. Los Alamos, New Mexico: Los Alamos National Laboratory. TIC: 235587. 8.2 CODES, STANDARDS, REGULATIONS, AND PROCEDURES AP-3.10Q, Rev. 2, ICN 0. Analysis and Models. Washington, D.C.: DOE OCRWM. ACC: MOL.20000217.0246. AP-SI.1Q, Rev 2, ICN 4. OCRWM Procedure: Software Management. Washington, D.C.: DOE OCRWM. ACC: MOL.20000223.0508. CRWMS M&O 1999a. Improve Documentation and Verification of Radionuclide Decay Model (Rev. 01). ID: U3050, Activity: SPP5190. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19990707.0132. ANL-NBS-HS-000026, Rev 00 84 April 2000 CRWMS M&O 1999b. Enhance Particle Tracking Features for TSPA (Rev. 01). ID: U3050, Activity: SPP5190. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19990707.0134. CRWMS M&O 1999c. Improve Matrix Diffusion Model (Rev. 01). ID: U3050, Activity: SPP5190. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19990707.0137. CRWMS M&O 1999d. Conduct of Performance Assessment. Activity Evaluation, September 30, 1999. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19991028.0092. DOE (U.S. Department of Energy) 2000. Quality Assurance Requirements and Description. DOE/RW-0333P, Rev. 9. Washington, D.C.: U.S. Department of Energy, Office of Civilian Radioactive Waste Management. ACC: MOL.19991028.0012. 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 Dr. J.R. Dyer (DOE/YMSCO) to Dr. D.R. Wilkins (CRWMS M&O), September 3, 1999, OL&RC:SB-1714, with enclosure, "Interim Guidance Pending Issuance of New NRC Regulations for Yucca Mountain (Revision 01)." ACC: MOL.19990910.0079. QAP-2-0, Rev. 5. Conduct of Activities. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19980826.0209. QAP-2-3, Rev. 10. Classification of Permanent Items. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19990316.0006. 8.3 SOFTWARE Los Alamos National Laboratory. Software Code: FEHM V2.10. 1999. V2.10. SUN Ultra Sparc, PC. 10086-2.10-00. Lawrence Berkeley National Laboratory. Software Code: FEHM V2.00. 2000. V2.00. SUN Ultra Sparc. 10031-2.00-00. 8.4 SOURCE DATA, LISTED BY DATA TRACKING NUMBER GS950608312231.008. Moisture Retention Data from Boreholes USW UZ-N27 and UE-25 UZ#16. Submittal date: 06/06/1995. GS960808312231.003. Moisture Retention Data for Samples from Boreholes USW SD-7, USW SD-9, USW SD-12 and UE-25 UZ#16. Submittal date: 08/30/1996. GS980908312242.039. Unsaturated Water Retention Data for Lexan-Sealed Samples from USW SD-6 Measured Using a Centrifuge. Submittal date: 09/22/1998. ANL-NBS-HS-000026, Rev 00 85 April 2000 LA0002PR831231.003. Probabilities from C-Wells Microsphere Data. Submittal date: 02/17/2000. LA0003MCG12213.002. Cumulative Probabilities for Colloid Transport Between One Matrix and Another Calculated from Interpolation of Pore Volume Data from Yucca Mountain Hydrologic (Stratigraphic) Samples. Submittal date: 03/10/2000. LB990501233129.001. Fracture Properties for the UZ Model Grids and Uncalibrated Fracture and Matrix Properties for the UZ Model Layers for AMR U0090, "Analysis of Hydrologic Properties Data". Submittal date: 08/25/1999. ANL-NBS-HS-000026, Rev 00 I-1 April 2000 ATTACHMENT I. SOFTWARE ROUTINE TEST Software routine TEST, Version 1.0 was used to compare number of decayed particles calculated based on the direct discrete approach and the integration approach. The algorithm used in the code is listed below: Input output file names for storing calculation data Input simulation start time and end time Input simulation time interval Input particle injection start time and end time Input particle half life Input number of particles per node Input total number of nodes Calculate number of simulation time steps (nt) and decay rate Calculate integration time interval Loop 1, for time step=1 to nt Calculate number of particles decayed per node using the integration approach in Eq. (2) Calculate total number of particles decayed over all nodes Loop 2, for all particles over a node Calculate number of particles decayed on the node using the discrete approach based on Eq. (1) End loop 2 Calculate total number of particles decayed over all nodes Calculate integration % error compared to discrete formula Record the maximum % error and the corresponding time Output the results at each time step Increase the time by the simulation time interval End loop 1 Output the maximum % error, the corresponding time, and the values calculated based on the discrete formula and the integration approach. End software routine TEST The software routine TEST is attached for reference. The software routine was compiled and run on a SunUltra Sparc work station with SunOS 5.5 C++ compiler (property tag: R431923, Sandia National Laboratories, Albuquerque, New Mexico). To comply with the requirements of AP-SI.1Q on software routines, routine TEST was checked with hand calculations. The test case had 1 node with 10 particles uniformly released during the time period of 0 to 9 years, the half life of the species was 10 years. Calculation results from TEST are attached. The hand calculation results at t=30 years were used to check the corresponding software routine results. ANL-NBS-HS-000026, Rev 00 I-2 April 2000 Table I-1 lists the results of hand calculations and the corresponding results of software routine TEST. A comparison between those results showed good agreement and proved that the software routine TEST was correct. Table I-1. Results From Hand Calculations and Software Routine TEST Time, t=30 Discrete Approach Integration Approach Particle Index Release Time, ti (years) Calculation 1-exp(-.(t-ti)) TEST Calculation {(t1-t2)+[exp(-.t1)-exp(-.t2)]/.}/.t TEST 1 0 0.875000 2 1 0.866028 3 2 0.856412 4 3 0.846106 5 4 0.835061 6 5 0.823223 7 6 0.810535 8 7 0.796936 9 8 0.782362 10 9 0.766741 t1=t-t1=30-0=30 t2=t-t10=30-9=21 .=0.06931472 .t=9/10=0.9 (t2-t1)=30-21=9 exp(-.t1)=0.12500 exp(-.t2)=0.23326 Total: sum=sum+0.5 8.758409 8.81819 8.76456 8.76463 int(sum) 8 8 8 8 ANL-NBS-HS-000026, Rev 00 I-3 April 2000 List of Software Routine TEST Output Results for the Test Problem input simulation start time and end time 0 40 input simulation time interval dt 10 The particle injection start time and end time 0 9 The half life 10 number of particles released per node 10 Total number of nodes 1 At time=0 Integration Decayed: 0 Discrete Decayed: 0 Integration Percent error: 0 At time=10 Integration Decayed: 3 Discrete Decayed: 3 Integration Percent error: 0 At time=20 Integration Decayed: 7 Discrete Decayed: 7 Integration Percent error: 0 At time=30 Integration Decayed: 8 Discrete Decayed: 8 Integration Percent error: 0 Over Results: Discrete == Integration ANL-NBS-HS-000026, Rev 00 I-4 April 2000 List of Software Routine TEST //start of code TEST #include #include #include #include #include main() { int i,j,k,l,ii,nt; char outfile[60],errfile[60]; double t1,t2,tt,tot,half,t,m,dt,time,ti2,ti1,np; double tot1,pert=0.,dtime,endtime,maxerror=0.0; double tmptime=0.,tmptot=0.,tmptot2=0.; cout<<"Output file name for storing calculation results"<>outfile; cout<<"Output file name for storing %error data"<>errfile; ofstream outFile (outfile, ios::out); if(!outFile) { cerr<<"Can not open file:"<>time>>endtime; cout<<"input time interval dt"<>dtime; cout<<"The injection start time and end time "<>t1>>t2; cout<<"The half life"<>half; cout<<"number of particles released per node"<>np; cout<<"Total number of nodes"<>m; half=0.6931472/half; //degration rate ANL-NBS-HS-000026, Rev 00 I-5 April 2000 nt=(endtime-time)/dtime; dt=(t2-t1)/np; //calculate # of particles decayed using the integration method for(ii=1;ii<=nt;ii++) { cout<<"At time="<=0.) tot +=(1-exp(-half*(time-t))); } tot=tot*m; tot=int(tot+0.5); cout<<"Discrete: tot="< maxerror) { maxerror=pert; tmptime=time; tmptot=tot1; ANL-NBS-HS-000026, Rev 00 I-6 April 2000 tmptot2=tot; } if(tot1 >=np*m) { cerr<<"Integration larger than the original"<>tracer solution started on time step 1 mptr read from optional input file: inp.ptrk Species: 1 at time: 9.993155373032168E-02 rest: 0.000000000000000 Species: 2 at time: 9.993155373032168E-02 rest: 0.000000000000000 Species: 3 at time: 9.993155373032168E-02 rest: 0.000000000000000 Species: 4 at time: 9.993155373032168E-02 rest: 0.000000000000000 Species: 5 at time: 9.993155373032168E-02 rest: 0.000000000000000 ********************************************************************* Time Step 1 Timing Information Years Days Step Size (Days) 0.999316E-01 0.365000E+02 0.365000E+02 Cpu Sec for Time Step = 0.4000E-01 Current Total = 0.4000E-01 Equation Performance Number of N-R Iterations: 1 Avg # of Linear Equation Solver Iterations: 2.0 Number of Active Nodes: 402. Total Number of Iterations, N-R: 1 , Solver: 2 Largest Residuals EQ1 R= 0.5335E-13 node= 240 x=0.1900 y= 0.000 z= 1.000 EQ2 R= 0.3113E-14 node= 54 x=0.2650 y= 1.000 z= 1.000 Node Equation 1 Residual Equation 2 Residual 201 -0.308559E-14 -0.364009E-15 Nodal Information (Water) source/sink source/sink Node p(MPa) e(MJ) l sat temp(c) (kg/s) (MJ/s) 201 1.000 0.11 1.000 25.000 0.500E-09 0.528E-10 Global Mass & Energy Balances Total mass in system at this time: 0.299168E+03 kg Total mass of steam in system at this time: 0.000000E+00 kg Total enthalpy in system at this time: 0.750546E+02 MJ Water discharge this time step: 0.315360E-02 kg (0.100000E-08 kg/s) ANL-NBS-HS-000026, REV 00 II-6 April 2000 Water input this time step: 0.315360E-02 kg (0.100000E-08 kg/s) Total water discharge: 0.315360E-02 kg (0.100000E-08 kg/s) Total water input: 0.315360E-02 kg (0.100000E-08 kg/s) Enthalpy discharge this time step: 0.333149E-03 MJ (0.105641E-09 MJ/s) Enthalpy input this time step: 0.333149E-03 MJ (0.105641E-09 MJ/s) Total enthalpy discharge: 0.333149E-03 MJ (0.105641E-09 MJ/s) Total enthalpy input: 0.333149E-03 MJ (0.105641E-09 MJ/s) Net kg water discharge (total out-total in): 0.464224E-11 Net MJ discharge (total out-total in): 0.487107E-12 Conservation Errors: 0.111155E-08 (mass), 0.139788E-09 (energy) ************************* Particle Tracking ==> Species: 1 Number Having Entered System: 10 Number Currently In System : 10 Number Having Left System : 0 Number Having Decayed : 0 Node Concentration # of Particles 201 0.000000E+00 0 ************************* Particle Tracking ==> Species: 2 Number Having Entered System: 0 Number Currently In System : 0 Number Having Left System : 0 Number Having Decayed : 0 Node Concentration # of Particles 201 0.000000E+00 0 ************************* Particle Tracking ==> Species: 3 Number Having Entered System: 0 Number Currently In System : 0 Number Having Left System : 0 Number Having Decayed : 0 Node Concentration # of Particles 201 0.000000E+00 0 ************************* Particle Tracking ==> Species: 4 Number Having Entered System: 0 Number Currently In System : 0 Number Having Left System : 0 Number Having Decayed : 0 Node Concentration # of Particles 201 0.000000E+00 0 ************************* Particle Tracking ==> Species: 5 Number Having Entered System: 10 Number Currently In System : 10 Number Having Left System : 0 Number Having Decayed : 0 Node Concentration # of Particles 201 0.000000E+00 0 Species: 1 ANL-NBS-HS-000026, REV 00 II-7 April 2000 at time: 0.1002053388090349 rest: 0.000000000000000 Species: 2 at time: 0.1002053388090349 rest: 0.000000000000000 Species: 3 at time: 0.1002053388090349 rest: 0.000000000000000 Species: 4 at time: 0.1002053388090349 rest: 0.000000000000000 Species: 5 at time: 0.1002053388090349 rest: 0.000000000000000 ********************************************************************* Time Step 36 Timing Information Years Days Step Size (Days) 0.990010E+04 0.361601E+07 0.109575E+06 Heat and Mass Solution Disabled ************************* Particle Tracking ==> Species: 1 Number Having Entered System: 500000 Number Currently In System : 0 Number Having Left System : 0 Number Having Decayed : 500000 Node Concentration # of Particles 201 0.000000E+00 0 ************************* Particle Tracking ==> Species: 2 Number Having Entered System: 500000 Number Currently In System : 433657 Number Having Left System : 56347 Number Having Decayed : 9996 Node Concentration # of Particles 201 0.410205E+04 1534 ************************* Particle Tracking ==> Species: 3 Number Having Entered System: 9996 Number Currently In System : 8598 Number Having Left System : 1390 Number Having Decayed : 8 Node Concentration # of Particles 201 0.125682E+03 47 ************************* Particle Tracking ==> Species: 4 Number Having Entered System: 8 Number Currently In System : 5 Number Having Left System : 3 Number Having Decayed : 0 Node Concentration # of Particles 201 0.000000E+00 0 ************************* ANL-NBS-HS-000026, REV 00 II-8 April 2000 Particle Tracking ==> Species: 5 Number Having Entered System: 500000 Number Currently In System : 441782 Number Having Left System : 58218 Number Having Decayed : 0 Node Concentration # of Particles 201 0.215799E+04 807 Species: 1 at time: 10200.10403832991 rest: 0.000000000000000 Species: 2 at time: 10200.10403832991 rest: 21456.00000000000 Species: 3 at time: 10200.10403832991 rest: 558.0000000000000 Species: 4 at time: 10200.10403832991 rest: 2.000000000000000 Species: 5 at time: 10200.10403832991 rest: 21762.00000000000 ********************************************************************* Time Step 70 Timing Information Years Days Step Size (Days) 0.199863E+05 0.730000E+07 0.680120E+05 Heat and Mass Solution Disabled ************************* Particle Tracking ==> Species: 1 Number Having Entered System: 500000 Number Currently In System : 0 Number Having Left System : 0 Number Having Decayed : 500000 Node Concentration # of Particles 201 0.000000E+00 0 ************************* Particle Tracking ==> Species: 2 Number Having Entered System: 500000 Number Currently In System : 0 Number Having Left System : 486844 Number Having Decayed : 13156 Node Concentration # of Particles 201 0.000000E+00 0 ************************* Particle Tracking ==> Species: 3 Number Having Entered System: 13156 Number Currently In System : 0 Number Having Left System : 13082 Number Having Decayed : 74 Node Concentration # of Particles 201 0.000000E+00 0 ************************* Particle Tracking ==> Species: 4 ANL-NBS-HS-000026, REV 00 II-9 April 2000 Number Having Entered System: 74 Number Currently In System : 0 Number Having Left System : 64 Number Having Decayed : 10 Node Concentration # of Particles 201 0.000000E+00 0 ************************* Particle Tracking ==> Species: 5 Number Having Entered System: 500000 Number Currently In System : 0 Number Having Left System : 500000 Number Having Decayed : 0 Node Concentration # of Particles 201 0.000000E+00 0 simulation ended: days 7.300E+06 timesteps 70 total N-R iterations = 2 total solver iterations = 4 total code time(timesteps) = 1061.690039 ****---------------------------------------------------------**** **** This program for **** **** Finite Element Heat and Mass Transfer in porous media **** ****---------------------------------------------------------**** **** Version : FEHM V2.00 99-05-03 **** **** End Date : 06/14/1999 **** **** Time : 09:47:58 **** ****---------------------------------------------------------**** ANL-NBS-HS-000026, REV 00 II-10 April 2000 Routine BKPM was run to extract particle breakthrough information at the outflow boundary. The default output files from BKPM were spec_1, spec_2, spec_3, and spec_4 for speciese 1, 2, 3, and 4, respectively. Each output file contains the following information: Time (year), mass flux (particles/year), fractions of particles left during the current time step, and normalized cummulative breakthrough data. The example below demonstrates the use of BKPM and the sample output results. picard< 41 > % bkpm Total number of particles: 500000 Change the Value by Typing in a New Number or 1 1 Calculate mass flow rate->1 or concentration->2 1 Cal. Normalized Cumulative Breakthrough Data (Y/N)? y Output Data Stored in File: spec_1 Output Data Stored in File: spec_2 Output Data Stored in File: spec_3 Output Data Stored in File: spec_4 Output Data Stored in File: spec_5 Finished BKPM extracted results for the 5 species are listed below. File spec_1: time (year) mass flux fraction_left normalized cumulative 0.499658E-01 0. 0.000000 0.000000 0.495010E+04 0. 0.000000 0.000000 0.149432E+05 0. 0.000000 0.000000 File spec_2: time (year) mass flux fraction_left normalized cumulative 0.499658E-01 0. 0.000000 0.000000 0.495010E+04 0.569162E+01 0.112694 0.112694 0.149432E+05 0.426818E+02 0.860994 0.973688 ANL-NBS-HS-000026, REV 00 II-11 April 2000 File spec_3: time (year) mass flux fraction_left normalized cumulative 0.499658E-01 0. 0.000000 0.000000 0.495010E+04 0.140404E+00 0.002780 0.002780 0.149432E+05 0.115921E+02 0.023384 0.026164 File spec_4: time (year) mass flux fraction_left normalized cumulative 0.499658E-01 0. 0.000000 0.000000 0.495010E+04 0.303030E-03 0.000006 0.000006 0.149432E+05 0.604787E-02 0.000122 0.000128 File spec_5: time (year) mass flux fraction_left normalized cumulative 0.499658E-01 0. 0.000000 0.000000 0.495010E+04 0.588061E+01 0.116436 0.116436 0.149432E+05 0.438006E+02 0.883564 1.000000 Hand calculations were performed to check the results from BKPM. The total mass input to the system is totmass=1 represented by a total of 500,000 (total) particles. The mass flow rate, fractions left, and normalized cumulative breakthrough data were calculated based on Equations (II-1), (II-2), (II-3), and (II-4), respectively. At time tc=0.999316E-1 years, tp=0.0 years, tc-tp=0.999316E-1 years . 5 , 4 , 3 , 2 , 1 , 0 500000 0 _ _ . 5 , 4 , 3 , 2 , 1 , 0 500000 0 0 _ _ . 5 , 4 , 3 , 2 , 1 , 0 10 999316 . 0 0 0 _ _ : 1 = = = = = - = = = × - = - i i cumulative normalized i i left fractions i i flux mass i Species ANL-NBS-HS-000026, REV 00 II-12 April 2000 At time tc=0.990010E4 years, tp=0.999316E-1 years, tc-tp=0.990000E4 years 116436 . 0 500000 58218 5 _ _ 116436 . 0 500000 0 58218 5 _ _ 10 117612 . 0 500000 1 00 . 9900 0 58218 5 _ _ 5 : 000006 . 0 500000 3 4 _ _ 000006 . 0 500000 0 3 4 _ _ 10 606061 . 0 500000 1 00 . 9900 0 3 4 _ _ 4 : 00278 . 0 500000 1390 3 _ _ 00278 . 0 500000 0 1390 3 _ _ 10 280808 . 0 500000 1 00 . 9900 0 1390 3 _ _ 3 : 112694 . 0 500000 56347 2 _ _ 112694 . 0 500000 0 56347 2 _ _ 10 113832 . 0 500000 1 00 . 9900 0 56347 2 _ _ 2 : 0 500000 0 1 _ _ 0 500000 0 1 _ _ 0 00 . 9900 0 0 1 _ _ 1 : 4 11 6 4 = = = - = × = × - = = = = - = × = × - = = = = - = × = × - = = = = - = × = × - = = = = = = - = - - - - cumulative normalized left fractions flux mass Species cumulative normalized left fractions flux mass Species cumulative normalized left fractions flux mass Species cumulative normalized left fractions flux mass Species cumulative normalized left fractions flux mass Species ANL-NBS-HS-000026, REV 00 II-13 April 2000 At time tc=0.199863E5 years, tp=0.990010E4 years, tc-tp=0.100862E5 years 0 . 1 500000 500000 5 _ _ 883564 . 0 500000 58218 500000 5 _ _ 10 876013 . 0 500000 1 2 . 10086 58218 500000 5 _ _ 5 : 000128 . 0 500000 64 4 _ _ 000122 . 0 500000 3 64 4 _ _ 10 120957 . 0 500000 1 2 . 10086 3 64 4 _ _ 4 : 026164 . 0 500000 13082 3 _ _ 023384 . 0 500000 1390 13082 3 _ _ 10 231842 . 0 500000 1 2 . 10086 1390 13082 3 _ _ 3 : 973688 . 0 500000 486844 2 _ _ 860994 . 0 500000 56347 486844 2 _ _ 10 853636 . 0 500000 1 2 . 10086 56347 486844 2 _ _ 2 : 0 500000 0 1 _ _ 0 500000 0 1 _ _ 0 00 . 9900 0 0 1 _ _ 1 : 4 7 5 4 = = = - = × = × - = = = = - = × = × - = = = = - = × = × - = = = = - = × = × - = = = = = = - = - - - - cumulative normalized left fractions flux mass Species cumulative normalized left fractions flux mass Species cumulative normalized left fractions flux mass Species cumulative normalized left fractions flux mass Species cumulative normalized left fractions flux mass Species Results from the hand calculations are listed in Table II-1 for comparison with the results from routine BKPM. The comparison results showed that software routine BKPM performed as designed and extracted correct results from the particle output file. ANL-NBS-HS-000026, REV 00 II-14 April 2000 Table II-1. Comparison of BKPM and Hand Calculation Results Time (year) Species Type of Cal. Mass Flux (particle/year) Fractions Left Normalized Cumulative BKPM 0.0 0.0 0.0 1 hand cal. 0.0 0.0 0.0 BKPM 0.0 0.0 0.0 2 hand cal. 0.0 0.0 0.0 BKPM 0.0 0.0 0.0 3 hand cal. 0.0 0.0 0.0 BKPM 0.0 0.0 0.0 4 hand cal. 0.0 0.0 0.0 BKPM 0.0 0.0 0.0 0.999316E-1 5 hand cal. 0.0 0.0 0.0 BKPM 0.0 0.0 0.0 1 hand cal. 0.0 0.0 0.0 BKPM 5.69162 0.112694 0.112694 2 hand cal. 5.69162 0.112694 0.112694 BKPM 1.40404E-1 0.00278 0.00278 3 hand cal. 1.40404E-1 0.00278 0.00278 BKPM 3.0303E-4 0.000006 0.000006 4 hand cal. 3.0303E-4 0.000006 0.000006 BKPM 5.88061 0.116436 0.116436 0.99001E4 5 hand cal. 5.88061 0.116436 0.116436 BKPM 0.0 0.0 0.0 1 hand cal. 0.0 0.0 0.0 BKPM 4.26818E+1 0.860994 0.97368 2 hand cal. 4.26818E+1 0.860994 0.973688 BKPM 1.15921 0.023384 0.026164 3 hand cal. 1.15921 0.023384 0.026164 BKPM 6.04787E-3 0.000122 0.000128 4 hand cal. 6.04787E-3 0.000122 0.000128 BKPM 4.38006E+1 0.883564 1.0 0.199863E4 5 hand cal. 4.38006E+1 0.883564 1.0 ANL-NBS-HS-000026, REV 00 II-15 April 2000 Listing of FEHM post processing software routine BKPM program bkpm C********************************************************************** C bkpm is a software routine used to extract breakthrogh cureves * C at the out flow boundary from the FEHM particle output file *.out * C The routine automatically read in FEHM control file: fehm.files * C and search for the particle output file, it then open the output * C and process the data to get the particle mass flow rate, the * C cumulative fraction of particle flow out of the system, and the * C number of particles flow out of the system at each time step. * C * C Chunhong Li, 6/14/1999 * C********************************************************************** character*80 line,prx*5,cidx*5,title*31,cbyn*1 character*60 outfile,fout*10,dfile*30,apix*10,ctmp real*8 time,time_old(30),reptime double precision conv,conc integer entered,cur_in,left,tot_entered,decayed integer left_old(30) character*6 spacer2 parameter(rate_const = 3.239e-7) data time_old/30*0.0/left_old/30*0/ conv=1.0 ifehm=2 apix='.btp' prx='spec_' title='Particle Tracking ==> Species: ' dfile='fehmn.files' open(ifehm,file=dfile,status='old',iostat=iofehm) if(iofehm.eq.0)then do ip=1,4 read(ifehm,*)outfile end do else write(*,*) 'Name of output file: ' read(*,'(a60)') outfile endif call prefix(outfile,leng) open(7,file=outfile) reptime=0. c determine number of particles in problem 1 read(7,'(a80)',end=99) line if (line(16:22).eq.'Entered') read(line,72) spacer2,tot_entered 72 format(1a30,i11) goto 1 ANL-NBS-HS-000026, REV 00 II-16 April 2000 99 write(*,*) 'Total number of particles: ',tot_entered total=tot_entered write(6,*)'Change the Value by Typing in a New Number or 1' read(5,*)indx if(indx.ne.1)total=indx write(6,*)'Calculate mass flow rate->1 or concentration->2' read(5,*)indx if(indx.eq.2)then write(6,*)'Conversion factor used for concentration calculation' read(5,*)conv endif write(6,*)'Cal. Normalized Cumulative Breakthrough Data (Y/N)?' read(5,*)cbyn rewind(7) write(6,*)'' ist=0 2 read(7,'(a80)',end=999) line if (line(12:16).eq.'Years') then read(7,*) time ist=ist+1 iout=20 idx=0 else if(line(24:54).eq.title)then read(line,59)ctmp,cidx 59 format(A54,a5) iout=iout+1 idx=idx+1 if(ist.eq.1)then do i=1,5 if(cidx(i:i).ne. ' ')jj=i enddo fout=prx//cidx(jj:5) open(iout,file=fout) write(6,*)'Output Data Stored in File: ',fout endif read(7,'(a80)') line if (line(16:22).eq.'Entered') then read(line,72) spacer2,entered read(7,'(a80)') line read(line,72) spacer2,cur_in read(7,'(a80)') line read(line,72) spacer2,left read(7,'(a80)') line read(line,72) spacer2,decayed left_this = left - left_old(idx) delta_time = time - time_old(idx) if(delta_time.ne.0.)then reptime=time_old(idx)+0.5*delta_time conc = real(left_this) cnorm=conc/total conc = conv*conc/delta_time left_old(idx) = left time_old(idx) = time ANL-NBS-HS-000026, REV 00 II-17 April 2000 if(cbyn.eq.'y')then write(iout,100) reptime,conc,cnorm,left/total else write(iout,110) reptime,conc endif endif endif endif endif goto 2 999 write(6,*)'' write(6,*)'Finished' 100 format(1x,E13.6,3x,E13.6,2(3x,f9.6)) 110 format(1x,E13.6,3x,E13.6) end subroutine prefix(fname,length) character*60 fname do length=1,60 if(fname(length:length).eq.'.')goto 100 end do 100 length=length-1 return end ANL-NBS-HS-000026, Rev 00 III-1 April 2000 ATTACHMENT III. SOFTWARE ROUTINE CHAIN Software routine CHAIN, Version 1.0, written in FORTRAN IV was acquired from U.S. Salinity Laboratory at 450 Big Springs Rd., Riverside, CA 92507. The routine was used to calculate analytical solutions for convective-dispersive transport of solutes involved in sequential firstorder decay reactions. The routine version number and subsequent changes is listed below. Routine Version number Changes made CHAIN May 1982. Code issuer did not have a version number. None The theory and applications of CHAIN were published by M. Th. van Genuchten (van Genuchten 1985). The paper documented in detail the theory, assumptions, and boundary conditions used in deriving the analytical solutions and the code CHAIN developed for calculating the analytical solutions. Two examples on the applications of CHAIN were given in the paper. The executable version of CHAIN was installed on a DELL OptiPlex GX1 PC with windows NT 4.0 operating system at Sandia National Laboratories in Albuquerque, New Mexico. The PC has 128 MB memory, 9 GB hard disk drive, and a CPU speed of 300 MHz. No modification was made to CHAIN after it was received or installed. Hand calculations were carried out to check against the results from CHAIN. A 4 species chain decay model in a semi-infinite system with 0 initial concentrations was used in the comparison. Only species 1 was released at x=0 from time 0 to t0 years at a constant concentration of B1. ( ) ( ) 4 , 1 0 , 0 , 2 , 0 ) , 0 ( , 0 0 , ) , 0 ( 0 0 1 1 = = = 8 . . = = . . . = i t t x c i t c t t t t B t c i i f p p The solutions for this system are (van Genuchten, 1985): .. .. . - - - = = 0 0 * 0 * 0 * ), , ( ) exp( ) , ( 0 ), , ( t t t t x c t t x c t t t x c c i i i i i f p . (Eq. III-1) where .i=Nb µi +.i, µi is the rate constanstant of species i for first order decay, .i is the source term release rate for species i, and Nb equals 1 or 0. For Nb=1, Bateman equations are used for input boundary conditions, for Nb=0, Bateman equations are not used, parameter Bi are read in, x is the distance from the origin of source release, and ci * is part of the analytical solution for species i and was given in appendix 1 in van Genuchten’s paper (van Genuchten, 1985) with ANL-NBS-HS-000026, Rev 00 III-2 April 2000 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ) ( ) ( ) ( ) ( ) ( ) ( 12 1 12 124 123 110 210 12 2 12 1 3 23 1 23 234 123 210 310 23 2 23 1 3 13 1 13 134 123 13 110 310 2 13 1 3 34 1 34 234 134 310 410 34 2 34 1 3 24 1 24 234 124 24 210 410 2 24 1 3 14 1 14 124 134 110 410 14 2 14 1 3 * 4 23 1 23 23 210 310 23 1 1 13 1 13 110 310 13 13 1 1 12 1 12 12 110 210 12 1 1 * 3 12 1 12 110 210 12 1 1 1 * 2 110 1 * 1 µ . µ . µ . µ . µ . µ . µ . µ . µ . µ . µ - + - + - + - + - - - + - + - + - - - + - + - = - - - + - + - + - - - = - + - = = R G G F F S R B A R G G F F S R B A R G G S F F R B A R G G F F S R B A R G G S F F R B A R G G F F S R B A c R S F F R B A R F F S R B A R S F F R B A c R F F S B R c F B c (Eq. III-2) With ( ) ( ) ( ) ( ) ( ) ( ) [ ]2 1 2 2 / 1 2 / 1 3 2 1 3 2 1 3 234 3 2 3 2 2 123 2 1 2 1 1 4 0 , 0 , 2 2 exp 2 1 2 2 exp 2 1 exp ijk i i ij ij j ijk iji jji ij i i i i ijk ijk ji k k ik j j kj i i ijk j i ij j i ij a DR v w k R k a wehre F F S t DR wt x R erfc D x w v t DR wt x R erfc D x w v t a F R R R A G R R A G R R A R R R R R R G R R R - + = .. .. . > = = - = .. .. . .. .. . . .. . . .. . + .. . .. . + + . .. . . .. . - .. . .. . - - = = = = + + = - = - = µ µ . µ µ µ µ µ µ µ µ µ µ µ µ µ In the above equations, v is average pore water velocity, and Ri is the ratardation factor of the i th species. Hand calculations were carried out for a 4 member sequential decay chain system to check the correctness of the CHIAN results under the same conditions. The half lives for species 1 through 4 were 1,000, 3,000, 20,000, and 4,000 years, respectively. Species 1 was released at x=0 from time 0 to 5,000 years at constant concentration of B1=1.0. Solute concentrations of each species were calculated at a distance of x=1.0 m. Other parameters used in the calculations are listed below. Nb=0, B1=1, B2= B3=B4=0, .i =0 (i=1,4), x=1.0 m, .i =Nb µi +.i,=0 (i=1,2,3,and 4) µ1=6.9315E-4/year, µ2=2.3105E-4/year, µ3=3.4657E-5/year, µ4=1.7329E-4/year, v=1.0519E-4 m/year, D=5.2595E-7 m2/year, R1=1.0, R2=1.0, R3=1.001, and R4=1.001. (Eq. III-3) ANL-NBS-HS-000026, Rev 00 III-3 April 2000 Software routine CHAIN was used to calculate the solutes breakthrough curves at x=1.0 m at a time interval of 300 years from 0 to 20,000 years. The CHAIN input and output files are attached with this document for referencing. Hand calculations were then carried out to randomly check the results from CHAIN. The following parameters were calculated based on the input parameters given above. R12=0, R13=-1.E-3, R14=-1.E-3, R23=-1.E-3, R24=-1.E-3, R34=0 µ12=4.6210E-4, µ13=6.5849E-4, µ14=5.1986E-4, µ23=1.9639E-4, µ24=0.5776E-4, µ34=-1.3863E-4 G123=4.6212E-7 G124=4.6212E-7 G134=-1.3878E-7 G234=-1.3878E-7 A1=0.3466 A3=5.5560E-12 The partial solutions ci * can be rewritten as ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 110 210 12 46 210 310 23 45 13 110 310 44 310 410 34 43 24 210 410 42 110 410 14 41 * 4 23 210 310 33 110 310 13 32 12 110 210 31 * 3 110 210 12 21 * 2 110 1 * 1 ) ( F F S C F F S C S F F C F F S C S F F C F F S C c S F F C F F S C S F F C c F F S C c F B c + - + + - + - - + + - + - - + + - = - - + + - + - - = + - = = (Eq. III-4) With To compare against the solutions from routine CHAIN at x=1.0 m, we needed to calculate functions Fijk at a given time t and solve for ci. In van Genutchen’s paper, approximation was made to calculate the products of exp(A)erfc(B) in Fijk using the following formula. 0000 . 0 ) ( 4412 . 0 ) ( 1316 . 0 ) ( 0000 . 0 ) ( 5000 . 1 ) ( 1666 . 0 ) ( 7649 . 1 5264 . 0 0000 . 0 5000 . 1 12 1 12 124 123 2 12 1 3 46 23 1 23 234 123 2 23 1 3 45 13 1 13 134 123 2 13 1 3 44 34 1 34 234 134 2 34 1 3 43 24 1 24 234 124 2 24 1 3 42 14 1 14 124 134 2 14 1 3 41 23 1 23 23 1 1 33 13 1 13 13 1 1 32 12 1 12 12 1 1 31 12 1 12 1 1 21 = - = - = - = - = - = = - = - = - = - = - = = - = = - = = - = - = - = µ . µ . µ . µ . µ . µ . µ . µ . µ . µ . µ R G G R B A C R G G R B A C R G G R B A C R G G R B A C R G G R B A C R G G R B A C R R B A C R R B A C R R B A C R B R C ANL-NBS-HS-000026, Rev 00 III-4 April 2000 ). 170 ( , 0 , , 0 , 0 ) ( ) exp( ) ( ) exp( ) exp( 2 ) ( ) exp( 0 )))))) 1 /( 5 . 2 /( 2 /( 5 . 1 /( 1 /( 5 . 0 ( ) exp( 1 ) ( ) exp( 5 . 2 061405429 . 1 , 453152027 . 1 421413741 . 1 , 284496736 . 0 254829592 . 0 , 3275911 . 0 1 1 ) )( exp( ) ( ) exp( 5 . 2 0 2 2 5 4 3 2 1 5 5 4 4 3 3 2 2 1 2 = .. .. . .. .. . > > - = > ˜ - - = < + + + + + + - ˜ > = - = = - = = + = + + + + - ˜ = = M B M B A or B M A for B erfc A and B erfc A A B erfc A B For B B B B B B B A B erfc A B For B where B A B erfc A B For p a a a a a t t a t a t a t a t a and Sij=0, for abs(1000xRij/Ri)<1 (van Genutchen, 1985) Since in our case, the absolute ratio of abs(1000xRij/Ri)<1, Sij=0 in the calculations. Hand calculations were carried out to check the results of CHAIN at selected time of t=3000, 8100, 12300, and 17400 years. ANL-NBS-HS-000026, Rev 00 III-5 April 2000 At time t=3000 years: Substituting parameters into the Equation (3), we calculated functions Fijk at the specified time. 0 10 4935 . 1 ) 0795 . 0 3217 . 1 ( ) 10 0519 . 1 10 1210 . 2 exp( ) 0795 . 0 6803 . 0 ( ) 10 0519 . 1 10 7206 . 1 exp( 5 . 0 10 2572 . 2 ) 0795 . 0 3176 . 1 ( ) 10 0519 . 1 10 1073 . 2 exp( ) 0795 . 0 6844 . 0 ( ) 10 0519 . 1 10 4634 . 3 exp( 5 . 0 10 4512 . 1 ) 0794 . 0 3224 . 1 ( ) 10 0519 . 1 10 1267 . 2 exp( ) 0795 . 0 6776 . 0 ( ) 10 0519 . 1 10 2857 . 2 exp( 5 . 0 10 6631 . 3 ) 0794 . 0 3357 . 1 ( ) 10 0519 . 1 10 1710 . 2 exp( ) 0795 . 0 6643 . 0 ( ) 10 0519 . 1 10 7170 . 6 exp( 5 . 0 34 24 14 23 13 12 34 6 4 6 6 410 34 6 4 6 6 310 34 6 4 6 6 210 35 6 4 6 6 110 = = = = = = × = . .. . . .. . × × + × × - = × = . .. . . .. . × × + × × - = × = . .. . . .. . × × + × × - = × = . .. . . .. . × × + × × - = - - - - - - - - - - - - - - - - - - - - S S S S S S erfc erfc F erfc erfc F erfc erfc F erfc erfc F As t<5,000 years, ci=ci *, From Equation (1), the normalized solute concentrations at t=3,000 years are 35 34 34 35 34 34 34 34 35 310 210 45 110 310 44 210 410 42 410 110 41 4 35 34 34 34 34 34 35 210 310 33 310 110 32 3 34 34 35 210 110 21 2 35 35 110 1 1 10 3111 . 2 ) 10 2572 . 2 10 4512 . 1 ( 4412 . 0 ) 10 6631 . 3 10 2572 . 2 ( 1316 . 0 ) 10 4512 . 1 10 4935 . 1 ( 5 . 1 ) 10 4935 . 1 10 6631 . 3 ( 1666 . 0 ) ( ) ( ) ( ) ( 10 2714 . 4 10 06 . 8 7649 . 1 10 8909 . 1 5264 . 0 ) 10 4512 . 1 10 2572 . 2 ( 7649 . 1 ) 10 2572 . 2 10 6631 . 3 ( 5264 . 0 ) ( ) ( 10 6273 . 1 ) 10 4512 . 1 10 6631 . 3 ( 5 . 1 ) ( 10 6631 . 3 10 6631 . 3 1 - - - - - - - - - - - - - - - - - - - - - × = × - × - × - × - × - × - × - × - = - + - + - + - = × = × × + × × - = × - × + × - × = - + - = × = × - × - = - = × = × × = = F F C F F C F F C F F C c F F C F F C c F F C c F B c ANL-NBS-HS-000026, Rev 00 III-6 April 2000 At time t=8100 years 0 10 5318 . 1 ) 1306 . 0 8670 . 1 ( ) 10 0519 . 1 10 1210 . 2 exp( ) 1306 . 0 1350 . 0 ( ) 10 0519 . 1 10 7206 . 1 exp( 5 . 0 10 4908 . 4 ) 1306 . 0 8558 . 1 ( ) 10 0519 . 1 10 1073 . 2 exp( ) 1306 . 0 1462 . 0 ( ) 10 0519 . 1 10 4634 . 3 exp( 5 . 0 10 9902 . 9 ) 1305 . 0 8706 . 1 ( ) 10 0519 . 1 10 1267 . 2 exp( ) 1305 . 0 1294 . 0 ( ) 10 0519 . 1 10 2857 . 2 exp( 5 . 0 10 8140 . 2 ) 1305 . 0 9065 . 1 ( ) 10 0519 . 1 10 1710 . 2 exp( ) 1305 . 0 0936 . 0 ( ) 10 0519 . 1 10 7170 . 6 exp( 5 . 0 34 24 14 23 13 12 2 6 4 6 6 410 2 6 4 6 6 310 3 6 4 6 6 210 4 6 4 6 6 110 = = = = = = × = . .. . . .. . × × + × × - = × = . .. . . .. . × × + × × - = × = . .. . . .. . × × + × × - = × = . .. . . .. . × × + × × - = - - - - - - - - - - - - - - - - - - - - S S S S S S erfc erfc F erfc erfc F erfc erfc F erfc erfc F For (t-t0)=8100-5000=3100 years 0 10 4818 . 1 ) 0808 . 0 3324 . 1 ( ) 10 0519 . 1 10 1210 . 2 exp( ) 0808 . 0 6696 . 0 ( ) 10 0519 . 1 10 7206 . 1 exp( 5 . 0 10 2703 . 2 ) 0808 . 0 3282 . 1 ( ) 10 0519 . 1 10 1073 . 2 exp( ) 0808 . 0 6738 . 0 ( ) 10 0519 . 1 10 4634 . 3 exp( 5 . 0 10 4237 . 1 ) 0808 . 0 3332 . 1 ( ) 10 0519 . 1 10 1267 . 2 exp( ) 0808 . 0 6668 . 0 ( ) 10 0519 . 1 10 2857 . 2 exp( 5 . 0 10 4337 . 3 ) 0808 . 0 3469 . 1 ( ) 10 0519 . 1 10 1710 . 2 exp( ) 0808 . 0 6531 . 0 ( ) 10 0519 . 1 10 7170 . 6 exp( 5 . 0 34 24 14 23 13 12 32 6 4 6 6 410 32 6 4 6 6 310 32 6 4 6 6 210 33 6 4 6 6 110 = = = = = = × = . .. . . .. . × × + × × - = × = . .. . . .. . × × + × × - = × = . .. . . .. . × × + × × - = × = . .. . . .. . × × + × × - = - - - - - - - - - - - - - - - - - - - - S S S S S S erfc erfc F erfc erfc F erfc erfc F erfc erfc F ANL-NBS-HS-000026, Rev 00 III-7 April 2000 From Equation (III-1), the normalized solute concentrations at t=8,100 years are [ ] [ ] 3 32 32 33 32 32 2 32 33 2 3 4 2 3 2 2 4 * 4 * 4 4 2 32 32 32 33 3 2 2 4 * 3 * 3 3 2 32 33 3 4 * 2 * 2 2 4 33 4 * 1 * 1 1 10 0463 . 4 )] 10 2703 . 2 10 4237 . 1 ( 4412 . 0 ) 10 4337 . 3 10 2703 . 2 ( 1316 . 0 ) 10 4237 . 1 10 4818 . 1 ( 5 . 1 ) 10 4818 . 1 10 4337 . 3 ( 1666 . 0 [ )] 10 4908 . 4 10 9903 . 9 ( 4412 . 0 ) 10 8140 . 2 10 4908 . 4 ( 1316 . 0 ) 10 9902 . 9 10 5318 . 1 ( 5 . 1 ) 10 5318 . 1 10 8140 . 2 ( 1666 . 0 [ ) 3100 , 1 ( ) 8100 , 1 ( 10 8135 . 3 )] 10 4237 . 1 10 2703 . 2 ( 7649 . 1 ) 10 2703 . 2 10 4337 . 4 ( 5264 . 0 [ )] 10 9902 . 9 10 4908 . 4 ( 7649 . 1 ) 10 4908 . 4 10 8140 . 2 ( 5264 . 0 [ ) 3100 , 1 ( ) 8100 , 1 ( 10 4563 . 1 ) 10 4237 . 1 10 4337 . 3 ( 5 . 1 ) 10 9902 . 9 10 8140 . 2 ( 5 . 1 ) 3100 , 1 ( ) 8100 , 1 ( 10 8140 . 2 10 4337 . 3 10 8140 . 2 ) 3100 , 1 ( ) 8100 , 1 ( - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - × = × - × - × - × - × - × - × - × - - × - × - × - × - × - × - × - × - = - = × = × - × + × - × - × - × + × - × = - = × = × - × - - × - × - = - = × = × - × == - = c c c c c c c c c c c c ANL-NBS-HS-000026, Rev 00 III-8 April 2000 At time t=12300 years 0 10 9434 . 1 ) 1609 . 0 3160 . 2 ( ) 10 0519 . 1 10 1210 . 2 exp( ) 1609 . 0 3140 . 0 ( ) 10 0519 . 1 10 7206 . 1 exp( 5 . 0 10 1668 . 7 ) 1609 . 0 2991 . 2 ( ) 10 0519 . 1 10 1073 . 2 exp( ) 1609 . 0 2971 . 0 ( ) 10 0519 . 1 10 4634 . 3 exp( 5 . 0 10 1362 . 1 ) 1609 . 0 3220 . 2 ( ) 10 0519 . 1 10 1267 . 2 exp( ) 1609 . 0 3220 . 0 ( ) 10 0519 . 1 10 2857 . 2 exp( 5 . 0 10 6850 . 1 ) 1609 . 0 3765 . 2 ( ) 10 0519 . 1 10 1710 . 2 exp( ) 1609 . 0 3765 . 0 ( ) 10 0519 . 1 10 7170 . 6 exp( 5 . 0 34 24 14 23 13 12 1 6 4 6 6 410 1 6 4 6 6 310 1 6 4 6 6 210 3 6 4 6 6 110 = = = = = = × = . .. . . .. . × × + - × × - = × = . .. . . .. . × × + - × × - = × = . .. . . .. . × × + - × × - = × = . .. . . .. . × × + - × × - = - - - - - - - - - - - - - - - - - - - - S S S S S S erfc erfc F erfc erfc F erfc erfc F erfc erfc F For (t-t0)=12300-5000=7300 years 0 10 3187 . 1 ) 1240 . 0 7815 . 1 ( ) 10 0519 . 1 10 1210 . 2 exp( ) 1240 . 0 2206 . 0 ( ) 10 0519 . 1 10 7206 . 1 exp( 5 . 0 10 5172 . 3 ) 1240 . 0 7714 . 1 ( ) 10 0519 . 1 10 1073 . 2 exp( ) 1240 . 0 2306 . 0 ( ) 10 0519 . 1 10 4634 . 3 exp( 5 . 0 10 0281 . 9 ) 1239 . 0 7846 . 1 ( ) 10 0519 . 1 10 1267 . 2 exp( ) 1239 . 0 2154 . 0 ( ) 10 0519 . 1 10 2857 . 2 exp( 5 . 0 10 4569 . 3 ) 1239 . 0 8169 . 1 ( ) 10 0519 . 1 10 1710 . 2 exp( ) 1239 . 0 1831 . 0 ( ) 10 0519 . 1 10 7170 . 6 exp( 5 . 0 34 24 14 23 13 12 3 6 4 6 6 410 3 6 4 6 6 310 4 6 4 6 6 210 5 6 4 6 6 110 = = = = = = × = . .. . . .. . × × + × × - = × = . .. . . .. . × × + × × - = × = . .. . . .. . × × + × × - = × = . .. . . .. . × × + × × - = - - - - - - - - - - - - - - - - - - - - S S S S S S erfc erfc F erfc erfc F erfc erfc F erfc erfc F ANL-NBS-HS-000026, Rev 00 III-9 April 2000 From Equation (III-1), the normalized solute concentrations at t=12,300 years are [ ] [ ] 2 3 4 5 3 4 3 3 5 1 1 3 1 1 1 1 3 * 4 * 4 4 1 4 3 3 5 1 1 1 3 * 3 * 3 3 1 4 5 1 3 * 2 * 2 2 3 5 3 * 1 * 1 1 10 2658 . 8 )] 10 5172 . 3 10 0281 . 9 ( 4412 . 0 ) 10 4569 . 3 10 5172 . 3 ( 1316 . 0 ) 10 0281 . 9 10 3187 . 1 ( 5 . 1 ) 10 3187 . 1 10 4569 . 3 ( 1666 . 0 [ )] 10 1668 . 7 10 1362 . 1 ( 4412 . 0 ) 10 6850 . 1 10 1668 . 7 ( 1316 . 0 ) 10 1362 . 1 10 9434 . 1 ( 5 . 1 ) 10 9434 . 1 10 6850 . 1 ( 1666 . 0 [ ) 7300 , 1 ( ) 12300 , 1 ( 10 8338 . 6 )] 10 0281 . 9 10 5172 . 3 ( 7649 . 1 ) 10 5172 . 3 10 4569 . 3 ( 5264 . 0 [ )] 10 1362 . 1 10 1668 . 7 ( 7649 . 1 ) 10 1668 . 7 10 6850 . 1 ( 5264 . 0 [ ) 7300 , 1 ( ) 12300 , 1 ( 10 6660 . 1 ) 10 0281 . 9 10 4569 . 3 ( 5 . 1 ) 10 1362 . 1 10 6850 . 1 ( 5 . 1 ) 7300 , 1 ( ) 12300 , 1 ( 10 6504 . 1 10 4569 . 3 10 6850 . 1 ) 7300 , 1 ( ) 12300 , 1 ( - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - × = × - × - × - × - × - × - × - × - - × - × - × - × - × - × - × - × - = - = × = × - × + × - × - × - × + × - × = - = × = × - × - - × - × - = - = × = × - × == - = c c c c c c c c c c c c ANL-NBS-HS-000026, Rev 00 III-10 April 2000 At time t=17400 years 0 10 9482 . 1 ) 1914 . 0 8612 . 2 ( ) 10 0519 . 1 10 1210 . 2 exp( ) 1914 . 0 8592 . 0 ( ) 10 0519 . 1 10 7206 . 1 exp( 5 . 0 10 1946 . 7 ) 1914 . 0 8373 . 2 ( ) 10 0519 . 1 10 1073 . 2 exp( ) 1914 . 0 8353 . 0 ( ) 10 0519 . 1 10 4634 . 3 exp( 5 . 0 10 1385 . 1 ) 1913 . 0 8701 . 2 ( ) 10 0519 . 1 10 1267 . 2 exp( ) 1913 . 0 8701 . 0 ( ) 10 0519 . 1 10 2857 . 2 exp( 5 . 0 10 6856 . 1 ) 1913 . 0 9472 . 2 ( ) 10 0519 . 1 10 1710 . 2 exp( ) 1913 . 0 9472 . 0 ( ) 10 0519 . 1 10 7170 . 6 exp( 5 . 0 34 24 14 23 13 12 1 6 4 6 6 410 1 6 4 6 6 310 1 6 4 6 6 210 3 6 4 6 6 110 = = = = = = × = . .. . . .. . × × + - × × - = × = . .. . . .. . × × + - × × - = × = . .. . . .. . × × + - × × - = × = . .. . . .. . × × + - × × - = - - - - - - - - - - - - - - - - - - - - S S S S S S erfc erfc F erfc erfc F erfc erfc F erfc erfc F For (t-t0)=17400-5000=12400 years 0 10 9445 . 1 ) 1616 . 0 3267 . 2 ( ) 10 0519 . 1 10 1210 . 2 exp( ) 1616 . 0 3247 . 0 ( ) 10 0519 . 1 10 7206 . 1 exp( 5 . 0 10 1729 . 7 ) 1616 . 0 3097 . 2 ( ) 10 0519 . 1 10 1073 . 2 exp( ) 1616 . 0 3077 . 0 ( ) 10 0519 . 1 10 4634 . 3 exp( 5 . 0 10 1368 . 1 ) 1615 . 0 3327 . 2 ( ) 10 0519 . 1 10 1267 . 2 exp( ) 1615 . 0 3327 . 0 ( ) 10 0519 . 1 10 2857 . 2 exp( 5 . 0 10 6851 . 1 ) 1615 . 0 3877 . 2 ( ) 10 0519 . 1 10 1710 . 2 exp( ) 1615 . 0 3876 . 0 ( ) 10 0519 . 1 10 7170 . 6 exp( 5 . 0 34 24 14 23 13 12 1 6 4 6 6 410 1 6 4 6 6 310 1 6 4 6 6 210 3 6 4 6 6 110 = = = = = = × = . .. . . .. . × × + - × × - = × = . .. . . .. . × × + - × × - = × = . .. . . .. . × × + - × × - = × = . .. . . .. . × × + - × × - = - - - - - - - - - - - - - - - - - - - - S S S S S S erfc erfc F erfc erfc F erfc erfc F erfc erfc F ANL-NBS-HS-000026, Rev 00 III-11 April 2000 From Equation (III-1), the normalized solute concentrations at t=17,400 years are [ ] [ ] 4 1 1 3 1 1 1 1 3 1 1 3 1 1 1 1 3 * 4 * 4 4 3 1 1 1 3 1 1 1 3 * 3 * 3 3 4 1 3 1 3 * 2 * 2 2 7 3 3 * 1 * 1 1 10 5845 . 3 )] 10 1729 . 7 10 1368 . 1 ( 4412 . 0 ) 10 6851 . 1 10 1729 . 7 ( 1316 . 0 ) 10 1368 . 1 10 9445 . 1 ( 5 . 1 ) 10 9445 . 1 10 6851 . 1 ( 1666 . 0 [ )] 10 1946 . 7 10 1385 . 1 ( 4412 . 0 ) 10 6856 . 1 10 1946 . 7 ( 1316 . 0 ) 10 1385 . 1 10 9482 . 1 ( 5 . 1 ) 10 9482 . 1 10 6856 . 1 ( 1666 . 0 [ ) 12400 , 1 ( ) 17400 , 1 ( 10 3878 . 2 )] 10 1368 . 1 10 1729 . 7 ( 7649 . 1 ) 10 1729 . 7 10 6851 . 1 ( 5264 . 0 [ )] 10 1385 . 1 10 1946 . 7 ( 7649 . 1 ) 10 1946 . 7 10 6856 . 1 ( 5264 . 0 [ ) 12400 , 1 ( ) 17400 , 1 ( 10 5425 . 2 ) 10 1368 . 1 10 6851 . 1 ( 5 . 1 ) 10 1385 . 1 10 6856 . 1 ( 5 . 1 ) 12400 , 1 ( ) 17400 , 1 ( 10 0 . 5 10 6851 . 1 10 6856 . 1 ) 12400 , 1 ( ) 17400 , 1 ( - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - × = × - × - × - × - × - × - × - × - - × - × - × - × - × - × - × - × - = - = × = × - × + × - × - × - × + × - × = - = × = × - × - - × - × - = - = × = × - × == - = c c c c c c c c c c c c ANL-NBS-HS-000026, Rev 00 III-12 April 2000 Table III-1. Comparison of CHAIN and Hand Calculation Results Species 1 Species 2 Species 3 Species 4 Time (year) Hand Calculation CHAIN Hand Calculation CHAIN Hand Calculation CHAIN Hand Calculation CHAIN 3000 3.6631E-35 3.663E-35 1.6273E-34 1.627E-34 4.2714E-35 4.2703E-34 2.3111E-35 2.312E-34 81000 2.8140E-4 2.814E-4 1.4563E-2 1.456E-2 3.8135E-2 3.814E-2 4.0463E-3 4.026E-3 12300 1.6504E-3 1.650E-3 1.6660E-1 1.666E-1 6.8338E-1 6.853E-1 8.2658E-2 8.240E-2 174000 5.000E-7 4.791E-7 2.5425E-4 2.571E-4 2.3878E-3 2.390E-3 3.5845E-4 3.608E-4 The hand calculation results compared well against the results from the software routine CHAIN. The small difference was due to truncation errors in hand calculations, as software routine CHAIN used double precision in the code. The comparison proved that CHAIN was correct. The analytical solution formulated in CHAIN had certain limitations on the range of parameters. Those limitations are listed below. (1) The number of species in the problem must be less or equal to 4. (2) The difference of (Rij .k- µij) does not equal 0. Where Rij=Ri-Rj and µij= µiRi- µjRj with i, j=1, 2, 3, and 4, but i does not equal j, and k=1, 2, and 3. This constraint prevents some of the parameters from being divided by 0. ANL-NBS-HS-000026, Rev 00 IV-1 April 2000 ATTACHMENT IV. SOFTWARE ROUTINE TEST_EBS_RANDOM The test routine, test_ebs_random.f Version 1.0, was written in Fortran 90. The core of the test routine is the EBS random release model with added code for generating random node locations within a given field, loading nodes into an array containing node coordinates, and randomly assigning each node an index. This attachment documents that the software routine performs as expected. LIST OF INPUTS AND SCREEN OUTPUTS OF THE TEST PROBLEM Listed below are the screen input and output of the test problem. Input total Number of Nodes 40 Input number of sub-regions 2 Input the seed for the random number geneator 13131 Input total number of M_fine failed packages 6 Input number of species 1 Released at this test (->1) or (->0) 1 For sub-region: 1 Input xmin and xmax 0 100 Input ymin and ymax 0 100 Input number of nodes in zone 1 20 How many M_fine failed package in this region 3 For sub-region: 2 Input xmin and xmax 120 250 Input ymin and ymax 120 250 Input number of nodes in zone 2 20 How many M_fine failed package in this regio 3 Node coordinates generated Input how many time steps to simulate 4 Time step: 1 Input data for zone 1 Input number of N_failed package 10 Input data for zone 2 Input number of N_failed package 10 Time step: 2 Input data for zone 1 Input number of N_failed package 15 Input data for zone 2 Input number of N_failed package 15 Time step: 3 Input data for zone 1 Input number of N_failed package 17 ANL-NBS-HS-000026, Rev 00 IV-2 April 2000 Input data for zone 2 Input number of N_failed package 17 Time step: 4 Input data for zone 1 Input number of N_failed package 20 Input data for zone 2 Input number of N_failed package 20 ANL-NBS-HS-000026, Rev 00 IV-3 April 2000 LIST OF OUTPUT RESULTS CONTAINED IN FILE NODE 40 !total_nodes, total number of nodes 2 !N_large, number of sub-regions 13131 !seed, for random number generator 6 !M_fine, number of M_fine failed packages 1 !nspecies, number of species 1 !irelease, Released at this test (->1 or ->0) !Input parameters for sub-region: 1 ( 0.00, 100.00) !xmin, xmax, X-range ( 0.00, 100.00) !ymin, ymax, Y-range 20 !number of nodes in zone 1 8 92.76 97.69 !node_index,(x,y) 3 58.12 70.45 !node_index,(x,y) 11 14.59 94.89 !node_index,(x,y) 13 78.99 58.85 !node_index,(x,y) 19 11.36 27.01 !node_index,(x,y) 1 30.52 80.25 !node_index,(x,y) 7 56.41 46.25 !node_index,(x,y) 24 35.86 40.80 !node_index,(x,y) 23 37.75 59.77 !node_index,(x,y) 20 40.77 74.52 !node_index,(x,y) 10 7.88 97.37 !node_index,(x,y) 5 98.07 5.48 !node_index,(x,y) 37 36.23 3.17 !node_index,(x,y) 6 44.49 44.11 !node_index,(x,y) 35 85.75 43.51 !node_index,(x,y) 16 3.88 54.84 !node_index,(x,y) 14 66.48 84.95 !node_index,(x,y) 18 11.24 99.72 !node_index,(x,y) 25 19.51 10.98 !node_index,(x,y) 15 22.46 11.69 !node_index,(x,y) 3 !number of individually failed package 93.58 48.74 !(x,y) of failed package 34.24 3.88 !(x,y) of failed package 24.39 75.84 !(x,y) of failed package !Input parameters for sub-region: 2 ( 120.00, 250.00) !xmin, xmax, X-range ( 120.00, 250.00) !ymin, ymax, Y-range 20 !number of nodes in zone 2 32 235.34 236.86 !node_index,(x,y) 40 128.89 174.41 !node_index,(x,y) 4 209.61 231.01 !node_index,(x,y) 29 213.41 182.47 !node_index,(x,y) 31 186.02 120.08 !node_index,(x,y) 9 190.41 188.59 !node_index,(x,y) 30 209.18 241.93 !node_index,(x,y) 39 223.25 163.77 !node_index,(x,y) 22 174.71 206.49 !node_index,(x,y) 27 136.03 219.27 !node_index,(x,y) 12 172.61 156.08 !node_index,(x,y) 26 136.88 147.70 !node_index,(x,y) 38 133.46 193.26 !node_index,(x,y) 33 181.27 144.80 !node_index,(x,y) 21 151.83 218.28 !node_index,(x,y) 17 186.94 184.39 !node_index,(x,y) 2 134.16 124.11 !node_index,(x,y) 34 196.44 167.73 !node_index,(x,y) 36 217.48 226.39 !node_index,(x,y) ANL-NBS-HS-000026, Rev 00 IV-4 April 2000 28 129.11 158.84 !node_index,(x,y) 3 !number of individually failed package 205.62 217.01 !(x,y) of failed package 147.41 152.02 !(x,y) of failed package 135.88 207.18 !(x,y) of failed package Node coordinates generated 4 !nsteps, number of time steps At time step: 1 !istp, time step 10 !Number of failed packages in sub-region: 1 10 !Number of failed packages in sub-region: 2 !locating node having least dist. to the given (x,y) 93.58 48.74 35 1 !(x,y),node index,sub-region 34.24 3.88 37 1 !(x,y),node index,sub-region 24.39 75.84 1 1 !(x,y),node index,sub-region 205.62 217.01 4 2 !(x,y),node index,sub-region 147.41 152.02 26 2 !(x,y),node index,sub-region 135.88 207.18 27 2 !(x,y),node index,sub-region 3 !Fine packages failed in N_large zone: 1 35 -1.00 !node index,pcnsk 37 -1.00 !node index,pcnsk 1 -1.00 !node index,pcnsk 10 !Packages failed in N_large zone: 1 25 -1.00 !node index,pcnsk 24 -1.00 !node index,pcnsk 20 -1.00 !node index,pcnsk 23 -1.00 !node index,pcnsk 15 -1.00 !node index,pcnsk 10 -1.00 !node index,pcnsk 13 -1.00 !node index,pcnsk 19 -1.00 !node index,pcnsk 18 -1.00 !node index,pcnsk 7 -1.00 !node index,pcnsk 3 !Fine packages failed in N_large zone: 2 4 -1.00 !node index,pcnsk 26 -1.00 !node index,pcnsk 27 -1.00 !node index,pcnsk 10 !Packages failed in N_large zone: 2 31 -1.00 !node index,pcnsk 30 -1.00 !node index,pcnsk 29 -1.00 !node index,pcnsk 9 -1.00 !node index,pcnsk 34 -1.00 !node index,pcnsk 39 -1.00 !node index,pcnsk 12 -1.00 !node index,pcnsk 32 -1.00 !node index,pcnsk 28 -1.00 !node index,pcnsk 17 -1.00 !node index,pcnsk At time step: 2 !istp, time step 15 !Number of failed packages in sub-region: 1 15 !Number of failed packages in sub-region: 2 ANL-NBS-HS-000026, Rev 00 IV-5 April 2000 !locating node having least dist. to the given (x,y) 3 !Fine packages failed in N_large zone: 1 35 -1.00 !node index,pcnsk 37 -1.00 !node index,pcnsk 1 -1.00 !node index,pcnsk 15 !Packages failed in N_large zone: 1 25 -1.00 !node index,pcnsk 24 -1.00 !node index,pcnsk 20 -1.00 !node index,pcnsk 23 -1.00 !node index,pcnsk 15 -1.00 !node index,pcnsk 10 -1.00 !node index,pcnsk 13 -1.00 !node index,pcnsk 19 -1.00 !node index,pcnsk 18 -1.00 !node index,pcnsk 7 -1.00 !node index,pcnsk 6 -1.00 !node index,pcnsk 16 -1.00 !node index,pcnsk 5 -1.00 !node index,pcnsk 14 -1.00 !node index,pcnsk 8 -1.00 !node index,pcnsk 3 !Fine packages failed in N_large zone: 2 4 -1.00 !node index,pcnsk 26 -1.00 !node index,pcnsk 27 -1.00 !node index,pcnsk 15 !Packages failed in N_large zone: 2 31 -1.00 !node index,pcnsk 30 -1.00 !node index,pcnsk 29 -1.00 !node index,pcnsk 9 -1.00 !node index,pcnsk 34 -1.00 !node index,pcnsk 39 -1.00 !node index,pcnsk 12 -1.00 !node index,pcnsk 32 -1.00 !node index,pcnsk 28 -1.00 !node index,pcnsk 17 -1.00 !node index,pcnsk 36 -1.00 !node index,pcnsk 40 -1.00 !node index,pcnsk 22 -1.00 !node index,pcnsk 38 -1.00 !node index,pcnsk 33 -1.00 !node index,pcnsk At time step: 3 !istp, time step 17 !Number of failed packages in sub-region: 1 17 !Number of failed packages in sub-region: 2 !locating node having least dist. to the given (x,y) 3 !Fine packages failed in N_large zone: 1 35 -1.00 !node index,pcnsk 37 -1.00 !node index,pcnsk 1 -1.00 !node index,pcnsk 17 !Packages failed in N_large zone: 1 25 -1.00 !node index,pcnsk 24 -1.00 !node index,pcnsk 20 -1.00 !node index,pcnsk ANL-NBS-HS-000026, Rev 00 IV-6 April 2000 23 -1.00 !node index,pcnsk 15 -1.00 !node index,pcnsk 10 -1.00 !node index,pcnsk 13 -1.00 !node index,pcnsk 19 -1.00 !node index,pcnsk 18 -1.00 !node index,pcnsk 7 -1.00 !node index,pcnsk 6 -1.00 !node index,pcnsk 16 -1.00 !node index,pcnsk 5 -1.00 !node index,pcnsk 14 -1.00 !node index,pcnsk 8 -1.00 !node index,pcnsk 3 -1.00 !node index,pcnsk 11 -1.00 !node index,pcnsk 3 !Fine packages failed in N_large zone: 2 4 -1.00 !node index,pcnsk 26 -1.00 !node index,pcnsk 27 -1.00 !node index,pcnsk 17 !Packages failed in N_large zone: 2 31 -1.00 !node index,pcnsk 30 -1.00 !node index,pcnsk 29 -1.00 !node index,pcnsk 9 -1.00 !node index,pcnsk 34 -1.00 !node index,pcnsk 39 -1.00 !node index,pcnsk 12 -1.00 !node index,pcnsk 32 -1.00 !node index,pcnsk 28 -1.00 !node index,pcnsk 17 -1.00 !node index,pcnsk 36 -1.00 !node index,pcnsk 40 -1.00 !node index,pcnsk 22 -1.00 !node index,pcnsk 38 -1.00 !node index,pcnsk 33 -1.00 !node index,pcnsk 21 -1.00 !node index,pcnsk 2 -1.00 !node index,pcnsk At time step: 4 !istp, time step 20 !Number of failed packages in sub-region: 1 20 !Number of failed packages in sub-region: 2 !locating node having least dist. to the given (x,y) 3 !Fine packages failed in N_large zone: 1 35 -1.00 !node index,pcnsk 37 -1.00 !node index,pcnsk 1 -1.00 !node index,pcnsk 17 !Packages failed in N_large zone: 1 25 -1.00 !node index,pcnsk 24 -1.00 !node index,pcnsk 20 -1.00 !node index,pcnsk 23 -1.00 !node index,pcnsk 15 -1.00 !node index,pcnsk 10 -1.00 !node index,pcnsk 13 -1.00 !node index,pcnsk 19 -1.00 !node index,pcnsk 18 -1.00 !node index,pcnsk 7 -1.00 !node index,pcnsk 6 -1.00 !node index,pcnsk 16 -1.00 !node index,pcnsk 5 -1.00 !node index,pcnsk ANL-NBS-HS-000026, Rev 00 IV-7 April 2000 14 -1.00 !node index,pcnsk 8 -1.00 !node index,pcnsk 3 -1.00 !node index,pcnsk 11 -1.00 !node index,pcnsk 3 !Fine packages failed in N_large zone: 2 4 -1.00 !node index,pcnsk 26 -1.00 !node index,pcnsk 27 -1.00 !node index,pcnsk 17 !Packages failed in N_large zone: 2 31 -1.00 !node index,pcnsk 30 -1.00 !node index,pcnsk 29 -1.00 !node index,pcnsk 9 -1.00 !node index,pcnsk 34 -1.00 !node index,pcnsk 39 -1.00 !node index,pcnsk 12 -1.00 !node index,pcnsk 32 -1.00 !node index,pcnsk 28 -1.00 !node index,pcnsk 17 -1.00 !node index,pcnsk 36 -1.00 !node index,pcnsk 40 -1.00 !node index,pcnsk 22 -1.00 !node index,pcnsk 38 -1.00 !node index,pcnsk 33 -1.00 !node index,pcnsk 21 -1.00 !node index,pcnsk 2 -1.00 !node index,pcnsk ANL-NBS-HS-000026, Rev 00 IV-8 April 2000 LIST OF TEST ROUTINE test_ebs_random.f program test_ebs_random !routine used to test the EBS random release model implicit none integer:: numinbuffs, numoutbuffs,idaughter,ithp,ip,num_left integer:: j,index_in,index_out,memleft,lns,lns1,tmp_ptindex integer:: M_fine, N_large, nodes_left,selected_node,irelease integer:: newly_failed,trans_node, in_buffer,M_add_i,M_add_N integer:: i,k,n,nodes, nzone, seed, itmp,total_nodes,pseudo_node integer:: ith_M, nodes_in_region,nspecies,current_node,temp_node integer:: M_failed,istp,nsteps,ioverlap,idrop,N_failed integer, allocatable, dimension(:)::ptindex,insnode,M_N_region integer, allocatable::M_N_failed_nodes(:),node_list(:),i_fine(:) real:: yearsleft,scaleconpart,sumdecayout,sumdecayin real:: decayin,decayout,rescale,conpart,conpartd4,adjconpart real:: x_fine, y_fine, x_dis, y_dis, xy_dis, xy_least_dis real:: xmin,ymin,xmax,ymax,xlen,ylen,begin_time,end_time real, allocatable, dimension(:,:):: cord real, allocatable, dimension(:):: pcnsk,in,t2sk,t1sk real*8:: pinmass,sumpinmass,tmp_pcnsk !set up parameters. cli_added 8/11/1999 real, parameter:: pthresh=0.5, pthresh2=0.25,add0_5=0.5 real, parameter:: yearindays=365.25, dayinsecs=86400. real, parameter:: large_dis=1.E9 !input parameters for generating coordinates open(7,file='nodes') print *, 'Input total Number of Nodes' read(5,*)total_nodes write(7,111)total_nodes 111 format(I7,33x,'!total_nodes, total number of nodes') print *,'Input number of sub-regions' read(5,*)N_large write(7,112)N_large 112 format(I7,33x,'!N_large, number of sub-regions') print *, 'Input the seed for the random number geneator' read(5,*)seed write(7,113)seed 113 format(I7,33x,'!seed, for random number generator') print *,'Input total number of M_fine failed packages' ANL-NBS-HS-000026, Rev 00 IV-9 April 2000 read(5,*)M_fine write(7,114)M_fine 114 format(I7,33x,'!M_fine, number of M_fine failed packages') print *,'Input number of species' read(5,*)nspecies write(7,116)nspecies 116 format(I7,33x,'!nspecies, number of species') print *,'Released at this test (->1) or (->0)' read(5,*)irelease write(7,117)irelease 117 format(I7,33x,'!irelease, Released at this test (->1 or ->0)'/) !allocate array and initialize the coordinate arrays total_nodes=total_nodes+1 if(.not.allocated(cord))then allocate(cord(total_nodes,2),ptindex(total_nodes), & pcnsk(total_nodes),insnode(total_nodes)) itmp=9+M_fine*2+N_large+M_fine*N_large allocate(in(itmp),t1sk(total_nodes),t2sk(total_nodes)) allocate(node_list(total_nodes)) endif total_nodes=total_nodes-1 do i=1,total_nodes node_list(i)=i enddo in(4)=M_fine in(5+M_fine*2)=N_large in(7+M_fine*2+N_large)=irelease cord(1:total_nodes,1)=-1. cord(1:total_nodes,2)=-1. !generate random node coordinates insnode(1)=0 itmp=0 current_node=0 do i=1,N_large print *, 'For sub-region: ',i print *, 'Input xmin and xmax' read(5,*)xmin,xmax print *,'Input ymin and ymax' read(5,*)ymin,ymax print *, 'Input number of nodes in zone',i read(5,*)nodes write(7,118)i write(7,119)xmin,xmax write(7,120)ymin,ymax write(7,121)nodes,i 118 format(/'!Input parameters for sub-region: ',I7/) 119 format('(',F10.2,',',F10.2,')',17x,'!xmin, xmax, X-range') ANL-NBS-HS-000026, Rev 00 IV-10 April 2000 120 format('(',F10.2,',',F10.2,')',17x,'!ymin, ymax, Y-range'/) 121 format(I7,33x,'!number of nodes in zone',I7/) xlen=xmax-xmin ylen=ymax-ymin do n=1,nodes num_left=total_nodes-current_node pseudo_node=num_left*ran(seed)+0.5 pseudo_node=current_node+pseudo_node if(pseudo_node == current_node)pseudo_node=current_node+1 current_node=current_node+1 temp_node=node_list(current_node) node_list(current_node)=node_list(pseudo_node) node_list(pseudo_node)=temp_node pseudo_node=node_list(current_node) itmp=itmp+1 cord(pseudo_node,1)=xmin+xlen*ran(seed) cord(pseudo_node,2)=ymin+ylen*ran(seed) write(7,122)pseudo_node,cord(pseudo_node,1), & cord(pseudo_node,2) 122 format(I7,3x,2(F10.2,3x), 4x,'!node_index,(x,y)') ptindex(itmp)=pseudo_node pcnsk(itmp)=-1 enddo insnode(i+1)=itmp print *,'How many M_fine failed package in this region' read(5,*)M_failed write(7,123)M_failed 123 format(/,I7,33x,'!number of individually failed package') do n=1,M_failed ith_M=ith_m+1 if(ith_M <= M_fine)then in(5+(ith_M-1)*2)=xmin+xlen*ran(seed) in(4+ith_m*2)=ymin+ylen*ran(seed) write(7,124)in(5+(ith_M-1)*2),in(4+ith_M*2) 124 format(2(F10.2,3x),14x'!(x,y) of failed package') else print *,'Maximum number of M_fine already reached' write(7,*)'Maximum number of M_fine already reached' exit endif enddo enddo print *,'Node coordinates generated' write(7,*)'' write(7,125) 125 format('Node coordinates generated') print *,'Input how many time steps to simulate' read(5,*)nsteps write(7,126)nsteps 126 format(/,I7,33x,'!nsteps, number of time steps'/) do istp=1,nsteps ANL-NBS-HS-000026, Rev 00 IV-11 April 2000 print *,'Time step: ', istp write(7,127)istp 127 format(/,'At time step:',I5,22x,'!istp, time step'/) ith_m=0 do j=1,N_large print *,'Input data for zone ',j print *,'Input number of N_failed package' read(5,*)N_failed write(7,128)N_failed,j 128 format(I7,33x,'!Number of failed packages in sub-region:',I7/) in(5+M_fine*2+j)=N_failed enddo write(7,129) 129 format('!locating node having least dist. to the given (x,y)'/) !get injection time in days t1sk(1) = begin_time/dayinsecs t2sk(1) = end_time/dayinsecs !get # of buffer infromation from RIP, cli_08/27/1999 if(in(7+M_fine*2+N_large) == 0)stop 'No mass will be released' !for random ebs release, get M_fine_group and N_large_group # M_fine=in(4) N_large=in(5+M_fine*2) !allocate arrays for EBS random release. cli_added 8/31/1999 if(.not. allocated(M_N_region))then M_add_N=M_fine+N_large allocate (M_N_region(M_add_N), M_N_failed_nodes(M_add_N)) allocate (i_fine(N_large)) i_fine=0 M_N_region=0 M_N_failed_nodes=0 endif !Select release nodes for random EBS release !Get M_fine group (x,y) and locate the node index and region if(istp == 1 )then in_buffer=4 do i=1,M_fine ioverlap=0 in_buffer=in_buffer+1 x_fine=in(in_buffer) in_buffer=in_buffer+1 y_fine=in(in_buffer) xy_least_dis=large_dis do j=1,N_large do k=insnode(j)+1,insnode(j+1) x_dis=cord(ptindex(k),1)-x_fine ANL-NBS-HS-000026, Rev 00 IV-12 April 2000 y_dis=cord(ptindex(k),2)-y_fine xy_dis=sqrt(x_dis*x_dis+y_dis*y_dis) if(xy_dis < xy_least_dis)then M_N_failed_nodes(i)=-k M_N_region(i)=j xy_least_dis=xy_dis endif enddo enddo !check for overlap do j=1,i-1 if(M_N_failed_nodes(i) == M_N_failed_nodes(j))then ioverlap=ioverlap+1 endif enddo if(ioverlap /= 0)then print *, 'Found',ioverlap,' overlaped locations' write(7,130) 130 format(/'!Found overlaped node'/) write(7,132)x_fine,y_fine,ptindex(-M_N_failed_nodes(i)), & M_N_region(i) write(7,131) 131 format(/'!Program terminated, reassign values'/) stop 'Program terminated, reassign values' endif !shuffle the nodes for M_fine group selected_node=insnode(M_N_region(i)+1)-i_fine(M_N_region(i)) trans_node=-M_N_failed_nodes(i) tmp_ptindex=ptindex(selected_node) ptindex(selected_node)=ptindex(trans_node) ptindex(trans_node)=tmp_ptindex tmp_pcnsk=pcnsk(selected_node) pcnsk(selected_node)=pcnsk(trans_node) pcnsk(trans_node)=tmp_pcnsk M_N_failed_nodes(i)=-selected_node i_fine(M_N_region(i))=i_fine(M_N_region(i))+1 write(7,132)x_fine,y_fine,ptindex(-M_N_failed_nodes(i)), & M_N_region(i) 132 format(2(F10.2),2(I7,3x),'!(x,y),node index,sub-region') enddo endif !For N_large groups: randomly select N_failed package nodes in_buffer=5+M_fine*2 do i=1, N_large in_buffer=in_buffer+1 N_failed=in(in_buffer) M_add_i=M_fine+i nodes_in_region=insnode(i+1)-insnode(i)-i_fine(i) ANL-NBS-HS-000026, Rev 00 IV-13 April 2000 newly_failed=N_failed-M_N_failed_nodes(M_add_i) nodes_left=nodes_in_region-M_N_failed_nodes(M_add_i) if(N_failed < nodes_in_region)then M_N_region(M_add_i)=i do j=1,newly_failed if(newly_failed > 1 )then do k=int(ran(seed)*nodes_left+add0_5) if(k /= 0)exit enddo else k=1 endif selected_node=insnode(i)+M_N_failed_nodes(M_add_i)+k trans_node=insnode(i)+M_N_failed_nodes(M_add_i)+1 tmp_ptindex=ptindex(trans_node) ptindex(trans_node)=ptindex(selected_node) ptindex(selected_node)=tmp_ptindex tmp_pcnsk=pcnsk(trans_node) pcnsk(trans_node)=pcnsk(selected_node) pcnsk(selected_node)=tmp_pcnsk nodes_left=nodes_left-1 M_N_failed_nodes(M_add_i)=M_N_failed_nodes(M_add_i)+1 enddo else M_N_failed_nodes(M_add_i)=nodes_in_region M_N_region(M_add_i)=i endif !output results write(7,135)i_fine(i),i 135 format(/,I7,3x,'!Fine packages failed in N_large zone:',I7/) itmp=insnode(i+1) do n=1,i_fine(i) write(7,136)ptindex(itmp),pcnsk(itmp) 136 format(I7,5x,F6.2,22x,'!node index,pcnsk') itmp=itmp-1 enddo write(7,137)M_N_failed_nodes(M_add_i),i 137 format(/,I7,3x,'!Packages failed in N_large zone:',5x,I7/) do n=insnode(i)+1,(insnode(i)+M_N_failed_nodes(M_add_i)) write(7,136)ptindex(n),pcnsk(n) enddo enddo enddo end program test_ebs_random