Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media Rev 00, ICN 00 ANL-NBS-HS-000022 July 2000 1 PURPOSE The purpose of this dispersion analysis is to provide estimates of the transverse and longitudinal dispersion that may occur at the sub-gridblock scale within the saturated zone (SZ) site-scale model. This analysis explicitly applies only to that part of the SZ flow path that is located within the fractured Tertiary volcanic rocks in the SZ down-gradient from beneath the potential repository. The amount of dispersion is characterized by values of dispersivity for each of the three principal directions. Estimates of dispersivity derived from this analysis are compared to estimates of dispersivities derived from the analytical solutions; these estimates can be used directly, or with modification, in performance assessment calculations. This analysis is designed to use existing modeling and analysis results as the basis for heterogeneous models of fracture permeability at the sub gridblock scale. The scope of this modeling exercise includes site scale model gridblocks containing both “background” fracture permeability and “enhanced” permeability: those dominated by large, throughgoing high permeability features. Two different source geometries are used: point and distributed planar. This analysis is governed by the CRWMS M&O "Analysis/Model Report (AMR) Work Direction and Planning Document" for, Geostatistical and Heterogeneity Analysis for the SZ Flow and Transport Model” (McKenna, 1999). The work presented in this AMR comprises a component of the SZ Process Model Report (PMR). 2 QUALITY ASSURANCE The activities documented in this Analysis/Model Report (AMR) were evaluated in accordance with QAP 2-0, Conduct of Activities, and were determined to be subject to the requirements of the U.S. DOE Office of Civilian Radioactive Waste Management (OCRWM) Quality Assurance Requirements and Description DOE/RW-0333P (QARD) (DOE 1998). This evaluation is documented in CRWMS M&O 1999a, b, Wemheuer, 1999 (activity evaluation for Work Package #WP1401213SM1). This AMR has been prepared in accordance with procedures identified in CRWMSM&O "Analysis and Modeling Report (AMR) Work Direction and Planning Document for Geostatistical and Heterogeneity Analysis for the SZ Flow and Transport Model ”(McKenna, 1999) and with procedure AP-3.10Q Analyses and Models. 3 COMPUTER SOFTWARE AND MODEL USAGE It has been determined that AP-SI.1Q applies to the software routines used in this work. 3.1 INDUSTRY STANDARD SOFTWARE The following industry standard software was used in this analysis and documentation: EXCEL 97-SR-1 Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 10 of 60 Used for spreadsheet calculations and calculation of basic statistics on the output of the calc_D_old and calc_cdf routines described below. SigmaPlot 4.0 Used for plotting and visualization of analysis results from calc_D_old (described below) in figures shown in this report. 3.2 PERMEABILITY MODELING The background permeability models are created with the indicator geostatistical simulation routine sisim. The enhanced permeability models are created by combining the background permeability models with the results of the ellip2 routine. Two different sets of UNIX shells are used to run the codes that build the permeability models and then accomplish post processing. One set is used for the background permeability models and one set for the enhanced permeability models. These shells are discussed below. 3.2.1 Geostatistical Simulations The routine sisim (version 2.0) was used to create the background fracture permeability models. This routine is included as part of the Geostatistical Software Library (GSLIB) (Deutsch and Journel, 1998). This subroutine is discussed in the Section 6.1.2.1, and the source code for this subroutine is listed in Attachment 4. A Pentium personal computer was used to run the sisim code. The routine modgeom (version 1.0) is used to change the random number seed after each call to the sisim or ellip2 routines. The seed is decremented by a constant value as set in the code. For the work reported herein, the seed is decremented by four each time modgeom is called. 3.2.2 Object Simulations The routine ellip2_MC (version 1.0) was used to create the disk shaped objects that are combined with the background fracture permeability models to create the enhanced permeability models. ellip2_MC is an extension of the GSLIB routine ellipsim. The feature added for this work is the ability to draw object orientations from a Gaussian distribution. The original ellipsim routine can only produce features with a constant orientation. The ellip2_MC subroutine is discussed in Section 6.1.2.2, and the source code for this subroutine is listed in Attachment 5. ellip2_MC was run on a Sun UltraSparc UNIX workstation. The enhanced fracture permeability models are created by combining the background permeability models with the disk-shaped objects created in ellip2_MC. Combining the two sets of simulation results is accomplished with the routine combine (version 1.0). This routine sets the log10 permeability to a constant value of –11.0 within each of the disk-shaped objects and writes out the log10 permeability field as an ascii file in GeoEAS format. This subroutine is discussed in Section 6.1.2.2, and the source code for this subroutine is listed in Attachment 6. combine was run on a SUN UltraSparc workstation. Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 11 of 60 3.2.3 Validation Checks Several postprocessing checks are made in an attempt to validate the fracture permeability models. These checks include the cumulative distribution function (cdf) and the mean spacing between high permeability features. These checks are made using the calc_cdf (version1.0) and gamfil (version 1.0) routines respectively. The use of these subroutines is discussed in Section 6.1.3.2, and the source code for both subroutines is listed in Attachments 7 and 9 respectively. 3.2.4 Reformatting The output files from either sisim or combine are reformatted to the input format required by FEHM by using the routine gs2fehm (version 1.0). This subroutine is discussed in section 6.1.3.2, and the source code for this subroutines is listed in Attachment 8. 3.3 FLOW MODELING The code FEHM Version 2.00 (Software tracking number 10031-2.00-00) (Zyvoloski et al. 1997a; 1997b). was obtained from the Software Congifuration Management and is used to solve for a steady state flow solution on the heterogeneous fields created with the sgsim and ellip2_MC codes. Solving for a steady-state solution is an appropriate use of FEHM and is within the range of FEHM. Different Unix shells, run_50_back.ksh and run_50_feat.ksh are used to run the flow model codes and do the particle tracking on the background permeability realizations and the enhanced permeability model realizations respectively. These shells are discussed in Section 6.2.2. The source codes of run_50_back.ksh and run_50_feat.ksh are listed in Attachment 11 and 12 respectively. The intrinsic macro flxo is activated within FEHM for these flow solutions. The flxo macro writes the internodal flux between each pair of selected nodes within the model. The maximum number of node pairs that can be written to the output file ( het_flow.out) is 500 The UNIX directories containing the files necessary for the flow and particle tracking as well as the calculation of the dispersivities are written to the same CD-ROM that contains the attachments as electronic files. Memory limitations do not allow for all flow and particle tracking output files to be retained. Therefore, these directories contain example input and output files for the flow and particle tracking solutions created with FEHM. The directory path for these files on the CD-ROM is: /AMR/FEHM/ for the flow solutions. FEHM is run on a Sun UltraSparc workstation to produce the flow solutions. 3.4 PARTICLE TRACKING The code FEHM Version 2.00 (Software tracking number 10031-2.00-00)is used to determine the streamlines (particle tracks) using the steady-state flow solution created in the previous step. There are two different particle-tracking routines within FEHM. This study uses the FEHM sptr macro for particle tracking. The sptr macro is an intrinsic macro within FEHM. Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 12 of 60 The flux weighted particle release locations for the planar source are determined with the routine cr8sptr (version 1.0).. This routine reads in the output of the flxo macro in FEHM (het_flow.out) and determines the flux weighted starting locations for the particles along the plane. The output of the cr8sptr routine is the file sptr.dat. FEHM uses the sptr.dat file as an input file to define the particle starting locations for the intrinsic sptr macro. The cr8sptr subroutine is discussed in Section 6.2.2, and the source code for this subroutine is listed as Attachment 13. The flux weighted particle release locations for the point source are determined with the routine cr8sptr_cell (version 1.0). This routine is a variation of cr8sptr used for the planar source. Input and output from this routine are the het_flow.out and sptr.dat files respectively. This subroutine is discussed in the section 6.2.2, and the source code for this subroutine is listed as Attachment 14. The directory path for the particle tracking runs on the CD-ROM is: /AMR/FEHM/track/. Example input and output files for the particle tracking are contained there. The particle tracking is run on a Sun UltraSparc UNIX workstation. 3.5 DISPERSION CALCULATIONS The dispersion calculations are done with the routine calc_D_old (version 1.0). The source code for this subroutine is listed as Attachment 15. The routine takes the het_flow.sptr2 output file from FEHM as input and calculates the longitudinal and transverse dispersivities at each collection plane. The “old” in the title of this routine refers to the current form of the output from FEHM contained in het_flow.sptr2. A new output format (“new”) in a developmental version of FEHM was also examined for these calculations, but due to excessive output, was not used for the final runs. Output from calc_D_old is written to four different files: disper_L.out, disper_X.out, disper_Z.out and debug.out. These contain dispersivity values for the three principal directions and the number of particles that reached each collection plane. These files are written as ascii files with nine values, corresponding to each of the nine collection planes written on each line. The new values for each realization are appended to the bottom of the existing file, such that at the completion of the run_50_back.ksh and run_50_feat.ksh shells, each of these files has a total of 50 lines. The dispersion calculations were run on a Sun UltraSparc UNIX workstation. The directory path for the dispersion calculations on the CD-ROM is: /AMR/FEHM/track/ and results of the dispersion calculations are stored in: /AMR/FEHM/track/RESULTS/ 4 INPUTS 4.1 DATA AND PARAMETERS No data were obtained from the TDMS for this AMR. All corroborative information necessary for the creation of the models in this AMR is derived from references. These references are cited Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 13 of 60 in the appropriate places in the text. A list of these three references and the corroborative information obtained from them is shown in Table 1 Table 1. Input Information Sources Data Set Data Description Ascession Number Fracture Length (Sweetkind and Williams-Stroud, 1996) Fracture Length Distributions TIC: 239254 Fracture Orientation (Sweetkind and Williams-Stroud, 1996) Fracture Orientation Distributions TIC: 239254 Fracture Permeability (CRWMS M&O, 1998a) Fracture Permeability Distributions MOL.19981008.0008 Feature Spacing (CRWMS M&O, 1999c) Feature Spacing Distributions Pending 4.2 CRITERIA There are no criteria at this time. 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)). 4.3 CODES AND STANDARDS There are currently no applicable codes or standards that pertain to this analysis. 5 ASSUMPTIONS There are several significant assumptions made in this study regarding the representation of fracture permeability within a numerical flow model. 5.1 EFFECTIVE CONTINUUM REPRESENTATION Numerical modeling of fracture properties is done in one of two ways: discrete fracture models or effective continuum models. The effective continuum model used in this study is a stochastic continuum model as multiple realizations of the permeability field are created. As used in this report, these terms apply solely to the permeability model and do not imply any connotations for the transport model. The transport used in this study is only streamline particle tracking in a single-porosity medium. Discrete fracture models represent each fracture as a distinct object within the modeling domain. While the discrete fracture approach retains the discrete nature of the observed fractures within the resulting model, the computational burden of calculating a flow solution on a discrete fracture model becomes extremely large for even a relatively simple fracture model. For this reason, and the fact that the Yucca Mountain Project does not have a qualified discrete fracture modeling software package, the effective continuum representation of fracture permeability is used in this study. Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 14 of 60 The main assumption inherent in using the effective continuum approach is that the permeability field arising from the discrete fractures can be accurately represented by a numerical model composed of continuum gridblocks. Under the continuum approach, the fractures do not appear as discrete entities, but rather they contribute to an effective permeability at each gridblock. The resulting models do not exhibit discrete fractures, but rather areas of higher and lower permeability. In order to retain as much of the discrete fracture behavior as possible, the numerical gridblocks are kept to a minimum size. For the work here, the gridblocks are 4 x 4 x 2 meters. This gridblock size is not inconsistent with the estimated widths of the fracture zones observed at Yucca Mountain (CRWMS M&O, 1999c). This assumption is used in Section 6.1 (see paragraphs 1 and 2) and also throughout Section 6.1.2. This assumption does not require further confirmation. 5.2 BOUNDARY CONDITIONS In this study, tracer transport is analyzed within a single site-scale model numerical gridblock that is 500 x 500 x 50 meters. This gridblock is part of a larger model and the boundary conditions for this gridblock within the site-scale model would be constant fluxes into and out of the neighboring gridblocks. However, for this study on a single site-scale gridblock, the boundary conditions are set to fixed head values on the two ends of the domain and to no-flow boundaries on the sides of the model domain. The head difference across the model domain is based on the average value of the gradient along the direction of the most likely flowpath in the SZ site-scale model. This gradient is 2.9x10-04 and is within the range of gradients presented by Luckey, et al. 1996, p. 27). The inherent assumption here is that these imposed boundary conditions do not adversely affect the dispersion calculations. This assumption is not explicitly restated in Section 6, but is used in Section 6.2.2. This assumption does not require further confirmation. 5.3 SOURCE SIZE The spatial distribution of the potential radionuclide plume within the unsaturated zone is a complex function of the local saturation, the canister failure history and the heterogeneous properties of the fracture system. This uncertainty in the unsaturated zone source history and transport properties leads to uncertainty in the shape of the potential radionuclide plume within the saturated zone. With regards to source size, it is assumed that each of two different conceptual models is a valid representation of the source geometry. The first conceptual model of the source is that used in the recent Total System Performance Assessment-Viability Assessment (TSPA-VA) and the second model is considered to be the most conservative case. The first conceptual model is that of planar source geometry. This model is adapted from the streamtube model of the contaminant plume employed in TSPA-VA (CRWMS M&O, 1998a, p. 8-40 and Figure 8-22) in which the radionuclides enter the saturated zone in one of six different streamtubes from the unsaturated zone. Each streamtube was associated with a different zone of the repository footprint (CRWMS M&O, 1998a, p. 8-40). The planar source implementation in this study assumes that the source zone is 400 meters wide by 10 meters high and is centered within the model domain. Because this study employs a heterogeneous representation of the permeability field, the particle release points are not uniformly distributed (the perfectly mixed Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 15 of 60 scenario used in TSPA-VA) but are distributed proportional to the groundwater flux through the 400 x 10 meter plane. The second conceptual model of the source is that of a point release. This source model corresponds to a single canister release in the unsaturated zone or, in the most conservative case, strong channeling (focusing) of flow within the unsaturated zone from a number of canister releases. It is assumed that under this scenario the potential radionuclide plume enters the saturated zone as a spatially compact source and therefore, under this scenario, all of the particle release points are within a single numerical gridblock. The gridblock chosen for the particle release points is that with the highest groundwater flux that also lies within the 400 x 10 meter plane used as the planar source. This choice produces the most conservative estimates of plume spreading. The assumed source geometries are based on the current understanding of UZ-SZ transport coupling, and the second of these two assumed geometries represents the most conservative case for this coupling. The two different source sizes are employed in Section 6.2.2. This assumption does not require further confirmation. 5.4 GEOSTATISTICAL REPRESENTATION It is assumed that geostatistical simulation can be used to adequately model the fracture lengths and the spacings between flowing intervals. This assumption is tested against the resulting models in Section 6.1.4. For this work, exponential variogram models are employed with the effective range set to be three times the desired feature, or spacing, length. By definition, for the exponential model, the effective range is three times the integral scale, and the integral scale is equal to the mean feature length (Deutsch and Journel, 1998, p. 25). This assumption is employed in Sections 6.1.4, 6.1.2.1 and 6.1.2.2. This assumption does not require further confirmation 5.5 ORIENTATION OF GROUNDWATER GRADIENT In this study, it is assumed that the higher permeability features created by the fracture patterns control the direction of groundwater flow (i.e., the hydraulic gradient is aligned with the principal direction of the permeability anisotropy ellipse). This assumption is supported by observations of the principal fracture orientation being generally north-south (Sweetkind and Williams-Stroud, 1996, p. 52) and the hydraulic gradient south of the repository to also be oriented in a north-south direction (see Luckey, et al. 1996). This assumption is used in Section 6.2.1, paragraph 1. This assumption does not require further confirmation. 5.6 CALCULATION OF DISPERSIVITY The analytical solution used here to compare to the numerical results is based on the assumption that the solute plume can be described by a Gaussian distribution of concentration in the three principal directions. This assumption does not hold up in the numerical results presented here and the differences are discussed in Section 6.2.3. This portion of the assumption underlies the analytical solution of equations 1, 2, 3 and 4 in Section 6.2.1. Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 16 of 60 For the calculations done here, the Gaussian assumption is also invoked to describe the distribution of log transformed travel times for calculation of the dispersivities. This assumption is made primarily to use the log-transform as a smoothing technique to decrease the effects of statistical outliers in the time or displacement distributions. Examination of the log transformed travel times showed the Gaussian distribution to be an adequate approximation. Also, by using a large number of particles (4000) the deviations from the Gaussian distribution are limited. This portion of the assumption is used in Equations 5, 6a and 6b in Section 6.2.1. This assumption does not require further confirmation. 5.7 FRACTURE FLOW For the models contructed here, it is assumed tha the longer length fractures dominate the flow sytstem and especiallly the relatively faster flow paths within the fractured rocks. This assumption is based on observations of rapid pressure responses in outlying wells during the Chole testing (Geldon et al., 1997) This assumption is also based on providing the most conservative (fastest) estimates of groundwater transport. Simple geometry shows that the larger fractures are better connected across the domain. If these fractures are also assigned relatively higher permeability values, flow will be faster than in the less connected and lower permeability fractures. This assumption is used in Section 6.1.2.1. This assumption does not require further confirmation. 6 ANALYSIS 6.1 REPRESENTATION OF FRACTURE PERMEABILITY The purpose of this portion of the modeling exercise is to create stochastic images of the fracture permeability within a 500 x 500 x 50 meter domain. These fracture permeability fields are used as input to a groundwater flow and particle tracking exercise being conducted to better define the amounts of dispersion at the sub gridblock scale for the Yucca Mountain Site Characterization Project (YMP) SZ site-scale model. The fracture permeability fields are unconditioned in the sense that there are no point measurements of permeability that must be reproduced within each realization. However, there are specific univariate distributions of fracture attributes observed on discrete fractures, as well as prevailing conceptual models that must be represented in the final fracture permeability models. This analysis explicitly applies only to that part of the SZ flow path that is located within the fractured Tertiary volcanic rocks in the SZ down-gradient from beneath the potential repository. It is proposed that the heterogeneous fields created through the techniques outlined in this report provide a meaningful representation of the uncertainty in particle tracking results arising from uncertainty in the spatial distribution of the fracture permeability. In contrast to discrete fracture models, which represent fractures as discrete objects within an impermeable matrix, an alternative technique for representing fracture permeability within a groundwater flow model is the stochastic continuum model (SCM). Stochastic continuum Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 17 of 60 models have been proposed as a means of representing the spatially heterogeneous nature of fracture permeability at the scale of continuum gridblocks (see Hsieh, 1998, pp. 336-338). In general, the discrete behavior of individual fractures is represented by an average property over a representative elementary volume (REV). The SCM approach has been used previously to model fracture systems at experimental facilities for nuclear waste programs in the U.S. (see Altman et al. 1996) and in Sweden (Follin and Thunvik, 1994; Tsang et al. 1996). In these studies, geostatistical simulation algorithms were used to model continuum properties of discrete fracture systems at the scale of numerical flow and transport model gridblocks (1’s to 10’s of meters). The geostatistical realizations of these properties were then used as input to flow and transport models to predict the advective travel time of particles through the heterogeneous systems. The current project also employs a SCM approach to modeling fracture permeability. Two techniques for creation of fracture permeability models are described. These techniques have been used to populate the three-dimensional model domain representative of a single gridblock within the site-scale saturated zone (SZ) model. The background pattern of grid-block transmissivity across the three-dimensional domain is defined using an indicator geostatistical simulation process. This spatial simulation algorithm populates the domain with a continuous distribution of log10 fracture permeability having specified ranges of spatial correlation. The distribution is non-parametric such that at predefined points on the cumulative distribution, the amount and anisotropy of the spatial correlation can be varied. This variability allows the higher permeability portion of the distribution to have an anisotropic and preferentially aligned orientation such as that observed in fractures and flowing features at Yucca Mountain. Additionally, some larger features that may dominate the flow and transport for a region of the saturated zone are simulated using a boolean simulation algorithm. The features produced by the boolean simulation are added to the background permeability pattern modeled by the indicator geostatistical simulations to produce the enhanced models of permeability. 6.1.1 Available Data The available data, listed in table 1 of Section 4.1, exist as both direct measurements of fracture permeability at the Yucca Mountain site as well as observations of discrete fractures made on outcrops and in the Exploratory Studies Facility (ESF). Two conceptual models of fracture permeability are developed (Section 6.1.2) from the interpretations of hydraulic testing done in both the saturated and unsaturated (air testing) zones as well as from observations of the fractures made in boreholes and on outcrops and within the ESF. The various sources of data are reviewed below. 6.1.1.1 Permeability Testing The dominant mechanism for water movement in the fractured volcanic rocks at Yucca Mountain is believed to be fracture flow, and consequently the origin, characteristics and occurrence of fractures have been the subject of several studies (Sonnenthal, et al. 1997; Sweetkind and Williams-Stroud, 1996). Hydraulic tests have been conducted using a variety of test methods and interpretive models. Single-well, double-packer tests have been conducted throughout the aquifer, and multi-well interference tests have been conducted in a number of the Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 18 of 60 identified flowing intervals (CRWMS M&O, 1998a, Figure 8-27). The log permeability histograms (Figure 1) of these tests overlap, suggesting that some of the single-well test results represent the more conductive flowing intervals (CRWMS M&O, 1998a; Figure 8-27). The natural-log variance of the data set in Figure 1 is 9.46 The minimum log permeability of the multi-well tests suggests that the log permeability threshold for flowing features is approximately –12.2 m2. All the permeability data shown in Figure 1 yields an empirical cumulative distribution (cdf) typical of a weakly bimodal distribution (Figure 2), similar to that used in the UZ studies of Tsang (1997). A specific set of tests has been conducted at the C-well complex southeast of the proposed repository. Interpretations of these tests and borehole televiewer logs by Geldon (1996) indicate that the highly transmissive zones are due to discrete zones of enhanced fracturing. Hydraulic testing at the C-well complex has also indicated the presence of large-scale (up to several kilometers) high-permeability features. In the vicinity of the C-well complex, these features are associated with observed faults (Geldon et al. 1997) 6.1.1.2 Borehole Flowmeter Testing A number of studies in fractured rock aquifers have shown that groundwater flow does not occur uniformly within all fractures, but instead often occurs in just a small fraction of the available fractures (see Hsieh, 1998, pp.348-350; Tsang, et al. 1996). More generally, the aquifer material may contain many different types of deformation structures: fractures, joints, shear zones, etc. that may be formed by different genetic processes (e.g., cooling or tectonic deformation). Given the various genetic histories, these different structures are often referred to using the general term “features”. The features that conduct the largest volume of groundwater are of most importance for groundwater transport. These features are often referred to as “flowing features”, “flowing intervals”, or “water conducting features”. A recent analysis of fracture and flowing feature data (CRWMS M&O, 1999c) provides information on the spacing between flowing features and fractures. Borehole flowmeter surveys indicate that discrete intervals of the aquifer have flow rates that are noticeably greater than the remainder of the aquifer. The spacing and thickness of these flowing intervals were analyzed (CRWMS M&O, 1999c), and it was found that the dip-corrected geometric mean flowing interval spacing is 19.7 m (Figure 3). It should be noted that there is no direct measurement of the length, width, structure or origin of the “flowing features” that create these flowing intervals. Only the spacing between these features is measured. 6.1.1.3 Fracture Observations Mapping studies in boreholes, and the ESF suggest a dense, well-connected set of fractures, with approximately 89% of the fractures in the range of 0.3 to 3m in length, and approximately 1% between 10 to 34m. Fracture spacings tend to be on the same order as fracture lengths, i.e., fractures of 10m length tend to be spaced 10m apart. Geometrically, this suggests that the fractures are well-connected over a range of length scales, contributing to a scale-dependent permeability (Sonnenthal et al. 1997, p 7-26 and p 7-34). Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 19 of 60 Figure 1. Histogram of permeability measurements made in the middle volcanic aquifer including the distribution (solid line) derived from an expert elicitation (CRWMS M&O, 1998b). The dotted and dashed lines are the log-normal distributions fit to the single-hole and multi-hole data respectively (CRWMS M&O, 1998a). This histogram is shown as a cumulative distribution function in Figure 2 (obtained from CRWMS M&O, 1998a) Single-Hole Multi-Hole PDF Experts Log10 (Permeability in m2) No. of Observations Probability 0.0 0.1 0.2 0.3 0.4 0.5 0 10 20 30 40 -20 -18 -16 -14 -12 -10 -8 Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 20 of 60 Figure 2. Cumulative distribution of observed log10 permeability measurements (CRWMS M&O, 1998a) as shown in Figure 1 and the model used to describe the cumulative distribution within the geostatistical simulation algorithm. The fractures appear to be of two origins, each with specific geometries and spacings. The first fracture type is cooling joints, which tend to be relatively smooth, subvertical fractures of approximately 4m length, oriented north-south, with a high proportion of blind terminations (Sonnenthal et al. 1997; Tsang, 1997). The second type of fracture is tectonic, induced by rock stress. These fractures have a wide range of lengths from the microscopic to many tens of meters. Tectonic fractures also tend to be subvertical and oriented north-south, as result of the reactivation of cooling joints. Borehole data and line surveys indicate that, after correcting for dip, the geometric mean fracture spacing is 0.25m for all fractures (Table 6; CRWMS M&O, 1999c). Sweetkind and Williams-Stroud (1996) conducted a synthesis of the fracture data available at Yucca Mountain. The results of this synthesis can be used to help constrain the models of fracture permeability created in this work. Specifically, information regarding fracture length and the variability in fracture orientation can be gained from the work of Sweetkind and Williams-Stroud (1996) and used directly in the creation of the feature-based “enhanced” conceptual model of fracture permeability as discussed below. Log10 Permeability (m2) -20 -18 -16 -14 -12 -10 -8 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 Observed Data Model Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 21 of 60 Figure 3. Distribution of the spacing between flowing features. Note that the X-axis is log scaled. Obtained from CRWMS M&O (1999c). Fracture length observations made by Sweetkind and Williams-Stroud (1996) were obtained on both surface outcrops and along a scanline in the ESF. The observations of surface length are truncated on both the short and long ends of the distribution. The truncation on the short end is necessary to make the mapping practical and the truncation on the long end is due to the finite extent of the outcrop exposure. Fracture mapping of the surface exposures employed a minimum trace length (truncation length) of 1.5 meters and the maximum observable length was approximately 15 meters. Given these constraints, the distributions of fracture length for the different outcrops appear to be either exponential or power-law (Sweetkind and Williams-Stroud, 1996, page 52 and Figure 21). Fracture mapping conducted within the ESF was done up to station 18+00 (north ramp). For this portion of the ESF, the detailed line survey was mapped with a truncation length of 0.30 meters. An additional high-resolution map along 60 meters of the ESF was completed with a truncation length 0.15 meters. The results of these two ESF fracture mappings, using considerably lower truncation lengths than were used on the outcrop surveys, also show the fracture lengths to be either exponentially or power-law distributed (Sweetkind and Williams-Stroud, 1996, p. 52 and Figure 22). 1 10 100 1000 Corrected Flowing Interval Spacing Data (m) 0 0.05 0.1 0.15 0.2 0.25 Relative Frequency Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 22 of 60 Orientations of fractures were observed within the non-lithophysal zone of the TSw (Sweetkind and Williams-Stroud, 1996). These orientations show a mainly northeast orientation with a steep dip to the west (Sweetkind and Williams-Stroud, 1996, Figure 17). Other orientations describing a significant proportion of the measured fractures are north and north-northwest. 6.1.2 Conceptual Models Two conceptual models of fracture permeability have been constructed in this analysis. The first model is referred to as the “background” model and the second model is the “enhanced” model. It is noted that the enhanced model is an extension of the background model. The background model is created solely with indicator geostatistical simulation to produce a series of high permeability features with the same orientation, anisotropy and length characteristics. This background model is based on a conceptualization of a single set of fractures as observed at Yucca Mountain. This model of fracture permeability is created such that the long-axis of the anisotropic permeability ellipse is aligned parallel to the imposed hydraulic gradient used in the flow and particle tracking model (Section 6.2). This modeling decision is made based on the assumption (Section 5.3) that the overall fracture permeability resulting from these high permeability features controls the direction of groundwater flow. Additionally, at Yucca Mountain, the primary fracture orientation is approximately north-south, and for the area south of the repository, this orientation is roughly parallel to the hydraulic gradient (see CRWMS M&O, 1998a, Figure 8-21). The enhanced model uses the initial permeability structure created for the background model and then overlays a set of discrete features onto the existing background model. These discrete features represent the large features observed in the field, and the larger unobserved features inherent in the suggested power-law distribution of fracture lengths, that may control the local flow patterns. Examples of such features are the high permeability connections observed between the C-well complex and boreholes USW H-4 and UE-25 WT#14 (Geldon et al. 1997, Table 9, p “Hydraulic Tests-37”). The middle volcanic aquifer consists of welded tuffs that include all or parts of the Prow Pass, Bullfrog and Tram formations of the Crater Flat Group (CRWMS M&O, 1998a, Table 8-11). This work is based on data for the middle volcanic aquifer and the conceptual models of fracture permeability are also created for the middle volcanic aquifer. However, these models are valid interpretations of any moderately to densely welded tuff in the saturated zone at Yucca Mountain. This interpretation is supported by results (CRWMS M&O, 1999c, Tables 8 and 9) that show the spacing of flowing features is not a function of the specific rock type or hydrogeologic unit. 6.1.2.1 Background Permeability Model The field observations and lithology of the middle volcanic aquifer indicate a medium dominated by fracture flow with strong anisotropy (major axis aligned north-south) and subvertical dip directions. The fracturing appears to be sufficiently dense and interconnected that it exhibits a weakly scale-dependent behavior, similar to a stochastic continuum. Relatively large flowing Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 23 of 60 features have been identified, and hydraulic tests define a bimodal distribution of the log permeability for the flowing features and the intervening rock. For the purposes of continuum modeling, these characteristics suggest simulating the middle volcanic aquifer as an anisotropic permeability field, with separate classes of spatial continuity to represent the flowing features and the surrounding rock, respectively. One method of geostatistical simulation, the indicator approach, is well adapted to such a bimodal simulation. The indicator approach uses the empirical cdf of the data to develop a discrete, nonparametric estimate of the true cdf of the simulated medium. Each cutoff level of the input discrete cdf has a separate spatial correlation function represented by a variogram, whose range represents the continuity of features at that cutoff level. For a preferentially fractured media, the major axis of variogram anisotropy corresponds to the fracture plane. This study uses the sisim (version 2.0) routine, from the GSLIB package (Deutsch and Journel, 1998) to create the background fields of log10 permeability. The input discrete cdf for log permeability is simply inferred from the pooled single-well and multiwell hydraulic test data (Figure 1). Cutoff values are chosen from the empirical cdf (Figure 2) to both realistically represent the observed data and create numerically stable simulations. The minimum of log10 permeability is –18.2 m2, and the maximum is –10.0 m2. The minimum log10 permeability of the multi-well hydraulic tests, –12.2 m2, is used as a cutoff value in the discrete cdf and denotes the threshold between background and feature permeability. Variogram models for each cutoff are inferred based on the conceptual model. In the ideal case, experimental variograms of the hydraulic test data could be used to infer variogram models, but there are too few tests to allow such direct inference. Furthermore, the relationship between discrete fracture characteristics and permeability are not well understood. As an indirect approach, this study uses the observed fracture lengths and flowing interval spacings to set variogram ranges, using the relationship that the integral scale of the exponential variogram model is equal to the average length of simulated features. For cutoffs below the critical value of –12.2 log10 m2, flow is assumed (see Section 5.7) to be dominated by the permeability resulting from the majority of fractures, whose length is 3 to 4 m or less (Sweetkind and Williams-Stroud, 1996). Thus in the plane of the fracture, the assumed variogram integral scale is 4 m (ie., an exponential variogram with effective range of 12 m). The thickness of such fractures (in the X-direction) is on the order of millimeters, suggesting that adjacent blocks of a finite element should not be correlated. This will be achieved if the effective range perpendicular to the fracture plane is less than the target grid scale of 4 m. Thus, below –12.2 log10 m2, all cutoffs will have an exponential variogram with 3:12:12 m (X,Y,Z) effective ranges in the east:north:vertical directions. The thresholds and the variogram parameters for the permeability model are shown in Table 2. The permeabilities above –12.2 log10 m2 are assumed to represent the flowing features. The lengths of these features are unmeasured, but, conservatively , it can be argued that they are comprised of higher permeability fractures, and as such, are planar features oriented north-south (Y-direction) with subvertical dips, with spacing approximately equal to length. This study uses an integral scale of 20m in the plane of the flowing features (i.e., an exponential model with Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 24 of 60 effective range of 60m). This integral-scale length is based on the interpretation of an exponential or power-law distribution to define the fracture lengths (Sweetkind and Williams- Stroud, 1996, p. 52 and Figure 21). Although distributions were not fit to the observed data, the choice of a 20 meter integral scale corresponds to the upper limit of these suggested distributions based on the data shown in Figure 21 of Sweetkind and Williams-Stroud (1996). The thickness of the measured flowing features is 8.3m, corresponding to an exponential variogram with an effective range of 25m perpendicular to the plane of the flowing feature. That is, the permeabilities –12.2 log10 m2 and above will have an exponential variogram with 25:60:60 m effective ranges in the east:north:vertical directions. Three views of an example background fracture permeability realization are shown in Figures 4 and 5. One limitation of the indicator method is that the frequency and size of features are related through the variogram, and thus the variogram range perpendicular to the plane of the flowing features may require calibration. This study thus calibrates the east-west variogram range for cutoffs of –12.2 log10 m2 and above by post-simulation audits of the permeability fields to verify that the flowing interval spacing is similar to that observed in the field. Table 2. Parameters describing the indicator variogram models used to build the background fracture permeability simulations. An exponential model is specified for all variograms. Threshold Permeability (log10 (m2)) cdf value at Threshold Nugget Sill Effective Range (Y-dir.) (meters) Anisotropy (X:Y:Z) -18.2 0.00 NA Minimum NA NA -15.7 0.06 0.0 0.056 12 3:12:12 -15 0.18 0.0 0.148 12 3:12:12 -14.5 0.30 0.0 0.210 12 3:12:12 -13.7 0.50 0.0 0.250 12 3:12:12 -13 0.65 0.0 0.228 12 3:12:12 -12.2 0.75 0.0 0.182 60 25:60:60 -11.8 0.85 0.0 0.131 60 25:60:60 -11.6 0.92 0.0 0.074 60 25:60:60 -10 1.00 NA Maximum NA NA Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 25 of 60 Figure 4. Example realization of the background fracture permeability conceptual model. The color legend indicates log10 permeability in m2. 0 100 200 300 400 500 Easting 400 300 200 100 0 Northing -16 -15 -14 -13 -12 -11 (m) (m) Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 26 of 60 Figure 5. Two images of a horizontal slice (layer 12) through the example background fracture permeability model. The upper image shows the full distribution of log10 permeability (m2). The lower image shows the locations of those features with log10 permeability greater than or equal to –12.2 m2. 0.0 100.0 200.0 300.0 400.0 500.0 0.0 100.0 200.0 300.0 400.0 500.0 Easting (m) Northing (m) -16 -15 -14 -13 -12 -11 0.0 100.0 200.0 300.0 400.0 500.0 0.0 100.0 200.0 300.0 400.0 500.0 Easting (m) Northing (m) Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 27 of 60 6.1.2.2 Enhanced Permeability Model Boolean simulation is defined as the superimposition of independent (non-interacting) random objects at random locations within the simulation domain (Schmitt and Beucher, 1997, p. 201). In contrast to true geostatistical models, the boolean model relies on complete spatial independence. In the context of fracture permeability modeling, boolean models can be used to place high permeability features at random locations within a domain. The random placement of these objects is derived from the often-used conceptualization of fracture locations as a spatial Poisson process. The shapes of these features need not be random and, for the case of fracture modeling, are often considered to be disks with radii and orientations drawn from specified distributions. In this current analysis, the boolean models are created with the code ellip2_MC (version 1.0). This code is an extension of the ellipsim code developed by Deutsch and Journel (1998) and allows the random placement of disks with radii drawn from user specified distributions including both exponential and power-law models. The orientations of the disks are also drawn from user specified distributions. The distributions for the orientations are either normal or triangular. For this exercise, the effect of a small number of domain spanning features is examined by simulating disks with a length that is equal to or larger than the domain length. The centers of the disks can lie outside the simulation domain such that only portions of the disks appear within the domain (see left front corner of Figure 6). The orientation of these features is drawn from a normal distribution with a mean orientation of 0 degrees (North) and a standard deviation of 15 degrees. The boolean simulation is constrained to comprise a total of two percent of the simulation domain. The resulting boolean simulations are combined with the background fracture permeability simulations using the routine combine. This combination is done by simply mapping the disks onto the background permeability simulation. The mapping is done such that the permeability assigned to the disk is retained at all points within the disk, regardless of the background permeability at that location. A constant log10 permeability of –11 is assigned to the disks. This permeability corresponds, approximately, to the 95th percentile of the measured permeability values as shown in Figure 2. A resulting combination of indicator geostatistical and boolean models is shown in Figure 7. Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 28 of 60 Figure 6. Example of large-scale features created using boolean simulation. These features are assigned a constant log10 permeability of –11.0 m2 and combined with an indicator geostatistical simulation such as that shown in Figure 4 to produce the final “feature” model of fracture permeability. Figure 7. Example realization of the feature-based conceptual model of fracture permeability. This model is created by combining the geostatistical simulation of background permeability in Figure 4 and the boolean simulation of features in Figure 6. 50 25 0 Elevation (m) 0 100 200 300 400 500 Easting (m) 500 400 300 200 100 0 Northing (m) -16 -15 -14 -13 -12 -11 50 25 0 Elevation (m) 0 100 200 300 400 500 Easting (m) 500 400 300 200 100 0 Northing (m) Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 29 of 60 6.1.3 Implementation Issues The size and discretization of the permeability model domain are set to be coincident with that of the FEHM flow and transport model. Several post-processing checks on the heterogeneous models are conducted prior to using them as input to the flow and particle-tracking model. 6.1.3.1 Simulation Domain The models of fracture permeability are created within a three-dimensional domain that is representative of a single gridblock (element) within the saturated-zone site-scale model. The dimensions of the domain are 500 x 500 x 50 meters. The size of this domain is near that of the finest resolution gridblock within the SZ site-scale model. For this study of the effects of sub gridblock heterogeneity on flow and transport parameters, the domain is discretized into 4 x 4 x 2 meter gridblocks. The number of elements in the X, Y and Z directions is 125, 125 and 25 respectively for a total of 390,625 elements. Due to the node-based structure of FEHM, the final dimensions of the flow model domain are actually 496.0 x 496.0 x 48.0 meters. This level of discretization within the domain represents a compromise between trying to capture the finescale heterogeneity resulting from variations in the permeability and connectivity of individual discrete fractures and still being able to run the multiple flow and transport models using FEHM in a reasonable time frame. This is the same discretization and domain size shown in the figures of this report. 6.1.3.2 Numerical Implementation A number of steps using different routines are required to complete the permeability modeling process. These steps are described below. The background permeability models are created using the GSLIB routine sisim. These models are created on a Pentium PC and then transferred to the NWMP UNIX network via ftp. A total of 50 realizations of the background permeability are created using sisim. The output files from sisim are named as: sisim.#.out, where “#” is the realization number. To facilitate efficient file transfer, these realizations are created in batches of 10 (each input file creates 10 realizations and stores them in the same output file). Electronic copies of the five input files are included on the CD-ROM accompanying this report. The directory containing the input and output files for the background permeability models is on the CD-ROM as: /AMR/SISIM/SET1/ The enhanced permeability models are created by running the UNIX shell file cr8_50_feat.ksh. These calculations are done on a Sun UltraSparc UNIX workstation. A listing of this shell file is included as Attachment 3. This shell file does the following: • uncompresses the corresponding background permeability realization and copies it to the file sisim.out in its current directory (/AMR/SISIM/SET1/ on the CD-ROM) Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 30 of 60 • runs the routine ellip2_MC routine to create a domain containing several discrete disk-shaped features. The input file for this routine is named ellip2.par. The output files from ellip2_MC are named ellip2.out and they are used as input to the routine combine. The routine ellip2_MC also creates a debugging file with the naming convention: ellip2.*.dbg, where the * indicates the realization number. • combines the background permeability model with the disk shaped features using the routine combine. The inputs to combine are the sisim.out and ellip2.out files discussed above. The output files from combine are named combine.out and used as input to the gs2fehm and gamfil routines. • checks the spacing of high permeability features (those with log10 permeability greater than –12.2 m2) using the routine gamfil. The input into gamfil is the combine.out file. The output from gamfil is named gam.*.res where the * indicates the realization number. All of the gam.*.res files are combined together using the UNIX cat command to create the gam_all.res file. This file is the source of the distributions shown in Figures 8 and 9. The source code for the gamfil routine is Attachment 9. • creates a new input file (ellip2.par) for the ellip2_MC routine by changing the random number generator seed in the old input file. This change of seed is accomplished by running the routine modgeom. The source for the modgeom routine is listed as Attachment 10. The above steps are done in a loop over all 50 realizations. Note that the feature spacing is also checked on the background permeability models in the same way as described above for the enhanced permeability models. The shell file to accomplish checking of the background permeability feature spacing is run_50.ksh in the directory: /AMR/SISIM/SET1/ on the CDROM. The source code for this shell is in Attachment 2. 6.1.4 Model Validation Model validation is completed to check the reproduction of the measured data (permeability distribution and feature spacing) in the resulting permeability files. Additionally, the reproduction of the assumed correlation lengths (variogram models) in the resulting permeability models are also checked. The criteria used in the validation are: 1) On average, the cdfs of the permeability realizations must match the cdf of the input permeability data. The quantitative criterion used here is to check the mean cdf value across all 50 realizations at each of the 8 threshold values. The mean threshold value from the 50 realizations must be within +/- 10% of the input threshold value. 2) The mean high permeability feature spacing is calculated for each of the 50 realizations and compared to the distribution of flowing interval spacings in Figure 3. The validation criterion is that the distribution of mean spacings must fall within the distribution of observed spacings as shown in Figure 3. 3) The variograms are calculated across each realization and are compared to the input variogram model. Unlike the permeability distribution and the feature spacing Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 31 of 60 distribution there are no direct observations of the variogram; the validation can only be made against the assumed variogram models used as input. The validation criterion is that the shape of the model variogram be reproduced within the extent possible given the ergodic hypothesis constraints as discussed below. 6.1.4.1 cdf Reproduction The reproduction of the measured distribution of permeability used as input to the fracture permeability models (Figures 1 and 2) is checked. This check is made by determining the value of the cdf at each of the eight thresholds defined in Table 1 for each of the 50 realizations for both conceptual models of fracture permeability. The results of these calculations are shown in Table 3 and Table 4. Table 3. Results of 50 realizations of the background permeability model compared to the input distribution at eight permeability thresholds. Threshold Log10 (m2) Input cdf Mean cdf Minimum cdf Maximum cdf -15.7 0.06 0.060 0.056 0.063 -15 0.18 0.178 0.173 0.183 -14.5 0.30 0.295 0.289 0.301 -13.7 0.50 0.481 0.472 0.488 -13 0.65 0.616 0.606 0.625 -12.2 0.75 0.822 0.810 0.834 -11.8 0.85 0.888 0.879 0.896 -11.6 0.92 0.938 0.931 0.947 Table 4. Results of 50 realizations for the enhanced conceptual model compared to the input distribution at eight permeability thresholds. Threshold Log10 (m2) Input cdf Mean cdf Minimum cdf Maximum cdf -15.7 0.06 0.059 0.054 0.062 -15 0.18 0.174 0.168 0.179 -14.5 0.30 0.287 0.281 0.294 -13.7 0.50 0.470 0.461 0.478 -13 0.65 0.600 0.588 0.612 -12.2 0.75 0.802 0.783 0.816 -11.8 0.85 0.865 0.854 0.878 -11.6 0.92 0.914 0.904 0.927 In general, the resulting background and enhanced permeability models reproduce the measured distributions of permeability. The resulting models do overpredict and underpredict some of the threshold cdf values; however, these discrepancies are always within 10% or less of the input values. The resulting models are deemed to be adequate representations of the input permeability cdf. 6.1.4.2 High Permeability Feature Spacing The resulting realizations of enhanced fracture permeability were checked to determine the spacing of high permeability features. A randomly located transect was drawn through each simulated field parallel to the X-axis (Easting axis in Figures 4 to 7). Along each transect, the simulated permeability values were compared against the assumed lower limit of permeability Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 32 of 60 for the flowing features (-12.2 log10 m2). If a simulated permeability exceeds that threshold, then the sampling routine counts that value as a flowing feature and continues on to the next simulated value. Contiguous values exceeding the flowing feature minimum are counted as a single flowing feature, and the statistics for the spacings and widths of such features are accumulated for each transect. Results of the spacing and width calculations are presented in Figure 8 and 9 as distributions calculated across all 50 realizations. Results of the transect sampling shown in the upper images of Figures 8 and 9 can be directly compared with the distribution of spacings between flowing features in Figure 3. Examination of Figures 8 and 9 shows that the geostatistical realizations produce average feature spacings that are within the distribution of measured feature spacings. These models are deemed to be valid representations of the observed feature spacing. The available data on flowing features do not indicate the width of these features. They could represent single discrete fractures or enhanced flow within fracture zones. The continuum model representation of these features is limited to a minimum feature width equal to that of the numerical gridblocks (4.0 meters). The distributions in the lower images of Figures 8 and 9 show that on average the mean feature width is approximately five meters. These mean values can be interpreted as every fourth feature being two gridblocks wide (8 meters) with the remainder of the features being one gridblock (4 meters) in width. 6.1.4.3 Variogram Reproduction Reproduction of the input variogram models is checked across the ensemble of realizations created for the background fracture permeability model. This validation step is conducted by calculating the variogram on each realization in the three principal directions (X,Y,Z) and comparing the results to the input variogram model. The variogram values calculated across the ensemble of realizations for any specific separation distance (or lag) will form a distribution. For the majority of the lag spacings, this distribution should capture the specified input model. The checks on the reproduction of the specified variogram model in the resulting permeability fields are shown in Figures 10, 11 and 12. Variograms are checked at two of the eight thresholds: threshold 3 at –14.5 log10 m2 and threshold 6 at –12.2 log10 m2. Variogram reproduction is excellent at the lower threshold (upper images, Figures 10, 11 and 12). The variograms at threshold 6 show the distinct effect of the finite domain size. Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 33 of 60 Figure 8. Distributions of high permeability feature spacing (upper image) and feature width (lower image) for the background fracture permeability models. The values shown are the mean values from one sampling transect for each of the 50 realizations. The heavy dashed line indicates the mean value of the distribution and the thin dashed line indicates the median value of the distribution (overlapping in these plots). 1.0 10.0 100.0 1000.0 0.0 0.05 0.10 0.15 0.20 Mean Feature Spacing (m) Relative Frequency 1.0 10.0 100.0 0.0 0.05 0.10 0.15 0.20 Mean Feature Width (m) Relative Frequency Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 34 of 60 Figure 9. Distributions of high permeability feature spacing (upper image) and feature width (lower image) for the enhanced permeability models. The values shown are the mean values from one sampling transect for each of the 50 realizations. The heavy dashed line indicates the mean value of the distribution and the thin dashed line indicates the median value of the distribution (overlapping in the upper plot). 1.0 10.0 100.0 1000.0 0.0 0.05 0.10 0.15 0.20 0.25 Mean Feature Spacing (m) Relative Frequency 1.0 10.0 100.0 0.0 0.05 0.10 0.15 0.20 0.25 Mean Feature Width (m) Relative Frequency Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 35 of 60 Figure 10. Results of east-west (X-direction) variograms calculated on all 50 realizations of the background fracture permeability model and compared to the input variogram model. The upper image shows the results for the –14.5 log10 m2 permeability threshold and the lower image shows the results for the –12.2 log10 m2 permeability threshold. The gray zones in both images indicate the range of variability in the variance of the realizations. Legend Variance Model Minimum Maximum Mean 0.0 10.0 20.0 30.0 40.0 50.0 0.00 0.05 0.10 0.15 0.20 0.25 Distance (meters) Indicator Semi Variance Legend Variance Model Minimum Maximum Mean 0.0 20.0 40.0 60.0 80.0 100.0 0.00 0.05 0.10 0.15 0.20 Distance (meters) Indicator Semi Variance Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 36 of 60 Figure 11. Results of north-south (Y-direction) variograms calculated on all 50 realizations of the background fracture permeability model and compared to the input variogram model. The upper image shows the results for the –14.5 log10 m2 permeability threshold and the lower image shows the results for the –12.2 log10 m2 permeability threshold. The gray zones in both images indicate the range of variability in the variance of the realizations. Legend Variance Model Minimum Maximum Mean 0.0 10.0 20.0 30.0 40.0 50.0 0.00 0.05 0.10 0.15 0.20 0.25 Distance (meters) Indicator Semi Variance Legend Variance Model Minimum Maximum Mean 0.0 20.0 40.0 60.0 80.0 100.0 0.00 0.05 0.10 0.15 0.20 Distance (meters) Indicator Semi Variance Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 37 of 60 Figure 12. Results of vertical (Z-direction) variograms calculated on all 50 realizations of the background permeability model and compared to the input variogram model. The upper image shows the results for the –14.5 log10 m2 permeability threshold and the lower image shows the results for the –12.2 log10 m2 permeability threshold. The gray zones in both images indicate the range of variability in the variance of the realizations. Legend Variance Model Minimum Maximum Mean 0.0 10.0 20.0 30.0 40.0 50.0 0.00 0.05 0.10 0.15 0.20 0.25 Distance (meters) Indicator Semi Variance Legend Variance Model Minimum Maximum lag variance Mean 0.0 10.0 20.0 30.0 40.0 50.0 0.00 0.05 0.10 0.15 0.20 Distance (meters) Indicator Semi Variance Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 38 of 60 Central to the theory underlying geostatistics is the ergodic hypothesis that allows the replacement of a spatial variable within a single unbounded domain by an average of the same variable over an ensemble of realizations (see Zhan, 1999, p. 113). In practical terms, the finite domain over which the ensemble of realizations are created needs to be roughly 10 correlation lengths or greater in each dimension to achieve accurate variogram reproduction. For the case considered here, the variograms at threshold 3 are reproduced accurately because the correlation lengths are all small enough (3:12:12 meters in the X, Y and Z directions) relative to the domain size (500 x 500 x 50 meters). The variogram ranges for threshold 12 are: 25:60:60. This violation of the ergodic hypothesis lowers the overall indicator variance of the simulated permeability fields and causes the variograms to fall below the input model (lower images, Figures 10, 11 and 12). The shape of the model variograms is preserved in the resulting realizations, but the sill (total variance) is less than that of the data set due to limitations imposed by the ergodic hypothesis. The shape of the input variogram model is maintained, albeit at a lower variance, thus these realizations are considered to be valid representations of the input geostatistical model. The implication of this lower sill value at the –12.2 m2 threshold is that the variability in the modeled permeability distribution at the higher permeability values is not as large as in the input data set (Figures 1 and 2). However, this decrease in variance is only slight as shown by the close reproduction of the input cdf by the models in the lower rows of Tables 3 and 4. Therefore, these lower sill values will have negligible effects on the dispersivity calculations. 6.2 CALCULATION OF DISPERSIVITIES Dispersion in groundwater transport systems can be defined as the occurrence and evolution of a transition zone between two domains of the fluid phase with different compositions. This definition does not imply any specific process through which this transition zone comes into being, nor does it imply a specific scale for this transition zone. Dispersivity is the rate of change of the second moment of the contaminant plume. Dispersivity is generally calculated along the principal axes of the flow system. Along each principal axis, the dispersivity, with units of length, can be considered as a characteristic length of the transition zone between the two domains of the fluid phase. 6.2.1 Approach For the calculation of the transverse dispersivities, the approach of Wen and Gomez-Hernandez (1998, pp. 144-145) is followed. Wen and Gomez-Hernandez (1998) examined two-dimensional transport in heterogeneous fields with particle tracking. This approach is outlined and expanded to the three-dimensional case below. This work represents the solute plume by a finite set of particles tracked along streamlines. In this situation, the dispersivities are given in terms of the change in second moments, X11, X22, X33, of the particle trajectories with respect to time (after Wen and Gomez-Hernandez, 1998): Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 39 of 60 dt dX U t dt dX U t dt dX U t V H L 33 22 11 2 1 ) ( 2 1 ) ( 2 1 ) ( = = = a a a Where t is time, U is the mean particle velocity, aL is the longitudinal dispersivity, aH is the horizontal transverse dispersivity and aV is the vertical transverse dispersivity. Equations 1a and 1b are taken directly from Wen and Gomez-Hernandez (1998, p. 133). Equation 1c is the same as 1a and 1b, but written for the third (Z) dimension. Apparent macrodispersivities are those determined from a single measurement of the particle trajectories at time t. These macrodispersivites are given, for an instantaneous injection boundary condition, as (after Wen and Gomez-Hernandez, 1998, p. 133): Ut t X t Ut t X t Ut t X t V H L 2 ) ( ) ( 2 ) ( ) ( 2 ) ( ) ( 33 22 11 = = = a a a Several stochastic theories defining macrodispersion have been developed (e.g., Dagan, 1988; Gelhar and Axness, 1983). For transport in a three-dimensional field with a flux weighted source boundary condition, Dagan (1988, eqns 24 and 38) gives the following expressions for the second moments of the plume as a function of normalized time: .. . .. . . .. . . .. . .. . .. . + + - + - = = .. . .. . . .. . . .. . .. . .. . + - + - - = - - ' 2 3 3 2 2 33 22 ' 2 3 2 2 11 ' 1 ' 4 ' 4 ' 4 ' 1 3 1 2 ) ' ( ) ' ( ' 1 1 ' 8 ' 8 ' 4 3 8 ' 2 ) ' ( t Y t Y e t t t t t t X t X e t t t t t t X . s . s (Eq. 1a) (Eq. 1b) (Eq. 1c) (Eq. 2a) (Eq. 2b) (Eq. 2c) (Eq. 3) (Eq. 4) Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 40 of 60 where sY 2 is the variance of the lnK field with isotropic spatial correlation of integral length . (practical range = 3.), and t’= X/. is the dimensionless time, or displacement. Combining Equations 2 a,b and c with Equations 3 and 4 gives expressions for the macrodispersivities as a function of time. Using the analytical solutions of Dagan (1988), both longitudinal and transverse dispersivity are calculated for a range of t’. These results are shown in Figure 13 for natural log transformed permeability variances of 1.0 and 9.46. These values of variance are the upper limit considered in the development of the theory and the variance of the saturated zone permeability data shown in Figure 1, respectively. The longitudinal dispersivity is seen to rise to an asymptotic value after a travel distance of several tens of correlation lengths (longitudinal correlation length = 20 m in this study). The two transverse dispersivities (i.e., horizontal and vertical) are equal in magnitude and initially increase in value and then decrease asymptotically to zero at infinite time. For the numerical modeling done in this study, the longitudinal and transverse dispersivities are calculated from the particle travel times and the transverse displacements of the particles. The transverse displacements are determined from the starting locations to the point where each particle crosses a “collection” plane downstream of the source location oriented normal to the mean flow direction. For this model setup and a flux-weighted, instantaneous injection of particles, Kreft and Zuber (1978) present an analytical solution for the longitudinal dispersivity: 2 ) ( ) ( 2 . .. . . .. . = x m x x t t L s a (Eq. 5) where mt(x) and st(x) are the mean and standard deviation of the travel times at travel distance of x. Other researchers (Desbarats and Srivastava, 1991; Moreno and Tsang, 1994) have found that the longitudinal dispersivity calculated using Equation 5 can fluctuate considerably from one realization to the next. These fluctuations have been attributed to outliers in the travel time distribution. Even with a flux weighted particle source term, there are some particles that travel through very low permeability material and take an excessively long time to reach the collection plane. These outliers give the travel time distribution an extremely long tail and lead to spurious values of the dispersivity. The problem of long time outliers in the calculation of the longitudinal dispersivity is handled here by using the log transform of the travel times. The mean and standard deviation of the log time distribution are determined and then converted back to raw space. Statistics based on the log values of the travel times are much less sensitive to outliers than are those calculated with the raw travel times. Inherent in this log transform filtering of the distribution is the assumption (first paragraph of Section 5.6) that the travel time distribution can be approximated as a lognormal distribution. Other researchers (Khaleel, 1994; Wen and Gomez-Hernandez, 1998) have succesfully applied this log transform approach to the determination of longitudinal dispersivity and we follow those examples in this work. The expressions for the mean and standard deviation of the travel time distribution are now determined as: Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 41 of 60 ] 1 )) ( )[exp( ( ) ( ] 2 / ) ( ) ( exp[ ) ( 2 ln 2 2 2 ln ln - = + = x x m x x x m x m t t t t t t s s s where mln t(x) and s2 ln t(x) are the mean and variance of the natural log travel times at a travel distance x. The transverse dispersivities are calculated as: x x x x x x V H 2 ) ( ) ( 2 ) ( ) ( 33 22 . a . a = = where .22(x) and .33(x) are the variances of the differences between the transverse particle travel position at travel distance x and the initial traverse coordinate at the source location for the transverse X (horizontal) and Z (vertical) positions respectively. These variances are used to approximate the transverse second-order spatial moments, X22 and X33, presented in equations 2b and 2c. This approximation has been used successfully to evaluate transverse dispersion in previous studies (Follin and Thunvik, 1994; Wen and Gomez-Hernandez, 1998). This section has presented an approach for the derivation of longitudinal and transverse dispersivities from the results of particle tracking simulations. This approach is based on the analytical solutions of Dagan (1988) and the assumption (second paragraph of Section 5.6) that a Gaussian distribution in each of the three principal directions can describe the solute plume. The growth of the second moment of the plume in the transverse directions is difficult to derive from the particle tracking results, especially in a strongly heterogeneous medium. Therefore, the variances of the transverse displacements of the particles are used as a proxy for the transverse second moments of the plume. The permeability fields used here differ from those specified in the development of the analytical solutions in that the natural log variance is higher than 1.0 at 9.46 and that the fields are not ergodic nor isotropic. Differences between analytical results and numerical results are compared at the end of this report. (Eq. 6a) (Eq. 6b) (Eq. 7a) (Eq. 7b) Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 42 of 60 Figure 13. Longitudinal (upper image) and transverse (lower image) dispersivity as a function of travel distance as calculated using the analytical solutions in Equations 2a and 2b along with Equations 3 and 4. The two lines correspond to the variances assigned to the theoretical limit of applicability (ln variance = 1.0) and the Yucca Mountain Project data (ln variance = 9.46). Travel Distance (meters) 0 200 400 600 800 1000 Transverse Dispersivity (meters) 0 1 2 3 4 5 6 7 ln Variance 1.0 ln Variance 9.46 Travel Distance (meters) 0 200 400 600 800 1000 Longitudinal Dispersivity (meters) 0 20 40 60 80 100 120 140 160 180 200 ln Variance = 1.0 ln Variance = 9.46 Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 43 of 60 6.2.2 Numerical Implementation Longitudinal and transverse dispersivities are calculated using Equations 5 (along with 6a and 6b) and 7a and 7b above. Input to these equations are the travel distances, the mean and standard deviation of the travel time distribution and the variances of the differences between transverse starting and arrival positions. These travel time and location statistics are determined by modeling streamlines through the heterogeneous fields. The modeling process is outlined as: • A heterogeneous permeability field is selected and converted from the GSLIB format to the FEHM format using the code gs2fehm. The permeability field is the output of either the sisim or combine routine depending on whether it is a background or an enhanced permeability model respectively. • A steady-state flow solution is obtained using FEHM (version 2.00) and the flxo macro within FEHM is activated to write out internodal fluxes within a 100 x 5 gridblock plane of the flow model. • The output of the flxo macro (het_flow.out) is used as input to the code cr8sptr or cr8sptr_cell. These two codes generate the FEHM sptr macro for the flux weighted source and the single grid block source respectively. This code reads in the internodal fluxes in flxo.out and determines the starting locations for the particles using a flux-weighted scheme. The flux weighting determines the number of particles to go into each cell. Within each cell, the particles are given random starting locations. The output of cr8sptr or cr8sptr_cell is the sptr macro for FEHM saved in the sptr.dat file. • A second FEHM run is done for the particle tracks. The sptr macro contained in sptr.dat is used for the particle starting locations. A single time step within FEHM is used and the solution from the previous FEHM run, stored in the het_flow.fin file, is copied to the initial conditions file, het_flow.ini, for this second run. • The FEHM output file het_flow.sptr2 contains the times and locations of all the particles as they move through the three-dimensional domain. This file is used as input to the code calc_D_old that calculates the longitudinal and transverse dispersivities. There are four outputs from calc_D_old: disper_L.out, disper_X.out, disper_Z.out and debug.out. The ”disper” files contain the longitudinal (L) and transverse (X and Z) dispersivities as calcuated in calc_D_old using Equations 5 and 7a and 7b above. The dispersivities are calculated for 9 different travel distances and each column in the output files of calc_D_old is for a different travel distance (increasing distance from leftmost to rightmost column). These five steps as listed in the bullets above are repeated for each realization of the permeability field. For the calculations presented below, 50 realizations are completed for both the background permeability and the enhanced models. The 50 runs are completed by using the UNIX K-shell scripts run_50_back.ksh and run_50_feat.ksh. Output files written by calc_D_old are appended for each realization such that at the end of the 50 runs, each dispersivity output has 50 lines of 9 columns each. Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 44 of 60 For the final analyses, the output files from calc_D_old are transferred to a PC via ftp and each file is read into Excel. Within Excel, the mean and median of each column are calculated. These values are the mean and median value of the dispersivity calculated across all 50 realizations for each of the nine travel distances. Graphs of the median dispersivity as a function of travel distance are shown in the Results section below. The number of particles required for accurate dispersivity calculations is not known a priori. The number used here was 4000 particles. This number was chosen by considering a number small enough to run in a reasonable time and also create output files of a manageable size (approximately 100MB), yet large enough to provide accurate results. The accuracy of the results were checked by selecting one permeability realization and running the dispersion calculations with the same flow solution, but with increasing numbers of particles for each calculation. For the realization examined, the calculated dispersivity became relatively stable when more than 2000 particles were used. Some variation in the number of particles necessary will occur from one realization to another, so 4000 particles were used to ensure that stable values of dispersivity would be calculated on all realizations. 6.2.3 Results Longitudinal, aL, transverse horizontal, aH, and transverse vertical dispersivities, aV, were determined on both the background and the enhanced fracture permeability models. For each fracture permeability model a total of 50 realizations were created and the particle tracking was completed using two different source terms: “point” and “plane”. The point source term begins all particles within the same numerical gridblock. The chosen gridblock is the one with the highest groundwater flux within the central region of the flow model at the same Y coordinate (20.0 meters downgradient of the upper fixed head boundary) as the plane source. The plane source, as described previously, assigns the particles in a flux weighted manner to the gridblocks within the central 400 x 10 m of the flow model at a location 20 meters downstream of the upper fixed head boundary. For either source term, the particles within any single gridblock are distributed randomly. Table 5. Travel distances for the dispersivity calculations. The relative distances are calculated as: (absolute travel distance ) /., where . is set to the model input value of 20.0 meters. Collection Plane Absolute Travel Distance (meters) Relative Travel Distance (Number of correlation lengths) 1 20.0 1.0 2 40.0 2.0 3 80.0 4.0 4 120.0 6.0 5 180.0 9.0 6 240.0 12.0 7 300.0 15.0 8 380.0 19.0 9 460.0 23.0 Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 45 of 60 For each directional dispersivity, the results are shown as an X-Y plot of dispersivity in meters as a function of absolute travel distance. The dispersivity shown in each plot is the median value determined across the 50 realizations. Travel distance is defined as the distance from the source to the collection plane. A total of 9 different collection planes were used at different distances downstream of the source. These distances are shown as both absolute and relative values in Table 5. 6.2.3.1 Longitudinal Dispersivity Longitudinal dispersivity is calculated from the coefficient of variation of the travel time distribution as measured at each collection plane (Equation 5). The numerical results for the two different conceptual models and the two different source terms are shown in Figure 14. Figure 14. Longitudinal dispersivity as a function of travel distance for both the background and enhanced permeability models with both point and plane particle sources. The most notable aspect of Figure 14 is the extreme values of dispersivity that occur with the enhanced permeability model and the plane source. These values are approximately an order of magnitude higher than the values for the background model with either source geometry or for the same enhanced model with the point source. Also, these values are considerably larger than the travel distance. Examination of the results from individual realizations shows that the calculated dispersivities vary drastically from one realization to the next. The distribution of dispersivity values across all 50 realizations exhibits a slightly multimodal behavior. Travel Distance (meters) 0 100 200 300 400 500 Longitudinal Dispersivity (meters) 0 500 1000 1500 2000 2500 3000 3500 Background, Point Source Enhanced, Point Source Background, Plane Source Enhanced, Plane Source Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 46 of 60 These extreme results are interpreted as being due to some realizations producing a strongly bimodal distribution of traveltimes while others do not. Two examples of high permeability features are shown below (Figure 15) as support for this interpretation. The two feature realizations in Figure 15 produce very similar longitudinal dispersivity values at 460 meter travel distance for the point source: realization 10 (upper image) produces a dispersivity of 233 meters while realization 7 (lower image) produces a dispersivity of 163 meters. However, with the planar source geometry, realization 10 (upper image) produces a longitudinal dispersivity of 392 meters while realization 7 (lower image) produces a value of 5,890 meters. In light of the images in Figure 15 and visual examination of other realizations, those realizations with high permeability features that can capture a significant fraction of the particles along the flowpath will have larger dispersivities than those that do not. The high permeability feature geometry that is best suited for capturing the flowpaths is a feature that has an oblique orientation to the principal flow direction. Features that are completely, or close to, parallel with the principal flow direction will not cause such a bimodal distribution of travel times. The lower image of Figure 15 shows a prime example of a high permeability feature that will capture all particles starting in between the two large features, but not capture particles starting outside of these two features. For the point source geometry, the high permeability features create channels of high flux and all of the particles are assigned starting locations within that channel. Dispersivities for the point source geometry are considerably smaller than those calculated for the planar source geometries. The results in Figure 14 also show decidedly non-Gaussian behavior in that dispersivity is highest at short travel distances and then asymptotically decreases to a final value at large travel distances. Compare this behavior to the analytical solution graphed in Figure 13 that shows a gradual rise in longitudinal dispersivity to an asymptotic value at large travel distance. Rather than spreading as a plume in which the longitudinal concentration profile could be described by a Gaussian distribution, as the analytical solution assumes, the travel times are much more variable at shorter travel distances. This high dispersivity at short distances is due to the fracture permeability models creating large spatial variability in groundwater flux across a plane normal to the mean flow direction. In flow through the fracture permeability models exhibiting high variability in log permeability, flux is partitioned into discrete zones much more so than in the lower variance Gaussian based permeability models inherent in the analytical solutions. This partitioning of flux into discrete zones creates large variability in groundwater travel times at short distances resulting in large calculated values of dispersivity. As the travel distance increases, the flowpaths get channeled into relatively higher permeability flowzones and the travel times become more homogeneous. The longitudinal dispersivity results obtained in this numerical study can be compared to results obtained from modeling tracer experiments at the C-wells complex. Three different sets of results have been obtained from tracer tests done between the same boreholes UE-25C #2 and UE-25C #3 (see Fahy, 1997; Geldon, et al. 1997). The straight-line distance between the boreholes is roughly 30 meters and dispersivity values of 1.9 to 2.6 meters (Geldon et al. 1997) 2.4 to 2.6 meters (Fahy, 1997) have been interpreted from matches to the observed tracer data. Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 47 of 60 These tracer-test derived values indicate longitudinal dispersivities on the order of 3.7 to 8.7 percent of the travel distance. In terms of absolute values and percentage of travel distance, the longitudinal dispersivities derived from the tracer tests are all less than the values derived from this numerical study. A reason for this difference may be that transport between the two C wells is thought to occur in only 2 or 3 discrete fracture pathways, while flow in the numerical studies presented here occurs throughout the fracture network. Another difference between the tracer test results and the result in this study is that the flow field for the tracer tests is radially convergent flow, while the numerical models use a uniform flow field. It is possible that the different flow fields also account for the differences in dispersivity values between the studies; however, analysis of these effects is beyond the scope of this study. Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 48 of 60 Figure 15. High permeability features in two example enhanced permeability model realizations that produce drastically different longitudinal dispersivity values. Particle tracking on realization 10 (upper image) produces a longitudinal dispersivity of 392 meters. Particle tracking on realization 7 (lower image) produces a longitudinal dispersivity of 5890 meters. The arrow indicates flow direction. 50 25 0 Elevation (m) 0 100 200 300 400 500 Easting (m) 500 400 300 200 100 0 Northing (m) 50 25 0 Elevation (m) 0 100 200 300 400 500 Easting (m) 500 400 300 200 100 0 Northing (m) Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 49 of 60 Figure 16. Expanded view of Figure 14 showing longitudinal dispersivity as a function of travel distance for the background permeability models with both point and plane particle sources, and the enhanced permeability model with a point source. Figure 16 shows the lower three plots in Figure 14 using an expanded vertical scale. From Figure 16, it is possible to see that the asymptotic values of longitudinal dispersivity for these combinations of permeability model and source geometry all fall within the 110 to 180 meter range. These dispersivities are roughly 24 to 39 % of the travel distance. From the results in Figure 16, it is unclear whether or not the longitudinal dispersivities have truly reached an asymptotic macrodispersivity. The results for the two different source geometries in the background permeability models appear to have reached an asymptote; however, the dispersivities calculated with the enhanced permeability model appear still be rising after 460 meters of travel. Results of the longitudinal dispersivity calculations also show that once the large-scale features are added into the permeability models, these features completely alter the behavior of the particle transport through the domain. Figure 17 shows cross-plots of the longitudinal dispersivity calculated with the enhanced permeability model as a function of the same dispersivity calculated with the background permeability model. For both the planar- and pointsource geometry, there is no correlation between the results obtained on the background permeability model and those obtained with the enhanced model. These results indicate that the nature of the transport process changes completely with the addition of large-scale highpermeability features. Travel Distance (meters) 0 100 200 300 400 500 Longitudinal Dispersivity (meters) 0 100 200 300 400 500 600 700 Background, Point Source Enhanced, Point Source Background, Plane Source Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 50 of 60 Figure 17. Scatterplots showing the longitudinal dispersivity value calculated at a 460 meter travel distance with the enhanced permeability models as a function of the value calculated with the background model. The upper image shows the point source results and the lower image shows the planar source results. Background Model Longitudinal Dispersivity (meters) 0 200 400 600 800 1000 1200 Enhanced Model Longitudinal Dispersivity (meters) 0 200 400 600 800 Background Model Longitudinal Dispersivity (meters) 0 50 100 150 200 250 300 350 400 450 500 Enhanced Model Longitudinal Dispersivity (meters) 0 1000 2000 3000 4000 5000 6000 7000 Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 51 of 60 6.2.3.2 Horizontal Transverse Dispersivity Horizontal transverse dispersivity is calculated from the variance of the horizontal displacements of the particles as measured at each collection plane (Equation 7a). The results for the two different conceptual models and the two different source terms are shown in Figure 18. Figure 18. Horizontal transverse dispersivity as a function of travel distance for both the background and enhanced permeability models with both point and planar particle sources. Similar to the results for longitudinal dispersivity shown in Figures 14 and 16, the results for horizontal transverse dispersivity in Figure 18 show that the largest values of dispersivity occur within the enhanced fracture permeability model when planar source geometry is assumed. Similar to the longitudinal dispersivity results, the planar source allows the particle to sample almost all possible flowpaths. When this source is coupled with an enhanced permeability model, the possibility of preferential channels developing with significantly different transverse flow directions is increased. Recall, that the large scale features have variable orientations about the main flow direction. As an example, consider the field shown in the lower image of Figure 15. While translation in the directions normal to the average flow direction does not cause transverse dispersion by itself, the translations created by the converging features shown in the lower image of Figure 15 will certainly enhance transverse variability in flow paths. The theoretical behavior for transverse dispersivity is that of a rapid rise to a peak value within a few correlation lengths of travel distance and then a gradual decrease to a dispersivity value of Travel Distance (meters) 0 100 200 300 400 500 Horizontal Dispersivity (meters) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 Background, Point Source Enhanced, Point Source Background, Plane Source Enhanced, Plane Source Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 52 of 60 zero. The planar source geometry in the background permeability model shown in Figure 18 reproduces the macroscopic dispersivity behavior predicted by theory (Equation 7a, compare to results in the lower image of Figure 13). It is difficult to determine whether or not results from the planar source geometry in the background permeability model have reached an asymptotic value or are still gradually decreasing (Figure 18). The final values of horizontal transverse dispersivity range from approximately 0.6 to 3.0 meters, or approximately 0.1 to 0.6 % of the travel distance. 6.2.3.3 Vertical Transverse Dispersivity Vertical transverse dispersivity is calculated from the variance of the vertical displacements of the particles as measured at each collection plane (Equation 7b). The results for the two different conceptual models and the two different source terms are shown in Figure 19. Figure 19. Vertical transverse dispersivity as a function of travel distance for both the background and enhanced permeability models with both point and planar particle sources. The vertical dispersivity values shown in Figure 19 demonstrate results that are different from the longitudinal and horizontal values shown above in Figures 14 and 18. For the vertical dispersivity calculations, the background permeability model, not the enhanced model, produces the largest values of dispersivity. Additionally, all combinations of conceptual model and source geometry converge to the same final dispersivity value of approximately 0.2 meters at a travel distance of 460 meters. The final value of transverse vertical dispersivity at the 460 meter travel distance (approximately 0.2 meters) is an order of magnitude smaller than that predicted by the analytical solution Travel Distance (meters) 0 100 200 300 400 500 Vertical Dispersivity (meters) 0.0 0.5 1.0 1.5 Background, Point Source Enhanced, Point Source Background, Plane Source Enhanced, Plane Source Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 53 of 60 (Equations 2c and 4). This result is due to the size of the domain used in the numerical models. The vertical dimension of the flow domain is 50 meters or about 2 to 3 vertical correlation lengths (see Figure 12). The analytical solutions assume an infinite unbounded domain. In practical terms, this means a domain of at least 10 to 20 correlation lengths at a minimum. While this relative size constraint is met in the longitudinal and horizontal transverse directions, it is certainly not met in the vertical transverse direction and thus the calculated dispersivity values are considerably smaller than those predicted by the analytical solution. 6.3 COMPARISON TO ANALYTICAL SOLUTIONS Two important conclusions are noted from a comparison of the analytical solutions to the numerical results presented here: 1) with the exception of one of the longitudinal dispersivity cases, the numerical results show dispersivity values that are similar to results derived from the analytical solutions. 2) while the dispersivity values at 460 meter travel distance are similar between analytical and numerical solutions, the values at shorter distances may be considerably different. Tables 6, 7 and 8 show the values of dispersivity calculated both numerically and analytically for the 460-meter travel distance. Note the lower values of vertical transverse dispersivity in the numerical results relative to the analytical results as discussed above. Table 6. Calculated longitudinal dispersivity values at 460 meter travel distance. For the numerical results, the median value across 50 realizations is given. Longitudinal Dispersivity Calculation Dispersivity after 460 meter travel distance (meters) Analytical (ln s2 = 9.46) 169 Background, Point source 115 Enhanced Point Source 140 Background, Planar source 175 Enhanced Planar Source 2380 Table 7. Calculated horizontal transverse dispersivity values at 460 meter travel distance. For the numerical results, the median value across 50 realizations is given. Horizontal Transverse Dispersivity Calculation Dispersivity after 460 meter travel distance (meters) Analytical (ln s 2 = 9.46) 2.4 Background, Point source 0.6 Enhanced Point Source 0.6 Background, Planar source 0.9 Enhanced Planar Source 3.0 Table 8. Calculated vertical transverse dispersivity values at 460 meter travel distance. For the numerical results, the median value across 50 realizations is given. Vertical Transverse Dispersivity Calculation Dispersivity after 460 meter travel distance (meters) Analytical (ln s 2 = 9.46) 2.4 Background, Point source 0.2 Enhanced Point Source 0.2 Background, Planar source 0.2 Enhanced Planar Source 0.2 Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 54 of 60 7 CONCLUSIONS Fracture data available from Yucca Mountain are used as input to create multiple effectivecontinuum realizations of two different conceptual models of fracture permeability. Both conceptual models are consistent with the available observations. The fracture permeability realizations honor the input cumulative distribution function and also reproduce the specified models of spatial correlation within the constraints imposed on the modeling by the finite size domain. Additionally, the spacing between high permeability features matches the spacing of flowing intervals measured at Yucca Mountain using borehole flow meters. These realizations are used as input to a groundwater flow model and the results of the flow model are used as the basis for streamline particle tracking. The temporal moments of the particle travel times at nine different collection planes downstream of the source area are used to determine the longitudinal dispersivity. The second moment of the transverse particle displacements is used to determine the horizontal and vertical transverse dispersivity values. These results are compared to results obtained with an analytical solution based on stochastic theory of groundwater transport. Results of this dispersivity study show: • Flow in highly heterogeneous fractured media is not dispersive in the classical sense of stochastic theory; rather the fracture permeability pattern leads to channelized flow. Solute plumes will not exist as elliptical bodies, but instead the solute will move as a “ribbon” along preferential pathways of relatively higher permeability. • In spite of the above conclusion, the analytical expressions for longitudinal and transverse dispersivity developed by Dagan (1988) do an adequate job of predicting the macrodispersivity after 460 meters, or roughly 23 correlation lengths, of travel distance. This result appears to be a case of the wrong physics giving the correct answers. With the exception of the vertical transverse dispersivity calculations, the numerical results suggest extremely non-Gaussian behavior at early times (short travel distances) which is contrary to the basis of the analytical solutions. The notable exception to this conclusion of correct results from the analytical solutions is the case of the longitudinal dispersivity calculated with the enhanced permeability model and a planar (distributed) source. • Large-scale, high-permeability features can dominate the solute transport behavior of the saturated zone. In the presence of continuous high-permeability features, using either a point or a distributed planar flux weighted source term, the background fracture permeability pattern has no relation to the transport results. • In general, the dispersivity values have not reached a final asymptotic value after 460 meters of travel distance (23 correlation lengths). The final values for the enhanced permeability model with the planar source appear to be the furthest from an asymptote after 460 meters for both the longitudinal and transverse dispersivity. The background model with the planar source gives the highest values of vertical transverse dispersivity. Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 55 of 60 7.1 RECOMMENDATIONS TO SITE CHARACTERIZATION Based on the conclusions presented above, several recommendations are made (to the site characterization team: Natural Environment Program Operations) regarding further characterization of the fractured rock portion of the saturated zone: • High permeability features control solute transport within fractured rocks. The locations, orientations, lengths and connectivity of these features needs to be defined. It may be possible to map some of these features deterministically using geophysical or drilling techniques; however many of these features will be too hard to detect. The majority of the additional information will most likely be statistical in nature. That information can be very useful in defining conceptual models upon which performance assessment (PA) can base abstractions. • The information on the frequency of flowing features obtained from borehole flowmeters is extremely useful in characterizing the fractured saturated zone. Flowmeter logging should be conducted in as many wells as possible. 7.2 RECOMMENDATIONS TO PERFORMANCE ASSESSMENT Based on the conclusions presented above, several recommendations are made to the performance assessment team regarding transport in the saturated zone: • Due to the channeling nature of the flow system, one-dimensional transport modeling along streamlines (the current approach to PA transport modeling) is strongly recommended. • The conceptual model of large streamtubes of perfectly mixed solute concentrations existing for considerable length along the saturated zone transport pathway, such as was used in the TSPA-VA, is unrealistic. However, it may still be possible to abstract the saturated zone transport model using these types of streamtubes. • The estimates of transverse dispersivities provided by the expert elicitation panel (CRWMS M&O, 1998b) are in general agreement with the results reported in this work. The results presented here do not extend to PA length scales, but all indications suggest that minimal transverse dispersion will occur along the transport pathway. In the vertical dimension, the water table will act upon a shallow plume as a no flow boundary and constrain the possible amount of vertical transverse dispersion. 7.3 UNCERTAINTIES/RESTRICTIONS There are several additional pieces of knowledge that, if obtained, could further the results presented here. Without this additional knowledge, there are some uncertainties in the results as discussed below: • Critical examination of the role of multiple fracture sets on the amount of transverse dispersivity. In this work, with the exception of the features in the enhanced model, all Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 56 of 60 fractures were oriented parallel to flow. Larger amounts of horizontal transverse dispersion may occur if there are multiple fracture sets oriented obliquely to the principal flow directions. Implementing this conceptual model in a continuum numerical model may prove difficult. • Doing a similar study using a discrete fracture network model could check the validity of the continuum model assumption employed herein. • The fractured rocks represent only a fraction of the saturated zone transport pathway. A dispersion study similar to the one presented here should be performed on the alluvial section of the saturated zone. The recent data coming out of the Nye County wells could be used in such a study. 8 REFERENCES 8.1 DOCUMENTS CITED: Altman, S.J.; Arnold, B.W.; Barnard, R.W.; Barr, G.E.; Ho, C.K.; McKenna, S.A.; and Eaton, R.R. 1996. Flow Calculations for Yucca Mountain Groundwater Travel Time (GWTT-95). SAND96-0819. Albuquerque, New Mexico: Sandia National Laboratories. ACC: MOL.19961209.0152. CRWMS M&O 1998a. "Saturated Zone Flow and Transport." Chapter 8 of Total System Performance Assessment-Viability Assessment (TSPA-VA) Analyses Technical Basis Document. B00000000-01717-4301-00008 REV 01. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19981008.0008. CRWMS M&O 1998b. Saturated Zone Flow and Transport Expert Elicitation Project. Deliverable Number SL5X4AM3. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19980825.0008. CRWMS M&O 1999a. M&O Site Investigations. Activity Evaluation, January 23, 1999. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19990317.0330. CRWMS M&O 1999b. M&O Site Investigations. Activity Evaluation, September 28, 1999. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19990928.0224. CRWMS M&O 1999c. Probability Distribution for Flow Interval Spacing. ANL-NBS-MD- 000003. Las Vegas, Nevada: CRWMS M&O. Submit to RPC URN-0183 Dagan, G. 1988. "Time-Dependent Macrodispersion for Solute Transport in Anisotropic Heterogeneous Aquifers." Water Resources Research, 24, (9), 1491-1500. Washington, D.C.: American Geophysical Union. TIC: 245369. Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 57 of 60 Desbarats, A.J. and Srivastava, R.M. 1991. "Geostatistical Characterization of Groundwater Flow Parameters in a Simulated Aquifer." Water Resources Research, 27, (5), 687-698. Washington, D.C.: American Geophysical Union. TIC: 245375. Deutsch, C.V. and Journel, A.G. 1998. GSLIB Geostatistical Software Library and User's Guide, 2nd Edition. New York, New York: Oxford University Press. TIC: 240101. DOE (U.S. Department of Energy) 1998. Quality Assurance Requirements and Description. DOE/RW-0333P, Rev. 8. Washington, D.C.: U.S. Department of Energy, Office of Civilian Radioactive Waste Management. ACC: MOL.19980601.0022. 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. Fahy, M.F. 1997. "Dual-Porosity Analysis of Conservative Tracer Testing in Saturated Volcanic Rocks at Yucca Mountain in Nye County, Nevada." International Journal of Rock Mechanics and Mining Sciences, 34, (3-4), 486. Amsterdam, The Netherlands: Elsevier Science Ltd. TIC: 237601. Follin, S. and Thunvik, R. 1994. "On the Use of Continuum Approximations for Regional Modeling of Groundwater Flow Through Crystalline Rocks." Advances in Water Resources, 17, 133-145. New York, New York: Elsevier Science. TIC: 245366. Geldon, A.L. 1996. Results and Interpretation of Preliminary Aquifer Tests in Boreholes UE- 25c #1, UE-25c #2, and UE-25c #3, Yucca Mountain, Nye County, Nevada. Water-Resources Investigations Report 94-4177. Denver, Colorado: U.S. Geological Survey. ACC: MOL.19980724.0389. Geldon, A.L.; Umari, A.M.A.; Fahy, M.F.; Earle, J.D.; Gemmell, J.M.; and Darnell, J. 1997. Results of Hydraulic and Conservative Tracer Tests in Miocene Tuffaceous Rocks at the C-Hole Complex, 1995 to 1997, Yucca Mountain, Nye County, Nevada. Milestone SP23PM3. Denver, Colorado: U.S. Geological Survey. ACC: MOL.19980122.0412. Gelhar, L.W. and Axness, C.L. 1983. "Three-Dimensional Stochastic Analysis of Macrodispersion in Aquifers." Water Resources Research, 19, (1), 161-180. Washington, D.C.: American Geophysical Union. TIC: 222815. Hsieh, P.A. 1998. "Scale Effects in Fluid Flow Through Fractured Geologic Media." Chapter 12 of Scale Dependence and Scale Invariance in Hydrology. Sposito, G., ed. 335-353. New York, New York: Cambridge University Press. TIC: 245367. Khaleel, R. 1994. "Scale and Directional Dependence of Macrodispersivities in Colonnade Networks." Water Resources Research, 30, (12), 3337-3355. Washington, D.C.: American Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 58 of 60 Geophysical Union. TIC: 245376. Kreft, A. and Zuber, A. 1978. "On the Physical Meaning of the Dispersion Equation and Its Solutions for Different Initial and Boundary Conditions." Chemical Engineering Science, 33, (11), 1471-1480. New York, New York: Pergamon Press. TIC: 245365. Luckey, R.R.; Tucci, P.; Faunt, C.C.; Ervin, E.M.; Steinkampf, W.C.; D'Agnese, F.A.; and Patterson, G.L. 1996. Status of Understanding of the Saturated-Zone Ground-Water Flow System at Yucca Mountain, Nevada, as of 1995. Water-Resources Investigations Report 96- 4077. Denver, Colorado: U.S. Geological Survey. ACC: MOL.19970513.0209. McKenna, S. 1999. Analysis and Modeling Report (AMR) Work Direction and Planning Document for Geostatistical and Heterogeneity Analysis for the SZ Flow and Transport Model, Saturated Zone Flow and Transport PMR. 14012031M3. Las Vegas, NV: CRWMS M&O. ACC: MOL.19990707.0298. Moreno, L. and Tsang, C.F. 1994. "Flow Channeling in Strongly Heterogeneous Porous Media: A Numerical Study ." Water Resources Research, 30, (5), 1421-1430. Washington, D.C.: American Geophysical Union. TIC: 222363. Schmitt, M. and Beucher, H. 1997. "On the Inference of the Boolean Model." Geostatistics Wollongong 1996, International Geostatistics Congress, Wollongong, Australia, September 1996. Pages 200-210. Dordrecht, The Netherlands: Kluwer Academic. TIC: 245400. Sonnenthal, E.L.; Ahlers, C.F.; and Bodvarsson, G.S. 1997. "Fracture and Fault Properties for the UZ Site-Scale Flow Model." Chapter 7 of The Site-Scale Unsaturated Zone Model of Yucca Mountain, Nevada, for the Viability Assessment. Bodvarsson, G.S.; Bandurraga, T.M.; and Wu, Y.S., eds. LBNL-40376. Berkeley, California: Lawrence Berkeley National Laboratory. ACC: MOL.19971014.0232. Sweetkind, D.S. and Williams-Stroud, S.C. 1996. Characteristics of Fractures at Yucca Mountain, Nevada: Synthesis Report. Administrative Report. Denver, Colorado: U.S. Geological Survey. ACC: MOL.19961213.0181. Tsang, C-F. 1997. Drift Scale Heterogeneous Permeability Field Conditioned to Field Data. Milestone SP331BM4. Berkeley, California: Lawrence Berkeley National Laboratory. ACC: MOL.19971125.0915. Tsang, Y.W.; Tsang, C.F.; and Hale, F.V. 1996. "Tracer Transport in an Stochastic Continuum Model of Fractured Media." Water Resources Research, 32, (10), 3077-3092. Washington, D.C.: American Geophysical Union. TIC: 245377. Wen, X-H. and Gomez-Hernandez, J.J. 1998. "Numerical Modeling of Macrodispersion in Heterogeneous Media: A Comparison of Multi-Guassian and Non-Multi-Guassian Models." Journal of Contaminant Hydrology, 30, 129-156. Amsterdam, The Netherlands: Elsevier Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 59 of 60 Science. TIC: 245368. YMP (Yucca Mountain Project) 1997. Q-List. YMP/90-55Q, Rev. 4. Las Vegas, Nevada: Yucca Mountain Site Characterization Office. ACC: MOL.19970714.0474. Zhan, H. 1999. "On the Ergodicity Hypothesis in Heterogeneous Formations." Mathematical Geology, 31, (1), 113-134. New York, New York: Plenum Publishing. TIC: 245383. 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. Zyvoloski, G.A.; Robinson, B.A.; Dash, Z.V.; and Trease, L.L. 1997. User's Manual for the FEHM Application - A Finite-Element Heat- and Mass-Transfer Code. LA-13306-M. Los Alamos, New Mexico: Los Alamos National Laboratory. TIC: 235999. 8.2 CODES, STANDARDS, REGULATIONS AND PROCEDURES AP-3.10Q, Rev. 2 , ICN 0 (C). Analyses and Models. Washington, D.C.: U.S. Department of Energy, Office of Civilian Radioactive Waste Management. ACC: MOL.20000217.0246 QAP-2-0, Rev. 5. Conduct of Activities. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19980826.0209 9 ATTACHMENTS The majority of the attachments are source listings for routines and shells. With the exception of Attachment 1, all attachments are included as electronic files. In most cases, the source file has the same name as the routine or shell with a *.doc extension. In some cases, the routine is made up of separate source and include files. For these routines, all files are listed as “Attachment #@”, where the # refers to the attachment number and the @ is the sub-attachment identifier. ATTACHMENT 1 OCRWM Document Input Reference Sheet ATTACHMENT 2 Listing of the run_50.ksh UNIX shell ATTACHMENT 3 Listing of the cr8_50_feat.ksh UNIX shell ATTACHMENT 4 Listing of the sisim FORTRAN Routine Attachment 4a acorni.doc Attachment 4b beyond.doc Attachment 4c chknam.doc Attachment 4d cova_3.doc Attachment 4e getindex.doc Attachment 4f ksol.doc Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media ANL-NBS-HS-000022 Rev 00 07/14/2000 60 of 60 Attachment 4g locate.doc Attachment 4h ordrel.doc Attachment 4i picksupr.doc Attachment 4j powint.doc Attachment 4k setrot.doc Attachment 4l setsupr.doc Attachment 4m sisim.doc Attachment 4n sortem.doc Attachment 4o sqdist.doc Attachment 4p srchsupr.doc ATTACHMENT 5 Listing of the ellip2_MC FORTRAN Routine Attachment 5a acorni.doc Attachment 5b chknam.doc Attachment 5c ellip2_MC.doc Attachment 5d expon.doc Attachment 5e gcum.doc Attachment 5f setrot.doc Attachment 5g sqdist.doc ATTACHMENT 6 Listing of the combine.c Routine ATTACHMENT 7 Listing of the calc_cdf.c Routine ATTACHMENT 8 Listing of the gs2fehm.c Routine ATTACHMENT 9 Listing of the gamfil FORTRAN Routine Attachment 9a CHKNAM.doc Attachment 9b GAMFIL1.doc Attachment 9c RAND.doc Attachment 9d gamfil.inc.doc ATTACHMENT 10 Listing of the modgeom.c Routine ATTACHMENT 11 Listing of the run_50_feat.ksh UNIX shell ATTACHMENT 12 Listing of the run_50_back.ksh UNIX shell ATTACHMENT 13 Listing of the cr8sptr.c Routine ATTACHMENT 14 Listing of the cr8sptr_cell.c Routine ATTACHMENT 15 Listing of the calc_D_old.c Routine AP-3.15Q.1 Rev. 06/30/1999 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev.: Change: Title: ANL-NBS-HS-000022 00 Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media (S0015) Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 2a 1 64 FR (Federal Register) 8640. Disposal of High-Level Radioactive Waste in a Proposed Radiologic Repository at Yucca Mountain. Proposed rule 10 CFR 63. Readily available. Referen ce Only Rules regulating the disposal of high-level waste in the United States N/A N/A N/A N/A 2 Altman, S.J.; Arnold, B.W.; Barnard, R.W.; Barr, G.E.; Ho, C.K.; McKenna, S.A.; and Eaton, R.R. 1996. Flow Calculations for Yucca Mountain Groundwater Travel Time (GWTT-95). SAND96-0819. Albuquerque, New Mexico: Sandia National Laboratories. ACC: MOL.19961209.0152. Pages 47-65 Referen ce Only GWTT-95 report cited as an example of stochastic continuum modeling N/A N/A N/A N/A 3 CRWMS M&O. 1998a. Total System Performance Assessment – Viability Assessment (TSPA-VA) Analyses Technical Basis Document. Chapter 8 Saturated Zone Flow and Transport. B00000000-01717- 4301-00008 Rev 01. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19981008.0008. DTN: SNT05082597001.003 Figure 8- 27 Referen ce Only Section 4, Inputs and Section 6, Analysis Citing distribution of permeabilities measured in the saturated zone N/A N/A N/A N/A 4 CRWMS M&O. 1998b. Saturated Zone Flow and Transport Expert Elicitation Project. SL5X4AM3. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19980825.008. Referen ce Only Results of expert elicitation used as a basis for dispersion values in TSPA-VA document N/A N/A N/A N/A 5 CRWMS M&O 1999a. Activity Evaluation of M&O Site Investigations. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19990317.0330. Referen ce Only Yucca Mountain Site Characterization Procedure N/A N/A N/A N/A AP-3.15Q.1 Rev. 06/30/1999 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev.: Change: Title: ANL-NBS-HS-000022 00 Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media (S0015) Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 6 CRWMS M&O 1999b. Activity Evaluation of M&O Site Investigations. Las Vegas, Nevada: CRWMS M&O. ACC: MOL. 19990928.0224 Referen ce Only Yucca Mountain Site Characterization Procedure N/A N/A N/A N/A 7 CRWMS M&O 1999c. Probability Distribution for Flowing Interval Spacing. ANL-NBS-MD-000003. Las Vegas, Nevada: CRWMS M&O ACC: Pending, DTN: SN990&T0571599.001. Figure 11 Referen ce Only Section 4, Inputs and Section 6, Analysis Citing distribution of spacings between flowing features N/A N/A N/A N/A 8 Dagan, G. 1988. “Time-Dependent Macrodispersion for Solute Transport in Anistropic Heterogeneous Aquifers.” Water Resources Research, 24 (9), 1491-1500. Washington, D.C.: American Geophysical Union. TIC: 245369. Referen ce Only Journal article presenting stochastic-analytical solution for dispersion in heterogeneous aquifers N/A N/A N/A N/A 9 Desbarats, A.J. and Srivastava, R.M. 1991. “Geostatistical Characterization of Groundwater Flow Parameters in a Simulated Aquifer.” Water Resources Research, 27 (5), 687-698. Washington, D.C.: American Geophysical Union. TIC: 245375. Referen ce Only Numerical study of dispersion in two-dimensional heterogeneous aquifer N/A N/A N/A N/A 10 Deutsch, C.V. and Journel, A.G. 1998. GSLIB: Geostatistical Software Library and User’s Guide. 2nd Edition. New York, New York: Oxford University Press. TIC: 224174. Pages 175-180, 182-183 Referen ce Only User's manual for some of the software routines used in this analysis N/A N/A N/A N/A AP-3.15Q.1 Rev. 06/30/1999 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev.: Change: Title: ANL-NBS-HS-000022 00 Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media (S0015) Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 11 DOE (U.S. Department of Energy) 1998. Quality Assurance and Requirements Description. DOE/RW-0333P, REV 8. Washington D.C.: U.S. Department of Energy, Office of Civilian Radioactive Waste Management. ACC: MOL.19980601.0022. Referen ce Only OCRWM Quality Assurance and Description N/A N/A N/A N/A 12 Dyer, J. R. 1999. "Revised Interim Guidance Pending Issuance of New U. S. Nuclear Regulatory Commission (NRC) Regulations (Revision 01, July 22, 1999), for Yucca Mountain, Nevada." Letter from J. R. Dyer (DOE) to D. R. Wilkins (CRWMS M&O), September 9, 1999, OL&RC:SB-1714, with enclosure, "Interim Guidance Pending Issuance of New U. S. Nuclear Regulatory Commission (NRC) Regulations (Revision 01)" . ACC: MOL.19990910.0079. Referen ce Only Interim guidance for work on Yucca Mountain Site Characterization Project N/A N/A N/A N/A 13 Follin, S. and Thunvik, R. 1994. “On the Use of Continuum Approximations for Regional Modeling of Groundwater Flow Through Crystalline Rocks.” Advances in Water Resources, 17 (3), 133-145. Southampton, U.K.: C.M.L. Publications; New York, New York: Springer-Verlag. TIC: 245366. Referen ce Only Example of stochastic continuum modeling in heterogeneous aquifers and example of numerical calculation of transverse dispersivity N/A N/A N/A N/A 14 Geldon, A.L. 1996. Results and Interpretation of Preliminary Aquifer Tests in Boreholes, UE-25c #1, UE-25c #2, and UE- 25c #3, Yucca Mountain, Nye County, Nevada. Water-Resources Investigations Report 94-4177. Denver, Colorado: U.S. Geological Survey. ACC: MOL.19980724.0389. Referen ce Only Results and interpretation of aquifer testing completed at Choles. N/A N/A N/A N/A AP-3.15Q.1 Rev. 06/30/1999 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev.: Change: Title: ANL-NBS-HS-000022 00 Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media (S0015) Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 15 Geldon, A.L.; Umari, A.M.A.; Fahy, M.F.; Earle, J.D.; Gemmell, J.M.; and Darnell, J. 1997. Results of Hydraulic and Conservative Tracer Tests in Miocene Tuffaceous Rocks at the C-Hole Complex, 1995 to 1997, Yucca Mountain, Nye County, Nevada. SP23PM3. Las Vegas, Nevada: U.S. Geological Survey, Yucca Mountain Project Branch. ACC: MOL.19980122.0412. Table 9, p. “Hydrauli c Tests- 37” Referen ce Only Results and interpretation of aquifer and tracer testing completed at C-holes. N/A N/A N/A N/A 16 Gelhar, L.W. and Axness, C.L., 1983, Three- Dimensional Stochastic Analysis of Macrodispersion in Aquifers, Water Resources Research, 19(1), pp. 161-180, Washington, D.C.: American Geophysical Union. TIC 222815. Referen ce Only Stochastic analytical solution for dispersion in heterogeneous aquifers. N/A N/A N/A N/A 17 Hsieh, P.A. 1998. “Scale Effects in Fluid Flow Through Fractured Geologic Media.” Scale Dependence and Scale Invariance in Hydrology. Pages 335-353. Editor: Sposito, G. Cambridge, United Kingdom; New York, New York: Cambridge University Press. TIC: 245367. Referen ce Only Review article on fractured rock hydrology. N/A N/A N/A N/A 18 Khaleel, R. 1994. “Scale and Directional Dependence of Macrodispersivities in Colonnade Networks.” Water Resources Research, 30 (12), 3337-3355. Washington, D.C.: American Geophysical Union. TIC: 245376. Referen ce Only Numerical study on dispersion in fractured rocks. Example of the log transform as a means of decreasing the significance of outliers in travel time distribution. N/A N/A N/A N/A AP-3.15Q.1 Rev. 06/30/1999 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev.: Change: Title: ANL-NBS-HS-000022 00 Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media (S0015) Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 19 Kreft, A. and Zuber, A. 1978. “On the Physical Meaning of the Dispersion Equation and Its Solutions for Different Initial and Boundary Conditions.” Chemical Engineering Science, 33 (11), 1471-1480. Oxford, United Kingdom; New York, New York: Pergamon Press. TIC: 245365. Referen ce Only Analytical solution for longitudinal dispersion in aquifer. N/A N/A N/A N/A 20 Luckey, R.R.; Tucci, P.; Faunt, C.C.; Ervin, E.M.; Steinkampf, W.C.; D’Agnese, F.A.; and PattersonG.L. 1996. Status of Understanding of the Saturated-Zone Groundwater Flow System at Yucca Mountain, Nevada as of 1995. Water Resources Investigations Report 96-4077. Denver, Colorado: U.S. Geological Survey. TIC: 227084. Referen ce Only Summary status report on saturated zone hydrology at Yucca Mountain. N/A N/A N/A N/A 21 Moreno, L. and Tsang, C-F. 1994. “Flow Channeling in Strongly Heterogeneous Porous Media: A Numerical Study.” Water Resources Research, 30 (5), 1421-1430. Washington, D.C.: American Geophysical Union. TIC: 222363. Referen ce Only Numerical study of dispersion in highly heterogeneous rock. N/A N/A N/A N/A 22 Schmitt, M. and Beucher, H. 1997. “On the Inference of the Boolean Model.” Geostatistics Wollongong ’96, International Geostatistics Congress, Wollongong, Australia, September 1996, 200-210. Baafi, E.Y. and Schofield, N.A., eds. Boston, Massachusetts: Kluwer Academic. TIC: 245400. Referen ce Only Conference proceeding article on the inference of Boolean models. N/A N/A N/A N/A AP-3.15Q.1 Rev. 06/30/1999 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev.: Change: Title: ANL-NBS-HS-000022 00 Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media (S0015) Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 23 Sonnenthal, E.L.; Ahlers, C.F.; and Bodvarsson, G.S. 1997. “Chapter 7: Fracture and Fault Properties for the UZ Site- Scale Flow Model.” The Site-Scale Unsaturated Zone Model of Yucca Mountain, Nevada, for the Viability Assessment. LBNL- 40376, 7-1 : 7-34. Bodvarsson, G.S.; Bandurraga, T.M.; and Wu, Y.S., eds. Berkeley, California: Lawrence Berkeley National Laboratory. ACC: MOL.19971014.0232. Referen ce Only Summary of fracture information completed for the Site-Scale unsaturated zone model. N/A N/A N/A N/A 24 Sweetkind, D.S. and Williams-Stroud, S.C. 1996. Characteristics of Fractures at Yucca Mountain, Nevada: Synthesis Report. Denver, Colorado: U.S. Geological Survey. TIC: 239254. DTN:GS960808314222.003 Page 52 and Figures 21 and 22 Referen ce Only Section 4, Inputs and Section 6, Analysis Citing distribution of mapped fracture lengths. N/A N/A N/A N/A 25 Sweetkind, D.S. and Williams-Stroud, S.C. 1996. Characteristics of Fractures at Yucca Mountain, Nevada: Synthesis Report. Denver, Colorado: U.S. Geological Survey. TIC: 239254. Figure 17 Referen ce Only Section 4, Inputs and Section 6, Analysis Citing distribution of measured fracture orientations. N/A N/A N/A N/A 26 Tsang, C.F. 1997. Drift-Scale Heterogeneous Permeability Field Conditioned to Field Data. Milestone SP331BM4. Berkeley, California: Lawrence Berkeley National Laboratory. ACC: MOL.19971125.0911. Corrected by ACC: MOL.19990622.0022. Referen ce Only Initial report on modeling of fracture permeability based on air testing in unsaturated zone. N/A N/A N/A N/A AP-3.15Q.1 Rev. 06/30/1999 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev.: Change: Title: ANL-NBS-HS-000022 00 Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media (S0015) Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 27 Tsang, Y.W.; Tsang, C.F.; Hale, F.V.; and Dverstorp, B. 1996. “Tracer Transport in a Stochastic Continuum Model of Fractured Media.” Water Resources Research, 32 (10), 3077-3092. Washington, D.C.: American Geophysical Union. TIC: 245377. Referen ce Only Journal article giving an example of a stochastic continuum model used for tracer transport. N/A N/A N/A N/A 28 Wen, X-H. and Gomez-Hernandez, J.J. 1998. “Numerical Modeling of Macrodispersion in Heterogeneous Media: A Comparison of Multi-Gaussian and Non-Multi- Gaussian Models.” Journal of Contaminant Hydrology, 30 (1-2), 129-156. Amsterdam, The Netherlands; New York, New York: Elsevier Science B.V. TIC: 245368. Referen ce Only Journal article describing the numerical modeling of solute transport in two-dimensional heterogeneous media. Example of numerical technique used to calculate transverse dispersivity. N/A N/A N/A N/A 30 YMP (Yucca Mountain Project) 1997. Q-List. YMP/90-5Q. Las Vegas, Nevada: Yucca Mountain Site Characterization Office. ACC: MOL.19970714.0474. Referen ce Only Yucca Mountain Site Characterization Project listing of Q data and documents N/A N/A N/A N/A 31 Zhan, H. 1999. “On the Ergodicity Hypothesis in Heterogeneous Formations.” Mathematical Geology, 31 (1), 113-134. New York, New York: Plenum. TIC: 245383. Referen ce Only Journal article defining the limits of the ergodic hypothesis N/A N/A N/A N/A 32 Zyvoloski, G.A., Robinson, B.A., Dash, Z.V. and Trease, L.L., 1997a, User’s Manual for the FEHM Application-A Finite- Element Heat and Mass-Transfer Code, LANL Report LA-13306-M. Los Alamos, New Mexico, Los Alamos National Laboratory. TIC: 2359999. Referen ce Only User's manual for the FEHM code N/A N/A N/A N/A AP-3.15Q.1 Rev. 06/30/1999 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev.: Change: Title: ANL-NBS-HS-000022 00 Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media (S0015) Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 33 Zyvoloski, G.A., Robinson, B.A., Dash, Z.V. and Trease, L.L., 1997b, Summary of the Models and Methods for the FEHM Application-A Finite-Element Heat and Mass-Transfer Code, LANL Report LA- 13307-MS. Los Alamos, New Mexico, Los Alamos National Laboratory. TIC 235587. Referen ce Only Background information for the FEHM code. N/A N/A N/A N/A 32 Software Code: run_50.ksh. STN: PENDING. Entire TBV- 3537 UNIX shell to check feature spacing of back ground permeability models. TBV indicates pending documentation of verification of routine. 1.0 X 33 Software Code: cr8_50_feat.ksh. STN: PENDING. Entire TBV- 3538 UNIX shell to create enhanced perm models. TBV indicates pending documentation of verification of routine. 1.0 X 34 Software Code: sisim. STN: PENDING. Entire TBV- 3539 FORTRAN routine to create indicator geostatical simulation of permeability. TBV indicates pending documentation of verification of routine. 1.0 X 35 Software Code: ellip2_MC. STN: PENDING. Entire TBV- 3540 FORTRAN Routine to create high permeability features for enhances model. TBV indicates pending documentation of verification of routine. 1.0 X 36 Software Code: combine.c. STN: PENDING. Entire TBV- 3541 Preprocessing C routine to combine the background perm models and high permeability features into the final enhanced permeability model. TBV indicates pending documentation of verification of routine. 1.0 X AP-3.15Q.1 Rev. 06/30/1999 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev.: Change: Title: ANL-NBS-HS-000022 00 Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media (S0015) Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 37 Software Code: calc_cdf.c. STN: PENDING. Entire TBV- 3542 Postprocessing C Routine to check the cdf of the permeability model at each indicator threshold. TBV indicates pending documentation of verification of routine. 1.0 X 38 Software Code: gs2fehm.c. STN: PENDING. Entire TBV- 3543 Preprocessing C routine to reformat the permeability models from GSLIB to FEHM formate routine. TBV indicates pending documentation of verification of routine. 1.0 X 39 Software Code: gamfil. STN: PENDING. Entire TBV- 3544 Postprocessing FORTRAN Routine to check the feature spacing in the permeability models. TBV indicates pending documentation of verification of routine. 1.0 X 40 Software Code: modgeom.c. STN: PENDING. Entire TBV- 3545 Preprocessing routine to change the random seed in the ellips_MC input files. TBV indicates pending documentation of verification of routine. 1.0 X 41 Software Code: run_50_feat.ksh. STN: PENDING. Entire TBV- 3546 UNIX shell to run the flow, particle tracking and dispersion calculations 50 times on the enhanced permeability models. TBV indicates pending documentation of verification of routine. 1.0 X 42 Software Code: run_50_back.ksh. STN: PENDING. Entire TBV- 3547 UNIX shell to run the flow, particle tracking and dispersion calculations 50 times on the background permeability models. TBV indicates pending documentation of verification of routine. 1.0 X AP-3.15Q.1 Rev. 06/30/1999 OFFICE OF CIVILIAN RADIOACTIVE WASTE MANAGEMENT DOCUMENT INPUT REFERENCE SHEET 1. Document Identifier No./Rev.: Change: Title: ANL-NBS-HS-000022 00 Modeling Sub Gridblock Scale Dispersion in Three-Dimensional Heterogeneous Fractured Media (S0015) Input Document 8. TBV Due To 2. Technical Product Input Source Title and Identifier(s) with Version 3. Section 4. Input Status 5. Section Used in 6. Input Description 7. TBV/TBD Priority Unqual. From Uncontrolled Source Unconfirmed 43 Software Code: cr8sptr.c. STN: PENDING. Entire TBV- 3548 Preprocessing C Routine to create flux weighted input file for FEHM sptr macro. TBV indicates pending documentation of verification of routine. 1.0 X 44 Software Code: cr8sptr_cell.c. STN: PENDING. Entire TBV- 3549 Preprocessing C Routine to create single cell input file for FEHM sptr macro. TBV indicates pending documentation of verification of routine. 1.0 X 45 Software Code: calc_D_old.c. STN: PENDING. Entire TBV- 3550 C Routine to calculate longitudinal and transverse dispersivities. TBV indicates pending documentation of verification of routine. 1.0 X