Main Document

Return to the Main Document

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" (CRWMS, 1999c). 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 "(CRWMS, 1999c) 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

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_MC routine. Two different UNIX shells are used to run the codes that build the permeability models and then accomplish post processing. One shell is used for the background permeability models, run_50.ksh (Attachment II) and one shell for the enhanced permeability models, cr8_50_feat.ksh (Attachment III). These shells are discussed below in section 6.1.3.2.

3.2.1 Geostatistical Simulations

The routine sisim (version 2.0, Attachment IV) 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). A Pentium personal computer was used to run the sisim code.

The routine modgeom (version 1.0, Attachment X) is used to change the random number seed after each call to the sisim or ellip2_MC 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 documentation for this subroutine is listed in Attachment V. 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. The combine subroutine is discussed in Section 6.1.2.2, and the documentation for this subroutine is listed in Attachment VI. combine was run on a SUN UltraSparc workstation.

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 VII and IX 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 documentation for this subroutines is listed in Attachment VIII.

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_feat.ksh and run_50_back.ksh are listed in Attachment XI and XII 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 contained in the output DTN (SN0005T0505500.001). 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 contained in the output DTN 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.

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 documentation for this subroutine is listed as Attachment XIII.

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 documentation for this subroutine is listed as Attachment XIV.

The directory path for the particle tracking runs in the output DTN 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 XV. 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 in the output DTN 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 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.

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.

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, 2000). 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 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.

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 constructed here, it is assumed that the longer length fractures dominate the flow system and especially the relatively faster flow paths within the fractured rocks. This assumption is based on observations of rapid pressure responses in outlying wells during the C-hole 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 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 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 m². 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 2000) 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 2000), 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).

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 2000).

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.

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).

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 2000, 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 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 m², and the maximum is –10.0 m². The minimum log10 permeability of the multi-well hydraulic tests, –12.2 m², 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 m², 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 m², 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 m² 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 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 m² 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 m² 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.

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 95^th 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.

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 fine-scale 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 in the output DTN submitted to the TDMS for this report. The directory containing the input and output files for the background permeability models is: /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 III. 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/ in the output DTN submitted to the TDMS)
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 m²) 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 IX.
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 X.

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/ of the output DTN submitted to the TDMS. The source code for this shell is in Attachment II.

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:

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.

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.
The variograms are calculated across each realization and are compared to the input variogram model. Unlike the permeability distribution and the feature spacing 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.

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). [Figure 4, Figure 5, Figure 6, Figure 7] Along each transect, the simulated permeability values were compared against the assumed lower limit of permeability for the flowing features (-12.2 log10 m²). 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 m² and threshold 6 at –12.2 log10 m². 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.

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 m² 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, X₁₁, X₂₂, X₃₃, of the particle trajectories with respect to time (after Wen and Gomez-Hernandez 1998):

		(Eq. 1a)
		(Eq. 1b)
		(Eq. 1c)

Where t is time, U is the mean particle velocity, _L is the longitudinal dispersivity, _H is the horizontal transverse dispersivity and _V 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):

		(Eq. 2a)
		(Eq. 2b)
		(Eq. 2c)

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:

		(Eq. 3)
		(Eq. 4)

where

_Y² 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:

(Eq. 5)

where m_t(x) and _t(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 log-normal distribution. Other researchers (Khaleel 1994; Wen and Gomez-Hernandez 1998) have successfully 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:

		(Eq. 6a)
		(Eq. 6b)

where m_{ln t}(x) and ²_{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:

		(Eq. 7a)
		(Eq. 7b)

where

₂₂(x) and ₃₃(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, X₂₂ and X₃₃, 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.

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.

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, _L, transverse horizontal, _H, and transverse vertical dispersivities, _V, 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.

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.

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.

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. 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.

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 point- source 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 high-permeability features.

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.

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 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.

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 (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.

7. Conclusions

Fracture data available from Yucca Mountain are used as input to create multiple effective-continuum 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.

Copies of the input and output electronic files used in this analysis have been submitted to the TDMS as Data Tracking Number (DTN): SN0005T0505500.001.

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.

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 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. (Q) Activity Evaluation, September 28, 1999. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19990928.0224.

CRWMS M&O 1999c. Geostatistical and Heterogeneity Analysis for the SZ Flow and Transport Model. Analysis and Modeling Report (AMR) Work Direction and Planning Document.. Las Vegas, NV: CRWMS M&O. ACC: MOL.19990707.0298.

CRWMS M&O 2000. Probability Distribution for Flow Interval Spacing. ANL-NBS-MD-000003 Rev. 00. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.20000602.0052

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.

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.