Rock Properties Model (RPM3.1) Rev 00, ICN 01 MDL-NBS-GS-000004 January, 2000 1. PURPOSE The purpose of this Analysis and Model Report (AMR) is to document Rock Properties Model (RPM) 3.1 with regard to input data, model methods, assumptions, uncertainties and limitations of model results, and qualification status of the model. The report also documents the differences between the current and previous versions and validation of the model. The rock properties models are intended principally for use as input to numerical physicalprocess modeling, such as of ground-water flow and/or radionuclide transport. The constraints, caveats, and limitations associated with this model are discussed in the appropriate text sections that follow. This work was conducted in accordance with the following planning documents: WA-0344, “3-D Rock Properties Modeling for FY 1998” (SNL 1997), WA-0358, “3-D Rock Properties Modeling for FY 1999” (SNL 1999), and RPM3.1, (CRWMS M&O 1999c). These work plans describe the scope, objectives, tasks, methodology, and implementing procedures for model construction. The constraints, caveats, and limitations associated with this model are discussed in the appropriate text sections that follow. The work scope for this activity consists of the following: 1. Conversion of the input data (laboratory measured porosity data, x-ray diffraction mineralogy, petrophysical calculations of bound water, and petrophysical calculations of porosity) for each borehole into stratigraphic coordinates 2. Re-sampling and merging of data sets 3. Development of geostatistical simulations of porosity 4. Generation of derivative property models via linear coregionalization with porosity 5. Post-processing of the simulated models to impart desired secondary geologic attributes and to create summary and uncertainty models 6. Conversion of the models into real-world coordinates. The conversion to real world coordinates is performed as part of the integration of the RPM into the Integrated Site Model (ISM) 3.1; this activity is not part of the current analysis. The ISM provides a consistent volumetric portrayal of the rock layers, rock properties, and mineralogy of the Yucca Mountain site and consists of three components: • Geologic Framework Model • RPM, which is the subject of this AMR • Mineralogic Model. The interrelationship of the three components of the ISM and their interface with downstream uses are illustrated in Figure 1. Figure 2 shows the geographic boundaries of the RPM and other component models of the ISM. Figure 1. Interrelationships Between Component Models, Integrated Site Model, and Downstream Uses Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 20 of 150 Hydrologic Flow and Radionuclide Transport Models Through Saturated Zone (SZ) Rock Properties Model (RPM) Total System Performance Assessment Geologic Framework Model (GFM) Mineralogic Model (MM) Repository Design Hydrologic Flow and Radionuclide Transport Models Through Unsaturated Zone (UZ) INTEGRATED SITE MODEL Figure 2. Index map showing the location of the Geologic Framework Model, the Rock Properties Model (this report), and the Mineralogic Model. Contour Interval 200 Feet Boundary of Rock Properties Model Boundary of Geologic Framework Model and Mineralogic Model Exploratory Studies Facility Cross-Block Drift Yucca Crest Tonsil Ridge Drill Hole NEVADA TEST SITE STATE OF NEVADA YUCCA MOUNTAIN SITE Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 21 of 150 Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 22 of 150 INTENTIONALLY LEFT BLANK Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 23 of 150 2. QUALITY ASSURANCE The modeling effort was evaluated in accordance with QAP-2-0, “Conduct of Activities,” and was determined to be quality affecting and subject to the requirements of the Quality Assurance Requirements and Description (DOE 1998). This evaluation is documented in Activity Evaluation of M&O Site Investigations (Q) (CRWMS M&O 1999a, 1999b). Accordingly, all efforts to construct the rock properties models have been conducted in accordance with approved quality assurance procedures. All work was performed in accordance with the Quality Assurance Requirements and Description under the auspices of the Civilian Radioactive Waste Management System Management and Operating contractor’s quality assurance program. Modeling work was performed in accordance with QAP-SIII-3, “Scientific Notebooks.” Prior to the effective date of QAP-SIII-3, AP-3.10Q, and AP-SIII.2Q, modeling work was performed in accordance with Sandia National Laboratories QAIP 20-2, “Scientific Notebooks.” The procedures that currently pertain to the rock properties modeling effort are listed in CRWMS M&O (1999c). This Analysis and Model Report was developed in accordance with AP-3.10Q, “Analyses and Models.” Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 24 of 150 INTENTIONALLY LEFT BLANK Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 25 of 150 3. COMPUTER SOFTWARE AND MODEL USAGE 3.1 SOFTWARE TRACKED BY CONFIGURATION MANAGEMENT The Rock Properties Model was constructed principally using geostatistical algorithms that are part of the public-domain GSLIB geostatistical subroutine library (Deutsch and Journel 1992, 1998). The major modeling codes subject to configuration management are listed below in Table 1, together with a brief description of their functionality. These computer software packages were obtained from configuration management and are judged appropriate for use in this type of modeling activity. The software was used within the range of calibration (to the extent that this statement applies to these codes). All codes listed in Table 1 run on Intel-type personal computers under the Microsoft Windows NT Version 4.0 operating system, except EARTHVISION, which is a Unix-based software product operating on a Silicon Graphics Octane workstation under IRIX Version 6.4 operating system. A number of additional “software routines” have been developed as part of the rock properties modeling process. These routines are listed in Table 2 and they are documented for quality assurance purposes in the scientific notebooks associated with this activity (Table 11); the use of scientific notebooks is explained in Section 6.3.2. Exceptions to this practice are included as attachments to this report (noted in tables that follow). With the exception of STRATC4 (which is a macro-type software routine), all software routines listed in Table 2 run as Microsoft Fortran Table 1. Software Tracked by Configuration Management Code Name Version STN Number Brief Description GSLIB V1.4 SGSIM 1.4 10110- 1.4MSGSI MV1.40-00 Generates conditional or unconditional gaussian simulations of a continuous variable; optional normal-score forward and back transformation (from the software library, GSLIB) GSLIB V2.0 IK3D 2.0 10122- 2.0MIK3D V2.0-00 Conducts indicator kriging of continuous or discrete variables using either hard or soft data or both (from the software library, GSLIB) GSLIB V1.4 NSCORE 1.201 10109- 1.4MNSC OREV1.20 1-01 Transforms a distribution of values to standard-normal form while preserving quantile relationships (from the software library, GSLIB) EARTHVISION 4.0 30035 V4.0 3-D earth science modeling package; used only to produce visualizations (report figures; e.g., Figure 38) of the Integrated Site Model 3.0/Rock Properties Model 3.0 for this report MATCHUP 12/5/98 10196- 12598-00 Software routine that creates a single file of depth-matched values from two columnar input files MATCH 12/5/98 10195- 12598-00 Modification of software routine MATCHUP for one GSLIB format file Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 26 of 150 PowerStation Version 4.0 and have been compiled and executed on Intel-type personal computers under the Microsoft Windows NT Version 4.0 operating system. Table 2. Software Routines Routine Name Version/ Date1 Reference Brief Description Acquired Software TRANS 1.3 Attachment III Transforms an input distribution to match the histogram of a target distribution (both in GSLIB format) while preserving the quantile relationships (from the software library, GSLIB) BICALIB 2.0 Rautman and McKenna 1998, p. 1143 Calibrates soft data (GSLIB format) for use as prior cdf values from matched hard-soft data pairs (from the software library, GSLIB) BACKTR 1.2 Rautman and McKenna 1998, p. 1127 Transforms a standard-normal distribution (GSLIB format) to match a reference histogram (from the software library, GSLIB; inverse of program NSCORE [Table 1]) Developed Software BUD 6/25/98 Rautman and McKenna 1998, p. 1013 Reads a “log-analysis standard” .las file and extracts specified well log traces TWOFOOT 6/25/98 Rautman and McKenna 1998, p. 1020 Reads an ascii file of well log data and resamples it on a user-specified spacing, averaging values over two-foot intervals COREGPC 5/22/99 Rautman 1999 p. 307 Post-processes a porosity simulation (GSLIB format) in standardnormal format and an unconditional simulation with the same spatial correlation structure, and generates a coregionalized simulation of a secondary property given a correlation coefficient between the two variables ETYPE 6/30/98 Rautman and McKenna 1998, p. 1134 Reads a set of simulation output (GSLIB format) files and creates summary E-type model; optional uncertainty output as standard deviation of simulations ZEOLITE5 9/10/98 Rautman and McKenna 1998, p. 1266 Post-processes a set of simulated Ks values (GSLIB format) together with an indicator model of “zeolite” alteration that is used as a template to insert random gaussian Ks values at altered grid notes; appropriate only for altered units VITROPHY RE 9/10/98 Rautman and McKenna 1998, p. 1274 Post-processes a set of simulated Ks values (GSLIB format) together with the corresponding porosity simulation and substitutes uniform random Ks values at grid nodes where the simulated porosity value is less than 0.05; appropriate only for TSw model unit. COORDS 12/17/98 Rautman and McKenna 1998, p. 1287 Adds Nevada state plane (x,y,z) coordinates to a GSLIB-format output file given a grid specification and a corresponding structure contour file of top elevations and thicknesses STRATC4 4 Attachment II Computes unit-specific stratigraphic coordinates using Equation 2 from Section 6.4.3.2 using a SIGMA PLOT transform (macro). Notes: 1. Formal “version numbers” were not used for software routines under SNL QAIP 19-1; instead, the “Date” is shown that the documentation associated with the routine was entered into the referenced scientific notebook. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 27 of 150 3.2 EXEMPT SOFTWARE Various nonqualified software products were also used during the course of this modeling activity to assist in organizing, managing, manipulating, comparing, and displaying data and information and computing standard statistical measures. This includes commercial off-the-shelf software and other industry-standard software. These software packages are used to perform support activities and are not used as the controlled source of information for the analysis, as defined in AP-SI.1Q. They are therefore not required to be qualified under these procedures. EXCEL95 and/or EXCEL97 were used to compile and organize rock properties measurements by drill hole and depth, and later by modeling unit prior to writing the final ASCII text files used as input to the modeling algorithms. Standard cell-based arithmetic operations of Excel were used on a temporary basis to check for and identify duplicate sample data (by depth), but the results of these calculations did not change any measurement values, and the columns involved were deleted prior to writing the final ASCII text input files. Standard cell-based arithmetic operations were also used to construct two new variables in selected data files related to hydrous-phase mineral alteration. One variable represents the “total hydrous mineral phase” content, and was constructed simply as the sum of the four relevant mineral fractions described in Section 6.4.3.1. The other variable is an indicator variable constructed using the logic of Equation 1, also described in Section 6.4.3.1. These computations are immediately verifiable by visual inspection. Cell-based arithmetic was also used to adjust the core bound-water contents for parity with the petrophysically based bound-water contents, as described in Section 6.4.7.1. This computation also is immediately verifiable by visual inspection. SIGMA PLOT, versions 2.0 and 5.0, was used to compile and organize rock properties measurements and to generate graphic displays of drill hole data for “sanity checking” that the selection and organization of the drill hole data was conducted properly. SIGMA PLOT arithmetic spreadsheet operations were also used to create the “bound-water” content variable described in Sections 6.4.3.1 and 6.4.7.1 as the simple difference between two original variables. This calculation is immediately verifiable by visual inspection. SIGMA PLOT was also used to generate various statistical measures, including histograms, cumulative distribution functions, crossplots, and linear regressions of selected data sets. EVS, version 3.75, was used to generate selected graphic visualizations (report figures, including block diagrams and cross-sectional views: e.g., Figures 38 [top half] and 39) of the completed rock properties models in stratigraphic coordinates only. The geostatistical and other modeling capabilities of EVS were not used for this analysis. VARIO, versions 1.16 and 1.20, from the geostatistical package UNCERT, was used to compute standard three-dimensional experimental variograms using the input data files. Proper operation of this software was verified by visual inspection of the results and by a simple manual verification documented in Attachment IV. VARIOFIT, version 1.20, from the geostatistical package, UNCERT, was used to fit standard variogram models to the experimental variograms produced with VARIO. Proper operation of this Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 28 of 150 software was verified by visual inspection of the results and by a simple manual verification documented in Attachment V. HISTPLT, version 1.2, from the geostatistical software library GSLIB (Deutsch and Journel 1992), was used to generate histograms and summary statistics for various data sets, including “sanity checks” of completed simulations and other models. Proper operation of this software was verified by visual inspection of the results and by comparison of results with similar values produced by other software (especially SIGMA PLOT, above). GAM3, version 1.2, from the geostatistical software library GSLIB (Deutsch and Journel 1992), was used to generate exhaustive three-dimensional variogram data for plotting as a “sanity check” of completed simulations and other models. Proper operation of this software was verified by visual inspection of the results. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 29 of 150 4. INPUTS 4.1 DATA AND PARAMETERS There are seven different classes of data used as input to rock properties modeling. Four of these categories (Sections 4.1.1 through 4.1.4) involve actual measurements of rock material properties, either in the laboratory or by calculation from geophysical measurements obtained in the field. A fifth class of property (Section 4.1.5) is also derived from in-situ geophysical measurements; it is sufficiently different in nature and used for a different purpose, therefore it is discussed separately. The remaining two classes of input data (Sections 4.1.6 and 4.1.7) consist of the stratigraphic contacts bounding the four rock properties modeling units. The sixth class consists of actually observed (measured) contacts. The seventh class involves interpreted or modeled contacts, used in cases where a given drill hole did not go deep enough to penetrate a unit completely. These seven types of data are described in the paragraphs that follow. Data tracking numbers (DTNs) associated with each type of data are provided in tables in the relevant section. The qualification (To Be Verified (TBV)) status of these data sources is provided in Attachment I, Document Information Reference System (DIRS). 4.1.1 Laboratory Core Porosity Data Laboratory measurements of porosity (measurement process described by Flint 1998, pp. 11 to 17) have been obtained for core samples taken from boreholes within both the potential repository block and immediately adjoining areas. DTNs associated with the laboratorymeasured core porosity values used in the Rock Properties Model are presented in Table 3. Flint (1998, pp. 11 to 17, Table 3) reports two porosity values for each sample: one obtained after drying the sample in a relative humidity (RH) oven at elevated humidity levels (60°C and 65 percent RH), and the other after drying in a 105°C oven (oven-dried) and ambient, but very low RH. The distinction is important, as Flint (1998, pp. 17 and 32 to 33, Figure 12) has demonstrated that matrix saturated hydraulic conductivity is better correlated with the RH porosity values. Because part of the intent of this rock properties modeling exercise is to model spatial variability in the hydraulic conductivity of the site, the RH core porosity values were used exclusively as input to the geostatistical modeling of porosity. Table 3. Source of Input Laboratory Core Porosity Data Data Source Description Reference RH core porosity data DTN: MO9911INPUTRPM.000 RH porosity, SD-6 DTN: GS980808312242.014 RH porosity, WT-24 DTN: GS980708312242.010 Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 30 of 150 4.1.2 Calculated Petrophysical Porosity Data Borehole petrophysical logs have been converted to quantitative estimates of porosity by CRWMS M&O (1999d, pp. 5-1 to 5-2). These data provide substantial information regarding spatial heterogeneity across the entire site area (Figure 2), particularly for regions remote from the potential repository block for which no core samples exist. Use of data from these regions outside the immediate rock properties model boundaries is important in identifying long-range spatial correlation patterns that reduce uncertainty within the modeled area. The DTN associated with the calculated petrophysical porosity data is presented in Table 4. CRWMS M&O (1999d, p. 5-1) generated two different porosity values using petrophysics. These values are identified in his files as POREF, which corresponds approximately to the RH laboratory core measurements, and POROTOT, corresponding approximately to the 105°C ovendried samples. Using the same logic that influenced the choice of the RH laboratory core measurements, POREF values have been used exclusively as input to porosity modeling for Integrated Site Model 3.0. The data files in Table 4 also include a calculated petrophysical value, VWC, representing the volumetric water content, derived from neutron logging (CRWMS M&O 1999d, p. 5-2). Values of VWC have been substituted for POREF in certain geologic intervals and under certain circumstances that are described in Section 6. 4.1.3 Laboratory Measured Secondary Property Data A second class of laboratory measurements performed on core samples consists of “secondary” material properties that generally are of greater interest in design or performance assessment modeling than simple porosity (completed models are referred to in Section 6 as “derivative” properties). Laboratory values of bulk density, matrix saturated hydraulic conductivity (Ks) (Flint 1998, pp. 11 to 17), thermal conductivity, and additional measurements of porosity (specifically those correlated with the secondary properties) were also measured for selected samples, in addition to the “primary” porosity data of Section 4.1.1. This measurement of multiple properties on the same physical specimen allows identification of statistically valid cross-variable correlations. Such correlations thus permit modeling the spatial variability of these secondary or derivative properties in regions where no actual measurements of these properties exist. DTNs associated with the secondary laboratory property data are presented in Table 5. 4.1.4 X-Ray Diffraction Indicators of Mineral Alteration The identification and modeling of “hydrous-phase mineral alteration” (see Section 6.2.4) as it affects matrix Ks is complicated by a general scarcity of both mineralogic and hydraulic conductivity data across the entire Yucca Mountain site area. As described in more detail Table 4. Source of Input Petrophysical Porosity Data Data Source Description Reference POREF, VWC DTN: MO9910POROCALC.000 Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 31 of 150 (Section 6.2.4), the approach is to use both direct (x-ray diffraction analyses [XRD]) and indirect (petrophysical data) in combination to model this type of alteration. Quantitative XRD mineralogic data (as mineral fractions of several zeolite minerals and smectite clays) are used as input to the alteration modeling portion of the rock properties modeling effort. The mineralogic data are both converted to hard (certain) indicators of mineral alteration and used to calibrate petrophysical data (Section 4.1.5 below) as soft (uncertain) indicators. Generation of both types of indicators is described in Section 6.4.7. The DTN associated with the XRD mineralogic analyses is presented in Table 6. 4.1.5 Petrophysical Indicators of Hydrous-Phase Mineral Alteration The use of both direct and indirect indicators to model hydrous-phase mineral alteration as it affects the magnitude and spatial distribution of matrix Ks was noted in Section 4.1.4 above. Petrophysically derived data, specifically the POROTOT and POREF calculated porosity values described in Section 4.1.2, are used to generate soft indicators of mineral alteration. The DTN associated with the petrophysical porosity data used to create the petrophysical indicators of hydrous-phase mineral alteration is presented in Table 7. The mechanics of generating and calibrating those indicators is presented in Section 6.4.7. 4.1.6 Observed (Measured) Lithostratigraphic Contacts Rock properties modeling has been conducted within stratigraphic coordinates, which represent the relative vertical position of each measured property value within a model unit. This stratigraphic coordinate conversion requires depth values for both the upper and lower contact of each aggregate model unit in each borehole, as described in Section 6.3.2. Typically, the required depth values were physically observed, either in core, in geophysical logs, or in downhole video Table 5. Source of Input Secondary Property Data Data Source Description Reference Porosity, bulk density, Ks data DTN: MO9911INPUTRPM.000 Porosity, thermal conductivity samples DTN: SNL01A050593001.007 Thermal conductivity, NRG- holes DTN: SNL01A05059301.005 Table 6. Source of Input X-ray Diffraction Mineralogy Data Data Source Description Reference Zeolite mineral contents DTN: LA9910DB831321.001 Table 7. Source of Input Petrophysical Data Used for Identification of Mineral Alteration Data Source Description Reference POREF, POROTOT, VWC DTN: MO9910POROCALC.000 Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 32 of 150 records. DTNs associated with the observed lithostratigraphic contact data are presented in Table 8. 4.1.7 Modeled Lithostratigraphic Contacts The vast majority of model-unit contacts required for the stratigraphic coordinate conversion are taken directly from a tabulation of observed contacts (Table 8). However, a small number of contacts were projected using Geologic Framework Model (GFM) 3.0 or GFM3.1 because the drill hole did not fully penetrate the model unit. This situation most frequently involves projecting the depth of the lower contact of a given modeling unit where the borehole bottomed within that unit. A few upper contacts had to be projected where the top of a particular model unit had been eroded at the borehole location. The stratigraphic coordinate conversion is also complicated for boreholes that intercept a fault. In this case, it has been necessary to use modeled isochore information to reconstruct the missing portion of the stratigraphic section and infer the “unfaulted depths” based on unit-thickness trends obtained using the GFM. These instances are documented in the scientific notebooks associated with this modeling activity. DTNs associated with the geologic framework models used to extract the projected lithostratigraphic contacts are presented in Table 9. Note that Rock Properties Model (RPM) 3.1 was created simply by inserting the new drill hole data for holes SD-6 and WT-24 into the data set from RPM3.0 (created using projected contacts from GFM3.0) and regenerating the model. 4.2 CRITERIA This Analysis and Model Report complies with the U.S. Department of Energy 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 geology of the site in support of the License Application (Subpart B, Section 21(c)(1)(ii)), and the definition of geologic parameters and conceptual models used in performance assessment (Subpart E, Section 114(a)). Table 8. Source of Input Observed Lithostratigraphic Contacts Data Source Description Reference Contacts for Geologic Framework Model 3.0 DTN: MO9811MWDGFM03.000 Contacts for SD-6 DTN: SNF40060298001.001 Contacts for WT-24 DTN: SNF40060198001.001 Table 9. Source of Input Projected Lithostratigraphic Contacts and Isochore Information Data Source Description Reference Geologic Framework Model 3.0 DTN: MO9804MWDGFM03.000 Geologic Framework Model 3.1 DTN: MO9901MWDGFM31.000 Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 33 of 150 4.3 CODES AND STANDARDS No codes or standards are applicable to the Rock Properties Model. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 34 of 150 INTENTIONALLY LEFT BLANK Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 35 of 150 5. ASSUMPTIONS The principal assumption involved in constructing the rock properties models involves the use of “porosity-as-a-surrogate” for modeling the spatial variability of other, “secondary” material properties that are typically of greater interest in performance modeling than porosity itself, but which are almost universally undersampled at Yucca Mountain. This concept of using abundant porosity data as a surrogate for other properties is not new. Rautman and McKenna (1997, p. 13) list a sampling of published references using the porosity-as-a-surrogate technique. For Rock Properties Model 3.0, porosity has been used to model the spatial distributions of bulk density, saturated hydraulic conductivity, and thermal conductivity. The concept of porosity as a surrogate is based on empirically observed correlations of porosity with undersampled secondary properties. A consequence of undersampling is that the spatial variability of the undersampled variable cannot be described confidently on a stand-alone basis, let alone such that the joint spatial continuity patterns of the two (or more) variables cannot be reproduced simultaneously. It is important to understand that modeling the spatial distribution of several material properties without considering the inter-variable correlations can lead to highly unrealistic input to physical-process modeling codes, which in turn can lead to highly unreasonable estimates of performance parameters. Simply sampling randomly from separate (univariate) probability density functions may easily produce such un-physical combinations as a low porosity–high thermal conductivity–high hydraulic conductivity tuff. The severity of neglecting cross-variable correlations in modeling spatially variable domains increases as physical-process modeling attempts to capture multiple coupled processes. Using porosity as a surrogate for various other materials in modeling Yucca Mountain is supported by consideration of the physics involved in the site-specific rock units being modeled. For example, for a given rock type, increasing the volume of pore space must decrease the bulk density of the rock mass. The part of the rock that “isn’t there” is available to hold fluids, but it contributes nothing to the total mass contained within a unit volume: the definition of bulk density. Again for a given rock type, the conduction of heat energy through the material is directly related to the density (or, inversely to the pore space) of the material. All else being equal, a higher porosity–lower density tuff will conduct heat less readily, leading to a lower thermal conductivity value. Note here that it is the total amount of void space in a rock that affects thermal conductivity, not simply the amount of pore space that is conducting water within the unsaturated zone. This fact has important implications to modeling of whole-rock thermal conductivity in the presence of large-diameter (to 1 meter) lithophysal cavities. As another example, although hydraulic conductivity is not generally well correlated with porosity across many classes of soils and/or rock materials, the empirical observation at Yucca Mountain is that this correlation is quite strong within limited groupings of rock types. Specifically, both welded and nonwelded lithologies appear to be associated with a continuum of saturated hydraulic conductivity values (see Section 6.4.5.1). Evidently, unless affected by some additional physical process (such as zeolitic or other similar alteration), there is a strong relationship between the progressive, overall reduction in porosity and the progressive reduction in the average diameter of the passages between interconnected pores (which is what exerts principal control on the flow of water through the existing pore space) across this continuum of Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 36 of 150 nonwelded to densely welded materials. Conversely, a genetic process that changes the diameter and/or geometry of the pore throats (e.g., zeolitization), while not commensurately filling in the total quantity of void space, can reduce the hydraulic conductivity of the rock by orders of magnitude while leaving the porosity essentially unchanged. In fact, both these cases are observed and modeled at Yucca Mountain. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 37 of 150 6. ANALYSIS/MODEL 6.1INTR ODUCTION A ND ROCK PROPER TIES MODEL OVER VIEW Yucca Mountain is located in the southwestern Nevada volcanic field and consists of faulted and tilted blocks composed of layered sequences of ash-flow, ash-fall, and bedded tuffs of Miocene age (Sawyer et al. 1994, pp. 1304 to 1318). The Rock Properties Model (RPM) is a 3-D (three-dimensional), discretized numerical representation of spatial variability and heterogeneity of a number of fundamental bulk and hydrologic material properties for the majority of the rocks within the unsaturated zone (UZ) at Yucca Mountain. The model divides the model volume into four internally lithologically similar model units. In descending stratigraphic order these are (1) the Paintbrush nonwelded (PTn) model unit, (2) the Topopah Spring welded (TSw) model unit, (3) the Calico Hills nonwelded (CHn) model unit, and (4) the Prow Pass Tuff (Tcp). The Tiva Canyon Tuff was not modeled for two primary reasons, which are discussed in detail in Rautman and McKenna (1997, p. 10). First, such a model is not needed because the Tiva Canyon Tuff is highly fractured, so that matrix properties are subordinate to the role of fractures in hydrologic flow. Second, because most boreholes are located in ravines eroded into the Tiva, few boreholes penetrate a meaningful stratigraphic section of the Tiva Canyon Tuff and data are accordingly sparse. Table 10 presents a correlation of the RPM model units with other project stratigraphies. For all four model units, the modeled material properties are: • Matrix porosity • Whole-rock bulk density • Matrix saturated hydraulic conductivity. For the TSw model unit, additional material properties are: • Lithophysal porosity • Whole-rock thermal conductivity. “Lithophysal” porosity and the two “whole-rock” properties in the above list are intended to represent the effect of the large (up to 1 meter in diameter) lithophysal cavities that are a very prominent feature of certain intervals within the Topopah Spring Tuff and the potential repository horizon. The rock properties model is tied geometrically to the bounding surfaces of model units within Geologic Framework Model (GFM) 3.1. There are three fundamentally different types of models included within the RPM. First, a suite of 50 simulated property models have been generated for each material property using Monte Carlostyle geostatistical conditional simulation techniques. These models reproduce both the measured data values and the full range of spatial heterogeneity and spatial continuity exhibited by those data. At locations other than those of the measured values, material properties are generated (simulated) randomly from the local conditional probability density functions, and the variability of these values is interpreted to represent the uncertainty associated with predictions Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 38 of 150 Table 10. Correlation Chart for Model Stratigraphy Note: 1. Shaded areas represent “header” lines for subdivided units. Source: DTN: MO9510RIB00002.004 Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 39 of 150 of material properties at these unsampled locations. The second type of model included in the RPM are summary expectation (“E-type”) models for each rock property. These E-type models are computed as the node-by-node mathematical expectation (average) across the 50 simulated realizations, and they represent the property values most likely to be encountered at each discretized location. The third type of model included in the RPM is also a summary-type model, and these are computed as the node-by-node standard deviations of the 50 individual simulated property models. This third type of model is intended to provide users with a first-pass estimate of geologic uncertainty throughout the potential repository site, where “geologic uncertainty” consists of that uncertainty that results from less-than-exhaustive site characterization. The two primary data sources used to model the spatial variability of rock properties at Yucca Mountain are laboratory-measured core samples and down-hole measured petrophysical property values. Porosity is the primary variable for which spatial variability and heterogeneity have been modeled directly from measured values. The other material properties are derivative properties that have been generated from the simulated porosity models using the geostatistical technique known as linear coregionalization. Additional information has been incorporated from x-ray diffraction (XRD) measurements of (principally) zeolite mineral content. This form of rock alteration is considered to be an important control on the variability of matrix saturated hydraulic conductivity. Correlation between XRD alteration mineralogy and down-hole petrophysical measurements has been used to extend modeling of this type of alteration into regions from which no mineralogical samples have been obtained. The rock properties models are intended principally for use as input to numerical physical-process modeling, such as of ground-water flow and/or radionuclide transport. In anticipation of such use on a site scale, the models are discretized uniformly in the horizontal plane on a 200-m by 200-m (656.2-ft by 656.2-ft) grid. Vertical discretization is unit-dependent at a nominal 2-m (7-ft) spacing in the PTn model unit and 5-m (16-ft) spacing elsewhere. The rock property models can be regenerated at virtually any desired level of discretization to meet the needs of downstream model users. 6.2 CHANGES BETWEEN MODEL VERSIONS Four significant changes in modeling methodology and data from those used in RPM2 (Rautman and McKenna 1997; note that RPM2 is the first “numbered” rock properties model) have been implemented for the current version of the rock properties model (RPM3). These changes are described briefly in this section. Expanded discussions are presented as appropriate in the sections on Data and Parameters (Section 4.1) and Development of the Models (Section 6.4). 6.2.1Cha nges in Model Un it Sub divisions and Modeled Region The first major change in the rock properties models of Integrated Site Model (ISM) 3.1 is that the combined Calico Hills-Prow Pass model unit of ISM2 (Rautman and McKenna 1997) has been subdivided into two new modeling units along the formation boundary at the top of the Prow Pass Tuff (Table 10). This revised subdivision scheme was adopted because the area of rock properties modeling for ISM3 is smaller than that for ISM2. The smaller model area is a consequence of increased focus on the immediate repository area by the YMP site-scale UZ hydrologic model. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 40 of 150 This decrease in the modeled volume resulted in an improved spatial distribution of boreholes within the remaining area and decreased the need to aggregate data in search of statistical mass, as that term was described by Rautman and McKenna (1997, pp. 10 to 11). 6.2.2 Adjustment of Model Unit Contacts to Reflect Mineralogy A second change is the adjustment of contacts between the welded interior core of the Topopah Spring Tuff (model unit TSw) and the respectively overlying and underlying PTn and CHn nonwelded units. For RPM3, these contacts were moved outward one lithostratigraphic interval, with the result that the breaks are now between the Tptrv3–Tptrv2 and the Tptpv2–Tptpv1 lithostratigraphic units (Table 10). The rationale for this contact adjustment is that the revised boundaries better separate rocks with the potential to have been altered to zeolite minerals from high-temperature-crystallized rocks that almost never undergo such alteration. The spatial distribution of one derivative property, matrix saturated hydraulic conductivity, is strongly affected by the presence or absence of zeolitic minerals; thus a more accurate representation of this particular transition is the objective. Note that such modeling of alteration at this time has been conducted only for the CHn and Tcp model units. 6.2.3 Revised Petrophysical Data The third and most significant, but perhaps least obvious change between the rock properties models of ISM2 and 3, is that all of the petrophysically based porosity data have been recalculated by CRWMS M&O (1999d, pp. 5-1 to 5-2) following identification of a number of problems with the earlier data sets (see Rautman and McKenna 1997, pp. 138 to 139). Specifically, the original petrophysical porosity data for ISM2 involved two different calculational approaches (described by Nelson 1996, pp. 17 to 23 and CRWMS M&O 1996, pp. 26 to 29), with the result that some spatial variability appeared to represent only the source of the calculated values (see discussion in Rautman and McKenna 1997, pp. 26 to 29 and 49 to 51). A separate issue involved the absolute magnitude of some of the porosity values calculated by Nelson (1996; the data values themselves are contained in GS960708312132.002) within selected stratigraphic intervals, principally the densely welded vitric subzones (Tptrv1, Tptpv3) of the Topopah Spring Tuff. As a consequence, it was unclear if some portion of the inferred spatial heterogeneity was related simply to the mathematical calculations and not to physical geology. For RPM3, CRWMS M&O (1999d, pp. 4-1 to 4-25) re-examined all of the underlying petrophysical logs, and processed the necessary fundamental petrophysical traces through a common data-reduction algorithm to calculate porosity. Specifically, use of the “balanced-water approach” using resistivity logs (CRWMS M&O 1996, p. 28) to calculate water saturation and volumetric water content was abandoned. Instead, porosity was calculated directly from the downhole density logs (separately above and below the water table), with volumetric water content derived directly from calibrated neutron logs (CRWMS M&O 1999d, p. 4-13 to 4-25 and 30 to 31). This recalculation and common processing history eliminate one (non-geologic) source of apparent variability. Recalculation of the petrophysical porosity data involved not only a common data-reduction algorithm, but also calibration of the input bulk density trace against a set of internal standards Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 41 of 150 (CRWMS M&O 1999d, pp. 3-1 to 3-2 and 4-2 to 4-5). Bulk density is the principal petrophysical measurement used in the porosity calculation. The internal standards are provided by four distinctive lithostratigraphic units, at least one of which is present in virtually each of the site boreholes. These units are the crystal-rich densely welded vitric (Tptrv1) and nonlithophysal (Tptrn) lithostratigraphic units, and the crystal-poor middle nonlithophysal (Tptpmn) and densely welded vitric (Tptpv3) units. Based on evidence from recovered drill core that these four intervals are of relatively uniform material, a set of density-correction factors was developed for each borehole (CRWMS M&O 1999d, pp. 6 to 9, Table 3). The density logs throughout each borehole were then adjusted proportionately (by individual logging run if appropriate) so that the average value in these distinctive zones indicated the correct bulk density. 6.2.4 Revised Approach to Identifying Hydrous-Phase Mineral Alteration The term “hydrous-phase mineral alteration” is used in this report to indicate changes in initial (volcanic) mineral compositions that collectively act to reduce the matrix hydraulic conductivity of laboratory specimens (and in-situ rocks) compared with equivalent unaltered rocks of approximately the same porosity (Rautman and McKenna 1997, pp. 35 to 40, Figures 25 and 26 and Table 8; Flint 1998, pp. 31 to 33 and Figure 12). Such rocks typically exhibit mineralogical alteration either to zeolite minerals, principally but not exclusively clinoptilolite and mordenite, or to expandable-layer clay minerals (“smectite” or montmorillonite; Flint 1998, p. 33). Within the two rock properties modeling units of principal interest (CHn and Tcp), hydrous-phase mineral alteration translates fairly precisely to zeolitic alteration. However, this distinction is maintained with respect to petrophysically based indicators of reduced Ks to emphasize that this approach is unable to distinguish zeolites and montmorillonitic clays in the downhole environment. A fourth change in rock properties modeling procedure involves the approach used to identify the spatial distribution of hydrous-phase mineral alteration, which affects the coregionalization of hydraulic conductivity with porosity in the lower two (CHn, Tcp) of the four properties modeling units. For version 2 of the rock properties models, each grid node was identified as either likely or unlikely to exhibit such alteration based on a spatial continuity model identified through petrophysical measurements (see Rautman and McKenna 1997, pp. 43 to 44 and 57 to 59). For version 3, altered grid nodes were modeled as a weighted combination of both hard (XRD analyses; Section 4.1.4) and soft (petrophysical; Section 4.1.5) data. Generation and use of the combined alteration data set is described in see Section 6.4.7. 6.3 PHILOSOPHICAL AND CONCEPTUAL MODELING FRAMEWORK This section is intended to provide a relatively high-level conceptual description of the geostatistical modeling approach, including a brief discussion of how stochastic material-property models such as these fit into the overall performance modeling of the Yucca Mountain site. The methodology has been described in some detail by Rautman and McKenna (1997, pp. 6 to 69). Additional information regarding geostatistical modeling techniques in general may be found by following the references cited in the Rautman and McKenna document. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 42 of 150 6.3.1Modelin g Phil osophy and Conceptual App roach Licensing of the Yucca Mountain site as a geologic disposal site for nuclear waste will require quantitative predictions of waste-isolation performance of both the rocks that form Yucca Mountain and the engineered barrier system for an extended period of time into the future. These predictions will require the use of numerical modeling techniques in an effort to capture critical aspects of highly complex physical processes; for example, the flow of ground-water and the transport of radionuclide contaminants under both unsaturated and saturated conditions. Modeling the performance of the site will be influenced both by present-day conditions and by future conditions that must account for perturbation of ambient conditions by the thermal pulse related to the presence of heat-generating radioactive materials and potential changes in climate. A fundamental principle involved in the numerical representation of real-world physical processes is that the relevant material properties of the modeled domain must be known at all positions within that domain. However, in contrast to this requirement for an “exhaustive” spatial description, the process of describing or characterizing a site invariably consists of collecting observations of properties or state variables at a limited number of locations, the exact positions of which are frequently determined by less-than-optimal external factors. This is particularly true for the 3-D characterization of a geologic site, such as Yucca Mountain. Because descriptive characterization is limited both by access (particularly to the subsurface) and by the availability of resources, that description is necessarily incomplete. Therefore, the exhaustive description of a geologic site for purposes of numerical modeling requires the prior assumption of some type of conceptual model for the site, which is then implemented to assign the necessary properties and other variables at every relevant point in space. Many types of conceptual models of varying complexity are available to guide an exhaustive representation of material properties at Yucca Mountain. The simplest possible conceptual model of rocks at Yucca Mountain is to assume that the material properties needed for numerical process modeling are homogeneous and uniform throughout the site. However, although such a model might possess some utility for rough, back-of-the-envelope calculations, even the most cursory inspection of Yucca Mountain indicates that such a representation is vastly oversimplified and of limited value in addressing regulatory-related issues. 6.3.1.1 Capturing Heterogeneity A slightly refined, but still quite simplified, conceptual model for the site makes use of the fact that, at Yucca Mountain (as in many other geologic environments), the lithologic deposits were produced by relatively widespread but temporally variable geologic processes. In particular, the volcanic activity responsible for the formation of Yucca Mountain was episodic in nature, with thick, widespread ash deposits produced by near-instantaneous (geologically speaking) eruptions separated by thin inter-eruption deposits that probably represent much longer intervals of time. If the time represented by a progressively accumulating geologic deposit is considered to be “frozen” and preserved in the vertical dimension, then the resulting conceptual model is one of successive subhorizontal layers, that may be broken and tilted, or folded and otherwise moved about at some later time. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 43 of 150 The GFM component of the ISM is such a layered representation of Yucca Mountain. However, without further refinement, a layered rock-properties model corresponding precisely to the GFM still relies upon the prior assumption that, within each originally subhorizontal layer, the material properties of interest are essentially uniform and homogeneous. Geologic studies of the volcanogenic rocks at Yucca Mountain and of similar deposits elsewhere in the world indicate that the geologic processes responsible for deposition of these materials vary both temporally and areally. For example, variations in cooling rates caused by local conditions affect the material properties in the resultant rocks. This spatial variation of process has produced spatial heterogeneity of material properties in all three dimensions. Yet despite the existence of vertical and lateral heterogeneity, the spatial distribution of material properties within geologic layers is not simply “random.” Knowledge of property values at one location imposes limits on the values of those properties likely to exist at “nearby” locations Therefore, an alternative conceptual model of “filling-in” a geologic framework with values randomly assigned from some inferred univariate distribution without regard for other nearby values (spatial correlation) is an unnecessary oversimplification (and potentially an unwarranted distortion) of the real world. 6.3.1.2 Heterogeneity vs. Uncertainty In contrast to heterogeneity, which is an objective feature of the real world, uncertainty is a knowledge-based concept. Distinguishing properly between uncertainty (as a state of imperfect information resulting from less-than-complete observation) and spatial heterogeneity (as a state of being, unaffected by the availability or lack of information) becomes critically important in the application of predictive engineering methods to the geologic environment. Natural earth materials are heterogeneous to a far greater degree than most materials of conventional engineering interest. Heterogeneity in rock properties affects the actual operation of physical processes, even though there is only one, unique heterogeneous reality. The impact of that heterogeneity on the numerical approximation (modeling) of those physical processes is only compounded by geologic uncertainty. Incomplete information must be accounted for in predictive modeling, as must the effects of material properties that are different in different locations. A key attribute, therefore, of the rock properties modeling activities, has been the description and quantification of the effects of geologic uncertainty on the physical description of the Yucca Mountain site. 6.3.1.3 Geostatistical Methods The current approach to modeling rock properties as part of the ISM attempts to strike a balance between the inherent simplification of a hollow-shell GFM and the near-infinite degree of complexity that is the real world. Selected major lithostratigraphic horizons are used as the constraining (framework) boundaries for a statistically based description of the measured rock material properties that have been sampled within those boundaries. Geostatistical methods have been used to create the exhaustive material property descriptions comprising the rock properties model. Geostatistical methods in general are one of a variety of methods for distributing isolated measurements of different attributes in space, and thus for modeling spatial heterogeneity. A fundamental principle underlying all geostatistical techniques is the quantification and use of some measure of spatial correlation, which may be defined informally as the degree to which samples “close” to one another resemble each other more than do samples “far” away from each Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 44 of 150 other, where “close” and “far” are defined from the data values. Furthermore, unlike many other methods for predicting the material property attributes of a large volume from direct observation of a relatively minuscule fraction of that volume, geostatistical methods offer a quantitative and more-or-less rigorous approach to the issues of knowledge-based uncertainty discussed in Section 6.3.1.2. Within the purview of geostatistical methods are two broad classes of algorithms for predicting attributes at unsampled locations constrained by some limited set of actual measurements: estimation and simulation. Geostatistical estimation is focused on the prediction of the attribute values most likely to be encountered at a given spatial position, and may be thought of as modeling the expected value of a variable of interest. Geostatistical estimation is most frequently described using the term, kriging, and it is simply a weighted-average interpolation method using some neighborhood of near-by relevant data. A common thread connecting all estimation methodologies (including non-geostatistical ones) is that they are interpolation techniques directed toward producing a model in which the estimated values grade progressively and generally smoothly away from the data locations and away from one another. The other broad class of geostatistical methods comprises a variety of simulation algorithms. In contrast with estimation, geostatistical simulation places principal emphasis on reproducing the input data values and the overall statistical character (including the spatial correlation characteristics) exhibited by the data ensemble, the total collection of input values. Models produced by geostatistical simulation typically do not grade smoothly between measured data values, but rather are more highly variable at the same time that they represent the broad heterogeneity structure of the measurements. These techniques are conceptually equivalent to the Monte Carlo simulation process frequently employed in engineering analyses. In common with other Monte Carlo simulation approaches, the emphasis is less on the specific predicted values, which are in effect simply the products of a random number generator with certain “desirable” properties, and much more on evaluation of the space of uncertainty associated with some performance measure computed to represent the behavior of the modeled system. 6.3.1.4 Evaluating the Consequences of Heterogeneity and Uncertainty A schematic diagram of the geostatistical process for combining statistical description with the Monte Carlo generation of multiple replicate models, which are then evaluated for performance consequences is presented in Figure 3 (see also Rautman and McKenna 1997, Figure 6). The logic is straightforward Monte Carlo, except that the need to capture both heterogeneity and spatial correlation jointly causes the geostatistical simulation process to consist of drawing entire “intact” property models (rather than single property values) from a conceptual distribution of alternative “realities.” As suggested by the upper row of panels in Figure 3, the individual measured values of a material property are combined with the overall statistical character of the complete data set (the data ensemble) and with the spatial correlation patterns exhibited by the data to produce replicate “exhaustive” models (center part of Figure 3) showing the distribution of that material property in space. Each replicate model reproduces the measured data at the data locations, and the overall variability of all the values in the model reproduces the histogram of the measurements. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 45 of 150 Figure 3. Conceptual representation of a Monte Carlo process incorporating geostatistical simulation techniques as the basis for assessing the impact of geologic uncertainty on a performance measure relevant to licensing of a geologic repository. A “consequence analysis” is any postsimulation mechanism for computing a measure of performance across the suite of replicate stochastic simulations (from Rautman and McKenna, 1997, Figure 6). Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 46 of 150 Additionally, the spatial correlation structure of the model, evaluated as a whole, approximates the spatial correlation among the input measurements. Completing the Monte Carlo methodology, each of these individual realizations (property models) may then be processed through some relevant transfer function or consequence analysis (for example, a radionuclide-transport computer code), and the likelihood of various acceptable vs. unacceptable performance responses evaluated. This final step is represented in the lower portion of Figure 3. 6.3.2 Methodology Overview Construction of the 3-D rock properties models has been guided by the philosophically separate but related concepts of heterogeneity and uncertainty described in the preceding section. The models attempt to make maximum use of “deterministic” genetic processes to constrain the modeling of measured property values away from the physical locations of those measurements. The effort to capture genetic processes as they relate to material properties has led to the separation of the geologic column of interest into several discrete geologic units, each of which is internally more “homogeneous” in some identifiable manner than subdivisions based on other criteria. An interpretation that the material properties of the rocks are controlled principally by the original genetic geologic processes and that much of the post-depositional alteration that has produced second-order variability in properties occurred before tectonic tilting and faulting suggests that the influence of such deformation should be discounted in the modeling process. A concept has been adopted of using a stratigraphic coordinate system during the modeling process, as distinct from a real-world coordinate system that describes the present-day location of points within the several geologic units. Both the measured data and the geostatistical models are described in stratigraphic coordinates throughout the modeling exercise. Models are back-transformed into real-world coordinates upon completion. Because measurements of most material properties of direct interest for use in performance calculations are quite limited in number and spatial distribution, the concept of porosity-as-asurrogate is adopted in order to use relatively abundant and widely distributed (in three dimensions) porosity data as a first approximation of the geologic heterogeneity of the site. The conceptual assumption (Section 5) is that the spatial variation and trends in many different properties should be roughly the same, as the entire ensemble of material properties is related to the geologic processes responsible for formation of these rocks. Sequential gaussian simulation techniques are used to generate the discretized, exhaustive, primary porosity models, conditioned to both laboratory core and downhole petrophysical measurements in an effort to integrate virtually all available data from the entire Yucca Mountain site. Conditional simulation methods produce models that effectively reproduce both the conditioning measurements and the overall (geo)statistical character of those measurements. A technique known as linear coregionalization is then used to produce similarly exhaustive models of selected derivative material properties. For certain material properties, additional postprocessing of these derivative models is required to induce desired geologic attributes in the final models. For example, an indicator-kriged model of hydrous-phase mineral alteration is used Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 47 of 150 to constrain the derivative models of saturated hydraulic conductivity for two units where the presence of this type of alteration is known. The modeling process has been conducted using a scientific notebook procedure (SNL QAIP 20-2, Rev. 02). Details of the modeling process, including expansion of topics discussed only in abbreviated or illustrative form in this report, may be found in the scientific notebooks listed in Table 11. 6.4 DEVELOPMENT OF THE MODELS 6.4.1Model Domain The location of the 3-D rock properties model is shown in relationship to the general Yucca Mountain area and to the Exploratory Studies Facility (ESF) in Figure 2. The modeled region extends from west to east from the vicinity of Fatigue Wash to the middle of Midway Valley, and from north to south from central Yucca Wash to the middle portion of Dune Wash. The model domain for the rock properties model is closely tied to the location of the YMP UZ site-scale model numerical process-model grid. 6.4.2 Separate Modeling of Distinctive, Aggregate Geologic Units Rock properties models have been created for four distinctly different aggregate geologic units: the upper Paintbrush nonwelded interval, the welded portion of the Topopah Spring Tuff, the Calico Hills nonwelded interval, and the Prow Pass Tuff. The relationship of these four aggregate modeling units to the detailed lithostratigraphic subdivision of Yucca Mountain and to other modeling units involved in the GFM and the mineralogical model of Yucca Mountain is shown in Table 10. Note that the aggregate geologic units selected for separate rock properties modeling efforts do not coincide in general with the breaks between genetic “packages” of rock, which at Yucca Mountain are typically collections of virtually coeval pyroclastic flow deposits associated with a major volcanic event such as a caldera-collapse sequence. However, as documented by Rautman and McKenna (1997, pp. 6 to 11, Figures 3 and 4), the available measurements of material properties indicate that the modeling units defined here are more “homogeneous” (consistent) internally than are the major genetic packages. Table 11. Scientific Notebooks for Rock Properties Models Scientific Notebooks Title Reference SNL-SCI-006 Scientific Notebook, Geostatistical Modeling of Porosity and Derivative Properties (FY 1998) Rautman and McKenna 1998 SNL-SCI-011 Scientific Notebook, Geostatistical Modeling of Porosity and Derivative Properties (FY 1999) Rautman 1999 Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 48 of 150 6.4.3 Available Data and Preliminary Processing The data used in modeling the spatial variability of material properties in the rock properties model were obtained from a number of sources, as described in Section 4, including both laboratory measurements of core samples (Sections 4.1.1 and 4.1.3) and down-hole petrophysical measurements (Sections 4.1.2 and 4.1.5) of in-situ rocks. Only surface-based drill holes have been used in constructing the rock properties model, as the large areal extent of the model and the requirement that the vertical positions of all data be expressed in stratigraphic coordinates effectively precludes the use of samples from the underground workings of the ESF. The locations of the various drill holes used in modeling each separate geologic unit are shown in Figures 4 through 7. Note that although there is major consistency of the drill hole coverage from unit to unit, the suite of holes that contain data relevant to any particular model unit is unique. Data tracking numbers associated with all data values are tabulated in Section 4. 6.4.3.1Data Compilation, Re sampling, an d G eneration of De rived Q uantities Compilation of porosity, secondary property data, and mineralogical analyses into initial working files, even though from a large number of sources, is relatively straightforward. Laboratory core data were used essentially as-is. Relevant petrophysical variables were extracted from much more massive source files, containing all extant downhole geophysical data, using software routine BUD (Table 2). Core and petrophysical porosity data were compiled on an individual drill hole basis in spreadsheets using the graphics package Sigma Plot (Section 3), prior to being combined on a byunit basis in Excel. Because the petrophysical data were recorded in the field on 0.5-ft depth spacings (some older holes were recorded on only 1-ft spacings), it was determined that direct use of all available petrophysical values would statistically overwhelm the available core porosity measurements (collected on a nominal 3-ft spacing) available for an equivalent length of drill hole. Therefore, to maintain approximate parity of statistical mass across the two different types of porosity data, the petrophysical values were resampled to a 3-ft depth spacing as well. In keeping with the practice of Rautman and McKenna (1997, pp. 29 to 30, Figure 14) in developing RPM2, the resampling algorithm also calculated a simple average of the adjoining measurements within plus and minus 1 ft of the nominal 3-ft depth value. This averaging process is necessary to deal with “missing” values present in the original petrophysical data (generally related to unacceptable borehole conditions). Both resampling and computation of the desired two-foot average values was accomplished using software routine TWOFOOT (Table 2). Relevant mineralogical analyses (hydrous-mineral phase “altered” mineral fractions of smectite, clinoptilolite, mordenite, and chabazite) were extracted from the source data files (Section 4.1.4) using the data manipulation capabilities of Excel. A new variable, total hydrous-phase mineral content, was created as the sum of these individual mineral fractions. A hydrous-mineral phase alteration-indicator variable, I, was also constructed for the mineralogical data set by the logical statement: (Eq. 1) I 1 Altered Minerals ( ) . 5-percent = . 0 Otherwise . . . . = Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 49 of 150 Figure 4. Drill holes used in modeling the PTn model unit. Note: some holes may include only partial penetrations of the unit. Nevada State Plane Easting, in feet 550000 560000 570000 580000 Nevada State Plane Northing, in feet 740000 750000 760000 770000 780000 790000 UZ Site-Scale Model Boundary RPM3 G-1 G-2 G-3 G-4 H-1 H-3 H-4 H-5 H-6 NRG-4 NRG-5 NRG-6 NRG-7 ONC-1 SD-7 SD-9 SD-12 UZ-1,14 UZ-4,5 UZ-6 UZ-7a UZ-16 WT-1 WT-2 WT-4 WT-7 WT-10 WT-11 WT-12 WT-13 WT-14 WT-15 WT-16 WT-17 WT-18 WT-24 UZN11 UZN31,32 UZN33,34 UZN37 UZN53,54,55 Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 50 of 150 Figure 5. Drill holes used in modeling the TSw model unit. Note: some holes may include only partial penetrations of the unit. Nevada State Plane Easting, in feet 550000 560000 570000 580000 Nevada State Plane Northing, in feet 740000 750000 760000 770000 780000 790000 UZ Site-Scale Model Boundary RPM3 G-1 G-2 G-3 G-4 H -1 H-3 H -4 H-5 H -6 NRG-4 NRG-5 NRG-6 NRG-7 ONC-1 p#1 SD-7 SD-9 SD-12core UZ-1 UZ-4,5 UZ-6 UZ-7a UZ-14 UZ-16 WT-1 WT-2 WT-3 WT-4 WT-7 WT-10 WT-11 WT-12 WT-13 WT-14 WT-15 WT-16 WT-17 W T-18 WT-24 UZN37 UZN53,54,55 UZN5 7,58 ,59 ,61 Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 51 of 150 Dr ill Holes in the CHn Nevada State P lane Easting, in feet 550000 560000 570000 580000 Nev ada State Plane Northing, in feet 740000 750000 760000 770000 780000 790000 UZ Site-Scale Model Boundary RPM3 G-1 G-2 G-3 G-4 H -1 H-3 H -4 H-5 H-6 NRG-7 ONC-1 p#1 S D -7 S D -9core SD-12 UZ-6 UZ-14 UZ-16 WT-1 WT-2 WT-3 WT-4 WT-7 WT-12 WT-14 WT-15 WT-16 WT-17 WT-18 WT-24 Figure 6. Drill holes used in modeling the CHn model unit. Note: some holes may include only partial penetrations of the unit. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 52 of 150 Figure 7. Drill holes used in modeling the Tcp model unit. Note: some holes may include only partial penetrations of the unit. Nevada State Plane Easting, in feet 550000 560000 570000 580000 Nevada State Plane Northing, in feet 740000 750000 760000 770000 780000 790000 UZ Site-Scale Model Boundary RPM3 G-1 G-2 G-3 G-4 H-1 H-3 H-4 H-5 H-6 p#1 SD-7 SD-9 SD-12 UZ-6 UZ-14 UZ-16 WT-2 WT-3 WT-7 WT-12 WT-17 WT-18 Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 53 of 150 In preparation for calibrating the soft indicators of hydrous-phase mineral alteration (Section 6.2.4), additional variables were created for both the core and petrophysical data sets. These variables, referred to as “bound-water,” were calculated simply as the arithmetic difference between the OD and relative humidity (RH) porosity values (for core; Section 4.1.1) or between the POROTOT and POREF porosity values (for petrophysics; Section 4.1.5). Further manipulation and use of the bound-water variables is discussed below in Section 6.4.7.1. 6.4.3.2 Stratigraphic Coordinate Conversion Each major lithologic interval selected for modeling has been modeled in a stratigraphic coordinate system that reflects the original, pre-faulting depositional continuity of these ash-flow and air-fall tuffaceous deposits (Figure 8). Stratigraphic coordinates use the same east-west and north-south coordinates (Nevada state plane coordinate system, defined in feet) as the drill hole from which the relevant data were obtained. However, the vertical coordinate of a sample is represented as the relative fractional position of that sample within the thickness of the entire unit at that horizontal location. The stratigraphic coordinate concept effectively removes the effect both of depositional thinning away from the source volcanic vent(s) and of post-depositional tilting and deformation, and it thus positions samples from equivalent portions of the overall unit at the same nominal internal position within a rectangular volume. As shown schematically in part (a) of Figure 8, regions of varying material properties are presumed to have been emplaced or otherwise formed by various alteration processes in an essentially stratiform manner. At Yucca Mountain, the volumetrically dominant rocks were formed by deposition by pyroclastic flows to form thick ash-flow sheets that thin laterally away from their source. Thus, there is a tendency for regions of somewhat similar material properties that were formed under somewhat similar pressure-temperature conditions to occupy roughly the same relative vertical position within a unit. Later faulting as part of Basin and Range tectonism disrupted the originally continuous volcanic rocks and tilted the rock units, with their contained material properties, toward the east, as indicated in part (b) of the figure. Modeling of those rock properties is illustrated in part (c) of Figure 8. The vertical locations of drill hole samples are specified within the stratigraphic coordinate system as a fractional distance where the base of the unit is assigned a distance of zero and the top of the unit is assigned a distance of one. Stratigraphic coordinates are thus dimensionless. As suggested by the mesh of intersecting dotted lines in the right-hand portion of Figure 8, a regular rectangular modeling grid is defined within each stratigraphic coordinate system. Because the various material property zones have been stretched or compressed vertically so that the overall stratigraphic thickness of the unit is constant, defining the modelling grid within this framework generally positions nodes with similar materials on a stratigraphically “horizontal” plane. This repositioning of similar materials in similar relative locations greatly simplifies the search for data in the neighborhood of an unsampled location, as shown conceptually by the search ellipse in part (c) of Figure 8. Although it is possible to rotate the principal direction of the search ellipse to match the overall tectonic dip of the unit (see part (b) of the figure), it is virtually impossible to modify the search strategy to account for vertical displacement of the material property zones by discrete faults. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 54 of 150 Figure 8. Conceptual illustration of the construction and use of stratigraphic coordinates. (a) Rock unit is formed by areally extensive volcanic (or sedimentary) processes. Zones of differing rock properties (shaded colors) are formed in a stratiform manner. (b) Tectonic deformation tilts and disrupts original stratiform continuity by faulting. (c) Modeling unit is returned to an approximation of original continuity in a rectangular coordinate system in which all vertical distances are measured as a fractional position measured from the top or bottom of the rock unit. (a) (b) (c) Original Depositional Continuity Conversion to Stratigraphic Coordinates Ash-Flow Deposition 0 1 Tectonic Tilting and Faulting 356' 932' Search Ellipse Structure- Contour Horizon Isochore Information 7356' 528' Drill Hole Search Ellipse Stratcoords.cdr Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 55 of 150 At the end of a modeling exercise, the transformation process between parts (b) and (c) in Figure 8 is reversed by assigning each grid node a computed vertical position derived from knowledge of the structure contour model for the top of each unit and the spatially varying thickness of each unit (software routine COORDS, Table 2). These values are obtained from the independently developed 3-D GFM (e.g., GFM3). In practice, implementation of the stratigraphic-coordinate concept is slightly more complicated than the idealized example of Figure 8. First, sample locations are typically specified in terms of their depth within a specific drill hole (the drilling process measures all locations from the collar of the hole regardless of the physical elevation of the hole and its contained samples). Thus, the measured depths were converted to stratigraphic depths initially, and only to stratigraphic elevations at the time of modeling. Second, for reasons involving principally numerical precision within the computer programs that implement the actual rock properties modeling algorithm(s), the fractional stratigraphic positions indicated in Figure 8 are multiplied by unit-specific scaling constant to obtain values that approximate the nominal thickness of the different units in the real world. Additionally, unlike the two-dimensional example shown in Figure 8, actual modeling of rock properties was conducted in full 3-D space. Finally, the issue arises regarding how to treat samples from a drill hole that fails to penetrate the entire thickness of the geologic unit in question (represented by the drill hole on the left-hand side of Figure 8). Clearly it is inappropriate to assign a stratigraphic elevation of zero to a sample obtained from the very bottom of the hole itself, as the materials at this elevation in general are not representative of materials at the base of the unit here or elsewhere. Yet, without drilling deeper, the distance between the foot of the hole and the true base of the unit is unknown. Such situations have been reconciled by inferring the base of the unit in question from the GFM (e.g., GFM3) and adjusting the fractional position accordingly. The presumption is that the base of the unit projected using the framework model is a reasonable approximation of the unknown true position at that location. Similar reconstructions are required in the case where erosion has removed the top of a particular modeling unit at a drill hole location or if part of a unit is missing from a drill hole because of faulting. The equation for calculating stratigraphic depth (StratDepth) is: , (Eq. 2) where SampleDepth is the depth of the relevant sample measurement in a given drill hole, UnitTop and UnitBottom are the measured or projected depths to the top and bottom contacts of that model unit at the drill hole location used to determine the thickness of the unit, and NominalThickness is the appropriate unit-specific scaling constant from Table 12. Conversion of drill hole depths to stratigraphic depths was accomplished using the transform capabilities of Sigma Plot (STRATC4.XFM; Table 2; Attachment II). Stratigraphic elevation is calculated simply as NominalThickness–StratDepth, and was calculated using the mathematical capabilities of Excel prior to writing final ASCII text files of the data values. StratDepth SampleDepth UnitTop – UnitBottom UnitTop – --------------------------------------------------------------- NominalThickness · = Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 56 of 150 6.4.4 Statistical and Spatial Description of Porosity Porosity is the primary material property used to generate the material properties models contained in RPM3. Even the models of derivative properties such as hydraulic conductivity are based ultimately on porosity. This section presents the statistical description (the ensemble characteristics; Sections 6.3.1.3 and 6.3.1.4) of the underlying porosity measurements that have been used to condition the stochastic models of each modeling unit. This presentation consists principally of histograms for the various model units and a summary of their descriptive statistics (Table 13). Experimental variograms calculated using the measured values and a fitted model for each unit are presented to represent the spatial continuity patterns expressed in the data. Finally, because of the unique geologic nature of the TSw model unit, two different types of porosity (matrix and lithophysal) are presented and justified. 6.4.4.1P Tn Model Unit Porosity Data–Porosity data obtained from the upper Paintbrush nonwelded model unit are portrayed in histogram format in Figure 9. A statistical summary of these data is given in Table 13. Spatial Continuity Data–Modeled variograms for porosity in the PTn model unit are presented in Figure 10. Parameters of the fitted variogram model are presented in Table 14. Details of the variogram modeling exercise for the PTn model unit may be found in the scientific notebook (Rautman and McKenna 1998, pp. 361 to 368). Table 12. Unit-Specific Scaling Constants for Nominal Model Unit Thickness Model Unit Nominal Thickness PTn 200 TSw 1000 CHn 400 Tcp 400 Table 13. Statistical Summary of Porosity Data Used in Modeling PTn TSw CHn Tcp Matrix Lithophysal Mean 0.342 0.122 0.146 0.254 0.229 Std.Dev. 0.122 0.061 0.076 0.075 0.090 Minimum 0.001 0.001 0.001 0.001 0.001 Maximum 0.750 0.519 0.551 0.511 0.494 N 2,300 11,052 11,796 2711 2603 Note: All values are porosity as a fraction (equivalent to m3/m3), except number of values (N). Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 57 of 150 6.4.4.2 TSw Model Unit Porosity Data–Porosity values obtained from the Topopah Spring welded (TSw) model unit are presented graphically in histogram format in Figure 11; the corresponding statistical summary of these data is presented in Table 13. Examination of the raw porosity data indicates that there are two different porosity values of interest in modeling rock material properties: matrix and lithophysal. Lithophysal porosity, as that term is used in this report, is taken to mean the porosity of volumes of rock many tens of centimeters in diameter, such that the porosity effect of large (decimeter scale and larger) lithophysal cavities is included. In contrast, the term matrix porosity is used in this report to refer to the porosity approximately equivalent to that measured for Figure 10. Experimental and fitted model variograms of matrix porosity in the PTn model unit in the (a) stratigraphic vertical and (b) stratigraphic horizontal directions. Separation Distance, in feet 0 50 100 150 200 Normal-Score Variogram 0.0 0.5 1.0 1.5 Data a priori Variance Model PTnVertVar.wmf (a) Separation Distance, h, in feet 0 5000 10000 15000 Normal-Score Variogram 0.0 0.5 1.0 1.5 Data, Az. = 135 Data, Az. = 45 a priori Variance Model, Az. = 135 Model, Az. = 45 PTnHorizVar.wmf (b) Porosity, as a fraction 0.0 0.2 0.4 0.6 0.00 0.02 0.04 0.06 0.08 0.10 0.0 0.2 0.4 0.6 0.8 1.0 Frequency Cumulative Frequency Figure 9. Histogram and cumulative distribution function of matrix porosity in the PTn model unit. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 58 of 150 laboratory core samples, in which the size of the matrix pores is small enough that water is held in them by capillarity under slightly unsaturated (negative pressure) conditions. The difference between the two types of porosity measurements is not trivial, as illustrated in Figure 12, a comparative down-hole plot of matrix and lithophysal porosity data from drill hole SD-7. Lithophysal porosity is indicated by the dark solid curve, whereas the matrix porosity values measured for core samples are indicated by the filled-circle symbols. Note the marked divergence of the porosity values indicated by these two sets of data in two vertical locations within the drill hole. In general, these zones of divergence correspond to the two lithophysal zones (upper and lower) defined by Buesch et al. (1996). However, the correspondence is not at all exact, and the petrophysical POREF (= lithophysal) porosity curve from downhole geophysics indicates that substantial lithophysal cavity development may extend significantly above and below the limits Table 14. Variogram Parameters for Spatial Continuity Model of Porosity in the PTn Model Unit Nest No. Model Type Range (ft) Sill Rotation Angle (degrees)1 Anisotropy Ratio1 Maximum (horizontal) Intermediate Minimum (vertical) 1 2 3 1 2 -- Nugget -- -- -- 0.01 -- -- -- -- -- 1 Spherical 300 50 20 0.20 135 135 135 0.167 0.067 2 Spherical 5,000 2,500 80 0.55 135 135 135 0.500 0.016 3 Spherical 10,000 7,500 100 0.24 135 135 135 0.750 0.010 Note: 1. Rotation angles 1–3 and anisotropy ratios 1–2 are keyed to input requirements of the SGSIM computer code, as documented in Deutsch and Journel 1992, pp. 167 and 22 to 29. Figure 11. Histograms and cumulative distribution functions of (a) matrix porosity and (b) lithophysal porosity in the TSw model unit. y Porosity, as a fraction 0.0 0.2 0.4 0.6 Frequency 0.00 0.05 0.10 0.15 0.20 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 TSwM_RHpor.wmf (a) Porosity, as a fraction 0.0 0.2 0.4 0.6 Frequency 0.00 0.05 0.10 0.15 0.20 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 TSwL_RHpor.wmf (b) Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 59 of 150 Porosity, fraction 0.0 0.2 0.4 0.6 Depth, in feet 400 500 600 700 800 900 1000 1100 1200 1300 1400 POREF VWC Core Porosity Tptrn Tptrl Tptpul Tptpmn Tptpll Tptpv3 Tptpln Tptpv1 Tptpv2 Tptbt1 Tptrv3 Figure 12. Comparison of different types of porosity data for drill hole USW SD-7. Horizontal dashed lines indicate lithostratigraphic units as defined by Buesch et al. (1996). Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 60 of 150 of the formally named lithophysal zones (Figure 12). A scatterplot of these lithophysal porosity values versus the depth-equivalent core RH porosity for the formally named lithophysal-zone intervals only is presented in Figure 13, part (a). The region of marked divergence from the 45° one-to-one correspondence line represents high lithophysal porosities matched on a nearestsample basis to the lower core porosity values. Note also the third set of data plotted in Figure 12; this curve is identified as volumetric water content (equivalent to the amount of water-filled porosity) and is shown by the dashed line without symbols. This third curve is observed essentially to overlie the true matrix (= core) porosity data throughout much of the drill hole. In locations where the water-filled porosity trace does not closely match the core values, it is always observed to indicate markedly lower values than the solid lithophysal porosity curve. The near-equivalence of core and VWC values is demonstrated in part (b) of Figure 13. The small set of mismatched pairs may be identified in Figure 12 as representing the interval between approximate 480 and 550 feet. Because many drill holes lack core data from which to obtain matrix porosity for modeling purposes, the practice has been adopted (for those non-cored holes only) of using volumetric water content as a surogate for matrix porosity. Outside these lithophysae-bearing intervals, matrix porosity is set equal to POREF porosity, whereas within these intervals of curve separation, matrix porosity is set equal to the VWC values. It is clear that this practice is merely a simple heuristic device, and that use of the volumetric water content as representing water-filled porosity unquestionably will underestimate the actual matrix porosity within the UZ for the simple reason that all the available matrix pore space in the UZ is not water filled (any large lithophysal cavities also will not be water filled). However, the approach is a reasonable approximation for several reasons. (1) Water saturation throughout much of the Topopah Spring Tuff is rather high, typically greater than about 80 percent. (2) The VWC p y y Core Porosity, as a fraction 0.0 0.1 0.2 0.3 0.4 Petrophysical Porosity (POREF), as a fraction 0.0 0.1 0.2 0.3 0.4 Figure 13. Crossplots of (a) petrophysical POREF porosity vs. core RH porosity and (b) petrophysical VWC vs. core RH porosity for formally defined lithophysal intervals in drill hole SD-7. Line indicates one-to-one correspondence. p y y Core Porosity, as a fraction 0.0 0.1 0.2 0.3 0.4 Petrophysical Porosity (VWC), as a fraction 0.0 0.1 0.2 0.3 0.4 (a) (b) Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 61 of 150 values in zones of significant lithophysal cavity development are much closer to the true matrix porosity than are the lithophysal porosity values, which can be observed from Figures 12 and 13 to be as much as double the matrix (core) porosity values for drill hole USW SD-7. And finally, (3) because one of the purposes of modeling the spatial heterogeneity of porosity is to attempt to model the variability of saturated hydraulic conductivity using porosity as a surrogate, the volumetric water content representing the water-filled porosity (at any in-situ saturation) most likely represents essentially all the pore space that is available for the transmission of water under unsaturated conditions. Note that it is possible to significantly underestimate the whole-rock saturated hydraulic conductivity of the Topopah Spring welded interval under true saturated conditions in this manner, because the effect of lithophysal cavities on the whole-rock permeability is neglected. However, such underestimation is likely to be a non-issue, as this relatively high stratigraphic interval almost invariably is present in the unsaturated zone throughout the entire modeled region. Additionally, under true saturated conditions, most ground-water flow through the TSw model unit would be through fractures (and lithophysae), and not through the matrix. Spatial Continuity Data–Modeled variograms for both matrix and lithophysal porosity in the TSw model unit are presented in Figure 14. Parameters for both fitted variogram models are presented in Table 15. Details of the variogram modeling exercise may be found in the scientific notebook (Rautman and McKenna 1998, pp. 437 to 444 and 459 to 464). 6.4.4.3 CHn Model Unit Porosity Data–Porosity values obtained from the Calico Hills nonwelded (CHn) model unit are presented graphically in Figure 15; the corresponding statistical summary of these data is presented in Table 13. Spatial Continuity Data–A modeled variogram for matrix porosity in the CHn model unit is presented in Figure 16. Parameters of the fitted variogram model are presented in Table 16. Details of the variogram modeling exercise may be found in the scientific notebook (Rautman and McKenna 1998, pp. 534 to 556). 6.4.4.4 Tcp Model Unit Porosity Data–Porosity values obtained from the Prow Pass Tuff (Tcp) model unit are presented graphically in histogram format in Figure 17; the corresponding statistical summary of these data is presented in Table 13. Spatial Continuity Data–A modeled variogram for matrix porosity in the Tcp model unit is presented in Figure 18. Parameters of the fitted variogram model for porosity in the Tcp model unit are presented in Table 17. Details of the variogram modeling exercise may be found in the scientific notebook (Rautman and McKenna 1998, pp. 704 to 726). Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 62 of 150 6.4.5 Statistical Description of Secondary Material Properties The secondary, or derivative, material properties considered in this analysis are matrix saturated hydraulic conductivity, thermal conductivity, and bulk density (Section 4.1.3). In general, there are insufficient individual measurements of these material properties to determine their spatial continuity patterns independently of porosity. Accordingly, they have been modeled under the presumption that each variable is coregionalized with porosity. This section contains the basic statistical description of each derivative material property, including the correlation of that property with corresponding porosity measurements, which is used in the coregionalization process. 6.4.5.1Matr ix S aturated H ydraulic Con ductivity The histogram of all available laboratory measurements of matrix saturated hydraulic conductivity is presented in Figure 19(a). A scatterplot of these same data against the Figure 14. Experimental and fitted model variograms for porosity in the TSw model unit. (a) Matrix porosity, stratigraphic vertical and (b) stratigraphic horizontal. (c) Lithophysal porosity, stratigraphic vertical and (d) stratigraphic horizontal. Separation Distance, in feet 0 100 200 300 400 500 600 Normal-Score Variogram 0.0 0.5 1.0 1.5 Data a priori variance Model TSwMVertVar.wmf (a) Separation Distance, h, in feet 0 5000 10000 15000 20000 25000 Normal-Score Variogram 0.0 0.5 1.0 1.5 Data N-S Data E-W a priori variance Model N-S Model E-W TSwMHorizVar.wmf (b) Separation Distance, in feet 0 100 200 300 400 500 600 Normal-Score Variogram 0.0 0.5 1.0 1.5 Data a priori variance Model TSwLVertVar.wmf (c) Separation Distance, h, in feet 0 5000 10000 15000 20000 25000 Normal-Score Variogram 0.0 0.5 1.0 1.5 Data N-S Data E-W a priori variance Model N-S Model E-W TSwLHorizVar.wmf (d) Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 63 of 150 corresponding laboratory porosity values is shown in part (b) of the figure, and a statistical summary of the data is presented in Table 18. Table 15. Variogram Parameters for Spatial Continuity Model of Porosity in the TSw Model Unit Nest No. Model Type Range (ft) Sill Rotation Angle (degrees)1 Anisotropy Ratio1 Maximum (horizontal) Intermediate Minimum (vertical) 1 2 3 1 2 TSw Model Unit, Matrix Porosity -- Nugget -- -- -- 0.01 -- -- -- -- -- 1 Spherical 250 250 20 0.19 0 0 0 1.000 0.080 2 Spherical 2,000 6,500 90 0.25 0 0 0 3.250 0.045 3 Spherical 29,000 12,500 350 0.55 0 0 0 0.431 0.012 TSw Model Unit, Lithophysal Porosity -- Nugget -- -- -- 0.01 -- -- -- -- -- 1 Spherical 600 400 15 0.19 0 0 0 0.667 0.025 2 Spherical 5,000 8,000 175 0.20 0 0 0 1.600 0.035 3 Spherical 38,000 17,000 200 0.60 0 0 0 0.447 0.006 Note: 1. Rotation angles 1–3 and anisotropy ratios 1–2 are keyed to input requirements of the SGSIM computer code, as documented in Deutsch and Journel 1992, pp. 167 and 22 to 29. Figure 15. Histogram and cumulative distribution function of matrix porosity in the CHn model unit. Porosity as a factor 0.0 0.2 0.4 0.6 Frequency 0.00 0.05 0.10 0.15 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 CHnRHpor.wmf Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 64 of 150 Figure 19(a) is distinctively multi-modal, and Figure 19(b) indicates four subpopulations that can be identified in the matrix hydraulic conductivity data set. Essentially, cluster “A” consists of samples for which the delta-porosity, computed simply as the arithmetic difference between the OD and RH laboratory-measured porosity values, is less than 5 percent (0.05). Cluster “A” corresponds to a continuum of hydraulic conductivity values that are rather strongly correlated to porosity, as presented separately in Figure 20. Rautman and McKenna (1997, pp. 37 to 40) present an expanded discussion of identifying hydrous-phase mineral alteration in core samples through the use of delta-porosity values. Cluster “B” in Figure 19(b) contains samples for which the delta-porosity between the OD and RH porosity measurements is greater than 5 percent. These are samples that exhibit significant hydrous-phase mineral alteration, as evidenced by the removal of moderately large quantities of water during drying at 105°C. This subpopulation is plotted separately in Figure 21, and it appears largely indistinguishable from a normally distributed (gaussian) population if the small group of Table 16. Variogram Parameters for Spatial Continuity Model of Porosity in the CHn Model Unit Nest No. Model Type Range (ft) Sill Rotation Angle (degrees)1 Anisotropy Ratio1 Maximum (horizontal) Intermediate Minimum (vertical) 1 2 3 1 2 -- Nugget -- -- -- 0.05 -- -- -- -- -- 1 Spherical 3,000 3,000 30 0.25 1.000 0.010 2 Spherical 12,000 5,000 400 0.70 0.417 0.033 Note: 1. Rotation angles 1–3 and anisotropy ratios 1–2 are keyed to input requirements of the SGSIM computer code, as documented in Deutsch and Journel 1992, pp. 167 and 22 to 29. Figure 16. Experimental and fitted model variograms for matrix porosity in the CHn model unit. (a) Stratigraphic vertical; (b) stratigraphic horizontal. Separation Distance, in feet 0 50 100 150 200 250 300 350 Normal-Score Variogram 0.0 0.5 1.0 1.5 Data a priori Variance Model CHnVertVar.wmf (a) Separation Distance, h, in feet 0 5000 10000 15000 Normal-Score Variogram 0.0 0.5 1.0 1.5 Data N-S Data E-W a priori Variance Model N-S Model E-W CHnHorizVar.wmf (b) Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 65 of 150 samples with conductivities greater than 10–6 m/sec are discounted as potentially misclassified samples. Cluster “C” of Figure 19(b) is represents a small subset of the entire population of matrix hydraulic conductivity measurements for which the corresponding porosity measurement is very low, and for which there appears to be no particular correlation of porosity with conductivity. These samples are almost exclusively from two vitrophyric lithostratigraphic units: Tptrv1 and Tptpv3, Figure 22(a) and (c). Experimentation indicates that almost the identical subset of samples can be extracted from the complete hydraulic conductivity data set simply as those samples whose measured porosity is less than 5 percent (0.05). These samples are presented in Figure 22(b) and (d). Flint (1998, p. 38) interpreted the relatively uniform-random distribution of these “matrix” hydraulic conductivities as caused by microfracturing within these densely welded rock types. Figure 17. Histogram and cumulative distribution function of matrix porosity in the Tcp model unit. Porosity as a factor 0.0 0.2 0.4 0.6 Frequency 0.00 0.05 0.10 0.15 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 TcpRHpor.wmf Figure 18. Experimental and fitted model variograms for matrix porosity in the Tcp model unit. (a) Stratigraphic vertical; (b) stratigraphic horizontal. Separation Distance, in feet 0 50 100 150 200 250 300 350 Normal-Score Variogram 0.0 0.5 1.0 1.5 Data a priori Variance Model TcpVertVario.wmf (a) Separation Distance, h, in feet 0 5000 10000 15000 Normal-Score Variogram 0.0 0.5 1.0 1.5 Data N-S Nata E-W a priori Variance Model N-S Model E-W TcpHorizVario.wmf (b) Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 66 of 150 Cluster “D” of Figure 19(b) represents those core samples for which the matrix saturated hydraulic conductivity was below the sensitivity limit of the laboratory permeameter. Note that these null values have been ignored in all statistical calculations for this modeling exercise. As shown in Figure 19, these “no-flow” samples have been arbitrarily set equal to 10–14 m/sec (log10 = –14.0) for plotting purposes, both here and elsewhere in this report. In contrast to the discussion of Figures 20 and 21, which involves more or less global characteristics of the measured matrix saturated hydraulic conductivity measurements from Yucca Mountain, a broader objective of this modeling exercise is to reproduce the material property characteristics appropriate for individual model units. The ability to discriminate and identify the hydraulic conductivity characteristics of specific vitrophyric rock units essentially by looking Table 17. Variogram Parameters for Spatial Continuity Model of Porosity in the Tcp Model Unit Nest No. Model Type Range (ft) Sill Rotation Angle (degrees)1 Anisotropy Ratio1 Maximum (horizontal) Intermediate Minimum (vertical) 1 2 3 1 2 -- Nugget -- -- -- 0.05 -- -- -- -- -- 1 Spherical 6,000 4,000 80 0.40 0.667 0.013 2 Spherical 15,000 10,000 800 0.55 0.667 0.053 Notes: 1. Rotation angles 1–3 and anisotropy ratios 1–2 are keyed to input requirements of the SGSIM computer code, as documented in Deutsch and Journel 1992, pp. 167 and 22–29. , log10 (Ksat, m/sec) -14 -12 -10 -8 -6 -4 Frequency 0.00 0.02 0.04 0.06 0.08 0.10 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 Ksat_alldata_hist.wmf Figure 19. (a) Histogram of laboratory-measured matrix saturated hydraulic conductivity from core samples and (b) crossplot of hydraulic conductivity vs. matrix porosity for the same samples. Porosity, as a fraction 0.0 0.2 0.4 0.6 log10(Saturated Hydraulic Conductivity, in m/sec) -15 -10 -5 (b) A B C D (a) Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 67 of 150 solely at their porosity characteristics (e.g., Figure 20 and discussion) is one example. Other relevant matrix hydraulic conductivity characteristics for purposes of this report are the by-unit statistical distributions of values. These target histograms for the individual simulated models of the four modeling units are presented in Figure 23 (a) through (d). 6.4.5.2 Bulk Density A histogram of more than 3,000 relative-humidity-oven-dried bulk density values from all four model units is presented in Figure 24(a). The composite histogram is distinctly multimodal, as might be expected from aggregation of densely welded tuffs with a collection of non- to partially welded ash-flow tuffs and nonwelded air-fall tuffaceous debris. Table 18. Statistical Summary of Matrix Saturated Hydraulic Conductivity Measurements All Data Unaltered Altered Mean -8.9553 -8.6536 -9.6805 Std.Dev. 1.7317 1.7248 1.5269 Minimum -11.7086 -11.3360 -11.7086 Maximum -4.6866 -4.6955 -4.6866 N 405 286 119 Note: All values in log10 m/sec except number of samples (N); “no-flow” samples omitted. Figure 20. (a) Histogram of laboratory-measured matrix saturated hydraulic conductivity from unaltered core samples (identified by delta-porosity) and (b) crossplot of hydraulic conductivity vs. matrix porosity for the same samples. Conductivity of non-flowing samples set to 10-14 m/sec for plotting only. log10 (Ksat, m/sec) -14 -12 -10 -8 -6 -4 Frequency 0.00 0.02 0.04 0.06 0.08 0.10 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 Ksat_unalt_hist.wmf (a) Porosity, as a fraction 0.0 0.2 0.4 0.6 Matrix Saturated Hydraulic Conductivity, in log10(m/sec) -14 -12 -10 -8 -6 -4 -2 Unaltered Regression (-14.0 = null) r2 = 0.604 r = +0.777 Ksat_unalt_xplot.wmf (b) Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 68 of 150 Despite the compositing of a group of diverse lithologies (and mineralogies) for Figure 24(a), the bulk density data appear to belong to a relatively consistent “family” of values. A scatter diagram showing the relationship between the bulk density values and RH-dried porosity is presented in Figure 24(b). Although there is some segregation of the cloud of data points into two overlapping subpopulations (which can be demonstrated to relate to hydrous-phase mineral alteration), the coefficient of determination (r2) value is a very strong 0.912, indicating only somewhat less than one-to-one correspondence. Given the generally low importance of bulk density in most performance assessment analyses, it has not been determined to be of sufficient importance to analyze altered and unaltered bulk densities separately. Unit-specific histograms that serve as the target distributions for the coregionalized bulk density models are presented in Figure 25. Model unit genesis is reflected clearly in the four parts of the figure: high-density densely welded tuffs in the TSw model unit (part (b)), low-density nonwelded and pumiceous rich materials in the PTn unit (part (a)), and intermediate-density nonwelded to only partially welded tuffaceous materials in the CHn and Tcp model units (parts (c) and (d)). 6.4.5.3 Thermal Conductivity Thermal conductivity and porosity data have been obtained from a modest number of samples collected from the PTn and TSw model units. Figure 26(a) presents a histogram and cumulative distribution function for these thermal conductivity values measured at 70°C and under 105°Cdried saturation conditions (note that thermal conductivity is a function of both absolute temperature and saturation state). The crossplot of these same values in Figure 26(b) clearly Figure 21. (a) Histogram of laboratory-measured matrix saturated hydraulic conductivity from altered core samples (identified by delta-porosity) and (b) crossplot of hydraulic conductivity vs. matrix porosity for the same samples. Conductivity of non-flowing samples set to 10-14 m/sec for plotting only. log10 (Ksat, m/sec) -14 -12 -10 -8 -6 -4 Frequency 0.00 0.05 0.10 0.15 0.20 0.25 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 Ksat_altd_hist.wmf Porosity, as a fraction 0.0 0.2 0.4 0.6 Matrix Saturated Hydraulic Conductivity, in log10(m/sec) -14 -12 -10 -8 -6 -4 -2 Altered -14.0 = null Ksat_altd_xplot.wmf (a) (b) Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 69 of 150 indicates an inverse relationship of thermal conductivity with porosity. A statistical summary of these available thermal conductivity data is presented in Table 19. There are two significant difficulties in using the existing thermal conductivity data, which is derived from laboratory testing on “matrix-sized” specimens. First, the in-situ bulk density of the Topopah Spring welded unit is most directly related to the lithophysal porosity of this unit, not to the matrix porosity. Open cavities too large to be measured as part of the matrix porosity laboratory procedure (or its petrophysical equivalent) will significantly reduce the whole rock bulk density, leading to lower effective thermal conductivity. Previous discussion of porosity data from the TSw model unit (Section 6.4.4.2) indicates that the lithophysal porosity values can be greater by a factor of two compared with the depth-equivalent matrix porosity values. Use of this factor-of-two porosity difference applied to the regression relationship shown in Figure 26(b) Figure 22. Scatterplots of matrix saturated hydraulic conductivity versus matrix porosity for (a) core samples from vitrophyric lithostratigraphic units Tptrv1 and Tptpv3 and (b) core samples with matrix porosity less than 0.05. (c) and (d) Histograms and cumulative distribution functions for samples presented in parts (a) and (b), respectively; no-flow samples omitted. Porosity, as a fraction 0.00 0.05 0.10 0.15 0.20 log10 Sat. Hydraulic Conductivity, in m/sec -14 -12 -10 -8 -6 -4 -2 (a) KsatVitroUnitXplot.wmf Porosity, as a fraction 0.00 0.05 0.10 0.15 0.20 log10 Sat. Hydraulic Conductivity, in m/sec -14 -12 -10 -8 -6 -4 -2 (b) KsatVitroPorXplot.wmf log10 (Ksat, m/sec) -14 -12 -10 -8 -6 -4 Frequency 0.00 0.05 0.10 0.15 0.20 0.25 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 (c) KsatVitroUnit.wmf log10 (Ksat, m/sec) -14 -12 -10 -8 -6 -4 Frequency 0.00 0.05 0.10 0.15 0.20 0.25 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 (d) KsatVitroPor.wmf Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 70 of 150 would lead to overprediction errors of roughly 30 to 50 percent for thermal conductivity, depending on the actual porosity level considered. Given that the thermal conductivity models will be coregionalized using lithophysal porosity in an effort to model the in-situ thermal conductivity of the rock mass, the phrase “whole-rock” thermal conductivity will be used in this report. A second, and rather severe difficulty with the available thermal conductivity data is that the density of sampling is not great (a maximum of 52 samples; see Table 19), and furthermore, those samples are highly biased both spatially and toward low-porosity materials, as follows: (1) Two of four drill holes that were sampled for thermal conductivity specimens, although located within the extended site area, are actually located some distance from the region transected by the main part of the ESF (NRG-4, NRG-5; see Figure 5). (2) The sampling vertically within a drill hole is not systematic, and in fact, the vertical distribution of samples cannot be considered particularly Figure 23. Matrix saturated hydraulic conductivity histograms and cumulative distribution functions for (a) the PTn, (b) TSw, (c) CHn, and (d) Tcp model units. Saturated Hydraulic Conductivity, in log10 m/sec -14 -12 -10 -8 -6 -4 Frequency 0.00 0.05 0.10 0.15 0.20 0.25 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 PTnKsatHist.wmf (a) Saturated Hydraulic Conductivity, in log10 m/sec -14 -12 -10 -8 -6 -4 Frequency 0.00 0.02 0.04 0.06 0.08 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 TSwKsatHist.wmv (b) Saturated Hydraulic Conductivity, in log10 m/sec -14 -12 -10 -8 -6 -4 Frequency 0.00 0.05 0.10 0.15 0.20 0.25 0.30 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 CHnKsatHist.wmf (c) Saturated Hydraulic Conductivity, in log10 m/sec -14 -12 -10 -8 -6 -4 Frequency 0.00 0.05 0.10 0.15 0.20 0.25 0.30 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 TcpKsatHist.wmf (d) Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 71 of 150 “representative” of the entire Topopah Spring welded model unit. (3) The only samples that represent the higher porosity values (needed to model the effect of lithophysal cavities on heat conduction) are taken from the PTn model unit, and thus represent nonwelded tuffs rather than the lithophysal portion of welded materials. With respect to the several biases known to exist in the thermal conductivity data, consider Figure 27 (a), which is a histogram of the 54 porosity values measured on the thermal conductivity test specimens (Table 19). Comparison of this figure with the histogram of all lithophysal porosity values measured from the TSw model unit (part (b) of the figure) clearly indicates the extent of the sampling bias. First, Figure 27(a) is clearly bimodal, representing as it does, samples from both welded and nonwelded rock types. If the group of samples with porosities higher than about 40 percent (presumably representing nonwelded tuffs) is discounted, the mode corresponding to the welded Topopah Spring samples is strongly skewed toward lower porosity values. For example, the approximate modal value in Figure 27(a) is 8 to 10 percent, whereas the modal value for the TSw model unit as a whole (Figure 27(b)) is 18 to 20 percent for lithophysal porosity. An attempt to reduce the impact of these sampling biases has been conducted in the following manner. First, the regression relationship presented in the scatterplot of Figure 26(b) is assumed to be a valid predictor of thermal conductivity across the range of porosity values appropriate for the TSw model unit. A systematic prediction of the nominal thermal conductivity of the TSw model unit is then generated using the systematically sampled (nominal 3-ft spacing) porosity data available for three drill holes located along the ESF main drift (SD-7, -9, and -12; Figure 29). These three sets of predicted values are then aggregated and the appropriate statistical quantities and histogram are computed (Figure 28; Table 20). The histogram of Figure 28 will be used as the target distribution for the coregionalized models. Note that although the mean thermal conductivity for both the 35 measured samples of the Topopah Spring welded unit and the Figure 24. (a) Histogram and cumulative distribution function of relative-humidity-oven-dried bulk density values from all four modeling units. (b) Scatterplot of bulk density as a function of matrix porosity for the same data. Porosity, as a fraction 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Bulk Density, in g/cm3 0 1 2 3 r ² = 0.912 r = -0.955 BulkDporXplot.wmf (b) log10 (Ksat, m/sec) 0 1 2 3 Frequency 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 BulkDhist.wmf (a) Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 72 of 150 systematically predicted thermal conductivity of the TSw model unit as a whole (Table 20) are remarkably similar at 1.2 W/m-K (rounded), the local thermal conductivity of major intervals within the unit differ markedly from one another (Figure 29). 6.4.6 Modeling Techniques This section contains a summarized yet logically accurate description of the geostatistical modeling techniques used to produce the rock properties models, including reference to the various computer codes (Table 1) and software routines (Table 2) used to generate the actual models. Sequential gaussian simulation is used to generate the primary porosity models, whereas linear coregionalization is used to generate the derivative models of secondary properties. Indicator kriging is used to produce the model of hydrous-phase mineral alteration that constrains the distribution of the derivative models of saturated hydraulic conductivity in two of the model units. Figure 25. Bulk density histograms for (a) the PTn, (b) TSw, (c) CHn, and (d) Tcp model units. Bulk Density, in g/cm3 0 1 2 3 Frequency 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 PTnBulkdHist.wmf (a) Bulk Density, in g/cm3 0 1 2 3 Frequency 0.00 0.05 0.10 0.15 0.20 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 TswBulkdHist.wmf (b) Bulk Density, in g/cm3 0 1 2 3 Frequency 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 CHnBulkdHist.wmf (c) Bulk Density, in g/cm3 0 1 2 3 Frequency 0.00 0.05 0.10 0.15 0.20 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 TcpBulkdHist.wmf (d) Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 73 of 150 6.4.6.1Discr etization of the Model Dom ain Details of the grid defined for all geostatistical modeling are presented in Table 21. The horizontal discretization of 200-m by 200-m is tied directly to the horizontal resolution of the unsaturatedzone site-scale flow model across its entire areal extent. Vertical resolution of the geostatistical grid is also tied to the unsaturated-zone flow model, particularly within the general repository volume. Additionally, the total number of grid nodes in the resulting suite of models was computationally tractable. Table 19. Statistical Summary of Measured Thermal Conductivity Data from Non-Zeolitic Rock Samples Thermal Conductivity (W/m-K) at 70°C and 105°C dry Porosity (105°C) All Data TSw Only Mean 1.054 1.241 0.197 Std.Dev. 0.516 0.427 0.156 Minimum 0.160 0.620 0.040 Maximum 2.200 2.200 0.610 N 52 35 54 Note: All values in Watts per meter-Kelvin (W/m-K) except number of samples. Figure 26. (a) Histogram and cumulative distribution function for measured thermal conductivity data. (b) Scatterplot of thermal conductivity as a function of matrix porosity. All thermal conductivities measured at 70°C and 105°C-dried conditions. Solid line in (b) is regression fit (used in Section 6.4.5.3); dashed lines are 95-percent confidence limits. Matrix Porosity, as a fraction 0.0 0.2 0.4 0.6 Thermal Conductivity, in W/m-K 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Kth = 1.5544 -2.4910 f r2 = 0.586 r = -0.766 ThKporosXplot.wmf (b) Matrix Thermal Conductivity, in W/m-K 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Frequency 0.00 0.05 0.10 0.15 0.20 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 ThKhisto.wmf (a) Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 74 of 150 6.4.6.2 Sequential Gaussian Simulation of Porosity Geostatistical simulation comprises a large class of modeling techniques that can produce very complex, and presumably therefore highly realistic numerical representations of spatially variable properties. Simulation may be thought of as “expanding” the actual information available in a stochastic manner that is also compatible with additional information derived from the data ensemble and the spatial context of those data. The process builds on the geologic intuition that Figure 27. Histograms and cumulative distribution functions for (a) total porosity measured for thermal conductivity test specimens and (b) lithophysal porosity from the TSw model unit (repeated from Figure 26(b)). Matrix Porosity, as a fraction 0.0 0.2 0.4 0.6 Frequency 0.00 0.05 0.10 0.15 0.20 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 ThKPorosHist.wmf Porosity, as a fraction 0.0 0.2 0.4 0.6 Frequency 0.00 0.05 0.10 0.15 0.20 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 TSwL_RHpor.wmf (b) (a) Figure 28. Histogram and cumulative distribution function for thermal conductivity systematically predicted from lithophysal porosity values in drill holes USW SD-7, SD-9, and SD-12 using the regression relationship from Figure 26(b). Thermal Conductivity, in W / m-K 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Frequency 0.00 0.05 0.10 0.15 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 predThKhist.wmf Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 75 of 150 unsampled locations nearby a known value will “tend” to resemble that value, whereas unsampled locations at increasing distances from a known measurement tend progressively to resemble that datum less and less. This intuition will be observed statistically across a suite of several equiprobable simulations. The philosophical framework of simulation is simple. Using concepts of random variables, one develops a model of the probability density function (pdf) for a material property of interest at all locations in space. By transforming the measured data to their respective positions on the probability density function and using simple kriging (Deutsch and Journel 1992, pp. 62 and 137), the desired pdfs can be made conditional to a set of measured values. Alternative realizations are simply generated by sampling from these pdfs. The variance of individual, location-specific, pdfs will vary with the amount of geologic uncertainty. Near conditioning data (Figure 30(c)), the pdf associated with an unsampled location will be relatively narrow. Where less information is known, such as away from data or in the vicinity of conflicting measurements, the pdf will be relatively broad (Figure 30(a-b)), leading to generation of a wide range of likely values across a suite of realizations. Because the underlying kriging algorithm used to derive the pdfs is an exact interpolator, the pdf degenerates to a spike with probability = 1 at a measured location (Figure 30(d)). Simulations may be conditional or unconditional. Conditional simulations are numerically anchored to a specific set of real-world data, and they exhibit three properties that adds to their usefulness in evaluating the effects of geologic uncertainty on physical process models. Specifically, conditional simulations: 1. Reproduce the known data values at the same locations within the model as represented by the real-world samples Table 20. Statistical Comparison of Measured and Predicted Thermal Conductivity Data for the TSw Model Unit (see text for discussion of prediction methodology) Measured Thermal Conductivity @ 70°C Predicted Thermal Conductivity at 70°C Mean 1.241 1.183 Std.Dev. 0.4271 0.1822 Minimum 0.620 0.676 Maximum 2.200 1.550 N 35 6063 Notes: Units are Watts/meter-Kelvin (W/m-K) except for number of values (N). 1. Also includes effect of measurement errors and lithologic variability. 2. Includes effect of lithologic variability only. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 76 of 150 2. Reproduce the full range of measurement variability, as represented by histogram and univariate descriptive statistics of the known data values 3. Reproduce the bivariate statistics, or two-point spatial correlation structure, of the input data. Unconditional simulations are similar, except that they are not conditioned to any particular spatially anchored data, and thus item 1 does not apply. As simulations with these three characteristics cannot be distinguished statistically from the ensemble of data used in their construction nor from each other, they serve as alternative, equally likely stochastic realizations of an incompletely sampled and measured reality shown conceptually in Figure 3. Simulations may be developed using parametric or nonparametric techniques for mechanically inducing the desired univariate (item 2 above) and bivariate (item 3) statistical properties. Parametric techniques rely upon the predictive power of well-understood multivariate probability Figure 29. Downhole variation in predicted thermal conductivity values based on systematically measured lithophysal porosity data for drill holes (a) USW SD-7, (b) SD-9, and (c) SD-12. Thermal conductivity key: predicted thermal conductivity at 70°C and 105°C-dried conditions. Porosity key: dark line—lithophysal porosity from petrophysical logs; light grey line with symbols—matrix porosity from core samples. USW SD-7 Po ros ity, frac tion 0.0 0.2 0.4 0.6 300 400 500 600 700 800 900 10 00 11 00 12 00 13 00 14 00 Th erm al Condu c tiv ity , W /m -K 0 1 2 Depth , feet (a) sd7predThK.wmf (b) d9 dThK f sd9predThK.wmf USW SD-12 Po rosity , fra ctio n 0.0 0.2 0 .4 0 .6 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 Thermal Conductiv ity, W /m -K 0 1 2 De pth , fee t (c) sd12predThK.wmf Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 77 of 150 functions, almost invariably the multivariate gaussian. A number of algorithms have been developed that implement gaussian-related simulation (for example, references in Deutsch and Journel, 1992). The sequential gaussian simulation program SGSIM (Table 1; see also Deutsch and Journel 1992, pp. 123 to 125 and 164 to 167) was used to generate 50 replicate models of porosity at all unsampled locations for each of the four model units, conditioned to the observed porosity data Table 21. Geostatistical Modeling Grid Specification Parameters Corresponding to the Southwest Corner of the Grid Grid Dimension Midpoint (ft/m) Spacing (ft/m) No. of Nodes Total Nodes Model X (Easting) 549,895.6 167,600.0 656.2 200.0 34 -- Model Y (Northing) 751,349.0 229,000.0 656.2 200.0 43 -- Model Z (Stratigraphic Vertical) PTn 3.3333 1.0159 6.667 2.032 30 43,860 TSw 8.2025 2.5000 16.405 5.000 61 89,182 CHn 8.4444 2.5399 16.666 5.080 24 35,088 Tcp 8.4444 2.5399 16.666 5.080 24 35,088 Probability Property Value (a) (b) (c) (d) Figure 30. Conceptual probability density functions representing the uncertainty associated with various unsampled locations. (a) Beyond the range of spatial correlation: pdf is virtually identical to the univariate histogram; essentially all that is known about the unsampled location is what is known about the population as a whole. (b) Far away from a sample, but within the range of spatial correlation: pdf is broad, indicating considerable uncertainty; distribution begins to focus on expected value. (c) Nearby a sample value: pdf is narrower indicating less uncertainty. (d) Immediately adjacent to a sample value: pdf is nearly a spike value corresponding to the adjacent sample datum. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 78 of 150 from the several drill holes. The sequential modeling process is relatively straightforward and is implemented as follows. 1. All data values are converted to positions on a univariate standard-normal (µ=0, s2=1) distribution using the graphical normal-score transform (Figure 31), implemented in program NSCORE (Table 1; see also Deutsch and Journel 1992, pp. 138 and 209 to 211). This transformation does nothing to the spatial correlation structure because the relative positions of all values with respect to each other are preserved (i.e., the transform is quantile-preserving). 2. The spatial correlation structure is identified using the normal-score transformed values and modeled using standard variography. 3. The transformed measured data are mapped into the model volume; samples located (only fortuitously) at a node in the stratigraphic-coordinate grid are assigned to that node and the node is not simulated. 4. A sequential random path is defined that will visit each unsampled node once and only once. 5. At each node along this path, a search is conducted for “nearby” data and any previously simulated grid nodes. The search parameters (anisotropic radii; number of data to use) are user specified. 6. The (user-specified) N closest data are identified and weighted by their “geological distance” (in contrast to simply their Euclidean distance), as defined in the stratigraphic coordinate system according to the mathematical formulation of the spatial continuity model (variogram). Because the normal-score transformed values are effectively relative positions on a cumulative distribution function, the resultant value is also effectively a relative position on the same cumulative distribution function. 0.0 0.2 0.4 0.6 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 -3 -2 -1 0 1 2 3 (b) Original Variable Transformed Variable (a) Figure 31. Graphical representation of the quantile-preserving normal-score transform process using cumulative distribution functions. A population with virtually any univariate distribution (a) can be transformed to any other univariate distribution (b) (here standard gaussian) in a manner represented by the arrows such that the quantile relationships among the data are preserved. The reverse transformation is also possible in the same manner. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 79 of 150 7. A value (in normal-score space) is drawn at random from the conditional probability distribution defined in step 6 and this value is assigned to represent the porosity at that point. The simulation process then moves to the next unsampled location along the random path defined in step 4 and the process is repeated beginning with step 5. 8. After all originally unsampled grid nodes have been simulated using the logic of steps 5 through 7, the resulting spatial array of normal-score values are back-transformed to the original porosity space using the inverse of the normal score transform of step 1 (software routine BACKTR, Table 2), and the simulation process is complete. Because porosity values are drawn at random for each unsampled grid node, the values obtained in different simulation runs will be different. Indeed, the weighting scheme used to develop the conditional expectation in each independent simulation will be different as well in that the same path through the 3-D grid is not used in successive simulations. Additionally, because the datasearch process considers previously simulated grid nodes as well as measured data (non-varying), the nearby values used to estimate the conditional expectation will also vary among simulation runs. At grid locations that are well constrained by consistent measured data, the variability of the simulated values across a suite of simulations will be small, as described by the spatial continuity model. However, at grid locations far from any conditioning measured data, or at grid nodes that are in the vicinity of conflicting measurements, the spread of porosity values that will be generated by the simulation algorithm across different computer runs will be quite broad, approaching the univariate variance of the data when considered without regard for spatial position. Uncertainty, measured by variability across the suite of simulations is small where much is known about the rock mass and progressively greater at longer distances from actual sampled values. 6.4.6.3 Indicator Kriging Using Uncertain Data The indicator methodology is most easily understood through the following definition of an indicator transformation of a spatially distributed variable, Z: , (Eq. 3) where I(x) is the indicator-transformed variable at spatial location x, Z(x) is the original variable also at location x, and Z* is a particular threshold value of interest. In other words, the original variable is replaced by the value 1 at all locations where it is less than or equal to the threshold and by the value 0 at all locations where it is greater. All of this presumes that there is precise knowledge of the original variable of interest, Z(x). Should knowledge of Z be less than 100-percent certain, however, it is possible to restate Equation 3 to account for the “soft” probability, p, that Z less than or equal to Z*. Thus: . (Eq. 4) I x ( ) 1 Z x ( ) Z* = . 0 Otherwise . . . . = I x ( ) p Z x ( ) Z* = . 1 p – Otherwise . . . . = Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 80 of 150 Note that it is perfectly legitimate to mix hard and soft indicators in the same analysis, and thus the methodology is suitable for combining the relatively definitive XRD indicators of hydrousphase mineral alteration with the less-certain-but-spatially-more-abundant petrophysical indicators (Sections 4.1.4 and 4.1.5, respectively). It remains to calibrate the appropriate value(s) of p, given the petrophysical estimates of bound-water content; this calibration effort is described below in Section 6.4.7.1. Implementation of indicator kriging process (program IK3D, Table 1; see also Deutsch and Journel, 1998, pp. 103 to 106) is straightforward enough in concept. 1. Hard XRD measurements of hydrous-phase mineral content are converted to a set of ones and zeros with respect to some threshold concentration judged relevant. 2. Soft petrophysical measurements of hydrous-phase mineral content are converted to a set of values ranging between zero and one in accordance with the calibration exercise (using program BICALIB, Deutsch and Journel 1998, pp. 235 to 236). 3. The spatial continuity structure of the composite hard and soft data set is determined through standard variography. 4. The transformed measured data are mapped into the model coordinate system. 5. At each desired grid location a search is conducted for nearby indicator data, I(x) (either hard or soft). 6. The nearest N values identified by the search are then weighted according to their relative geologic distances, as computed from the variogram model. The weighted average is assigned as the appropriate indicator value at that location, and the estimation process moves to the next grid node. Steps 5 and 6 are iterated until the grid is complete. Note that implicit in Equations 3 and 4 is the probabilistic interpretation that p is the likelihood of Z(x) being less than the threshold, Z*. If one is dealing with a binary classification in which the rock present at location x is either “altered” versus “unaltered,” it is a logical conclusion that: 7. All grid nodes with a modeled probability of p < 0.5 may be considered unaltered, and that those with p > 0.5 may be deemed altered. 6.4.6.4 Coregionalization Modeling of Derivative Properties There are two frequently used methods that have been used to incorporate cross-variable correlations of material properties into rock properties models in the presence of undersampling of one variable. First, one may assume a coefficient of determination (r2) equal to one and simply apply the empirically determined regression equation to predict the secondary (undersampled) variable. A second method is to model a randomly distributed error (“noise”) about the regression line. However, neither of these two alternatives is particularly satisfying. In many instances, r2 = 1.0 implies a substantially stronger relationship than exists in fact. In the second case, a cross-plot of the two resulting variables reproduces the desired cross-variable relationship, but any spatial correlation exhibited in nature by the secondary variable is effectively destroyed by the addition of spatially uncorrelated noise. Cokriging and cosimulation (Journel and Huijbregts 1978, pp. 324 to 326; Deutsch and Journel 1992, pp. 121 to 123) are wellestablished algorithms for producing models that reproduce both the observed correlation Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 81 of 150 between two variables and the observed spatial correlation structure of each. However, in the presence of severe undersampling of the secondary variable, such as is the case at Yucca Mountain, there are simply insufficient data from which to infer the necessary auto- and crosscorrelation relationships. A practical alternative mechanism to cokriging or cosimulation in the presence of undersampling involves the assumption of a linear model of coregionalization (Journel and Huijbregts 1978, pp.171 to 175 and 516 to 517; Luster 1986; Altman et al. 1996, pp. 54 to 59), which effectively states that the spatial continuity of the secondary variable, and of the joint (cross-variable) correlation structure, is presumed to be approximately identical to that of the primary variable. Although the mathematical relationships involved in the derivation are tedious (e.g. Rautman and McKenna 1997, pp. 63 to 66, and references therein), it can be shown that it is possible to generate a pair of cross-correlated simulations as a weighted linear combination of two independent standard-normal spatially correlated simulations where the weighting factor is effectively the cross-variable correlation coefficient, r. What this translates to in practical terms, is to take a conditional porosity simulation in standard-normal form, such as is generated in step 7 of Section 6.4.6.2, and combine it (software routine COREGPC, Table 2) with an identically generated but independent unconditional simulation (these are paired in this activity by run number), also in standard-normal form, using r and as the appropriate weighting coefficients. The resulting linear combination is also essentially in standard normal form, and it can be back-transformed to match the desired distribution of the target secondary variable (software routine TRANS, Table 2). If appropriate, the back transformed coregionalized property simulation is post-processed further to impart additional relevant characteristics as described below in Section 6.4.8. 6.4.7 Modeling of Hydrous-Phase Mineral Alteration Volcanic glass within the lower two modeling units (CHn and Tcp) has been variably altered to (dominantly) zeolite minerals throughout a major portion of the model area. These altered rocks exhibit markedly reduced saturated hydraulic conductivity by comparison with unaltered materials of approximately the same porosity. The rock properties modeling effort attempted to include not only the XRD mineralogic data (described in Section 4.1.4), which provide virtually 100-percent certain identification of hydrous-phase mineral alteration, but also less accurate but more abundant and widely distributed petrophysical indicators of such alteration (Section 4.1.5). 6.4.7.1Calib rating Soft Ind icators of Hyd rous-Phase M ineral A lteration Because the petrophysical data available provide a less-than-100-percent certain identification of hydrous-phase alteration, it is necessary to calibrate these values to account for the added uncertainty. The calibration effort involved 334 samples from the CHn and Tcp model units for which depth-matched pairs (generated using software routine MATCHUP, Table 2) of both XRD mineral analyses and petrophysical bound-water contents could be obtained. Petrophysical bound-water content was calculated simply as the difference between the POROTOT and POREF log traces (see also Section 4.1.5). Core bound-water content was calculated as the difference between the OD and RH porosity values (Section 4.1.1), and is initially identical to the deltaporosity value described in Section 6.4.5.1. Comparison of depth-matched core and petrophysical bound-water data (software routine MATCH, Table 2) indicate that the laboratory core measurement 1 r2 – Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 82 of 150 process gave a bound-water content of approximately twice that indicated by the down-hole petrophysical measurements; see Figure 32. Although a precise explanation for the discrepant measurements is uncertain, the empirical relationship can be used to adjust the core measurements to provide a common basis with the petrophysical values. Adjusted values for core bound water content have been used in the calibration work that follows. Figure 33 presents a scatterplot of total hydrous-phase mineral content vs. adjusted bound-water content. In general, an increase in adjusted bound-water content corresponds directly to an increase in the total hydrous minerals. The calibration consists of cross-tabulating (software routine BICALIB, Table 2) the number of pairs in each of the categories: • Hydrous-phase mineral content: – Greater than 5 percent – Less than or equal to 5 percent. • Adjusted bound-water content: – Less than or equal to 0.03 – Greater than 0.03 to less than 0.04 – Greater than 0.04 to less than 0.05 – Greater than 0.05. Figure 32. Scatterplot of core vs. petrophysically derived bound-water content for 354 depth-matched pairs of samples. Regression line (solid, with 95-percent confidence limits, dotted) has been forced through the origin and has a slope of 1.74; unconstrained regression has r2 = 0.67 (r-squared for a constrained line is not meaningful). Petrophysical Bound-water Fraction 0.00 0.02 0.04 0.06 0.08 0.10 Core Bound-water Fraction (delta f) 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 1:1 2:1 Petrophysical Bound-water Content Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 83 of 150 The cross-tabulated counts of soft-value pairs were converted to a decimal proportion (software routine BICALIB), and these values were taken as the prior probability of obtaining the specified total hydrous-phase mineral content given a specified adjusted bound-water content. These soft, prior-probability values are presented in Table 22. In contrast, XRD hydrous-phase mineral contents were coded as hard probability values of zero or one. 6.4.7.2 Indicator Kriging of Hydrous-Phase Mineral Alteration Both hard and soft indicators of hydrous-phase mineral alteration were combined and supplied as input to an indicator kriging exercise (program IK3D; Deutsch and Journel, 1998, pp. 103 to 106). Indicator kriging (Section 6.4.6.3) is a variant of ordinary kriging, in which the variable of interest is estimated as a weighted linear combination (average) of the available observed values within some local neighborhood of influence. In common with ordinary kriging, the weights applied to the observed values are calculated in accordance with the spatial continuity model developed from the combined indicator data set. Necessarily, alteration in the CHn and Tcp model units was Table 22. Prior Probability of Altered vs. Unaltered Rocks as a Function of Bound-Water Content Hydrous-Phase Mineral Content Bound-Water Content, Cumulative Prior Probability < 0.03 0.03 – 0.04 0.04 – 0.05 > 0.05 < 5 percent 0.2474 0.6122 .7778 0.8625 > 5 percent 0.7526 0.3878 0.2222 0.1375 Figure 33. Scatterplot of total hydrous-phase mineral content vs. adjusted bound-water content for 334 depth-matched paired samples. Interior lines correspond to the various threshold values used in Table 22. Petrophysical Bound Water Fraction 0.00 0.05 0.10 Total Hydrous Minerals, percent 0 20 40 60 80 100 Core Data Non-core Data Petrophysical Bound-water Content Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 84 of 150 modeled separately, even though the two units were combined for purposes of estimating the prior probability values discussed above. CHn Model Unit–The variogram model for hydrous-phase mineral alteration indicators in the CHn model unit is presented in Figure 34, and the parameters of the fitted variogram model are in Table 23. Details of the variogram modeling exercise may be found in the scientific notebook (Rautman and McKenna 1998, pp. 663 to 695). Table 23. Variogram Parameters for Spatial Continuity Model, Alteration in the CHn and Tcp Model Units Nest No. Model Type Range (ft) Sill Rotation Angle (degrees) Anisotropy Ratio Maximum (horizontal) Intermediate Minimum (vertical) 1 2 3 1 2 CHn Model Unit -- Nugget -- -- -- 0.005 -- -- -- -- -- 1 Spherical 4000 1500 30 0.010 0 0 0 0.3750 0.0075 2 Spherical 7000 4000 500 0.035 0 0 0 0.5714 0.0714 Tcp Model Unit -- Nugget -- -- -- 0.010 -- -- -- -- -- 1 Spherical 2500 2500 150 0.025 0 0 0 1 0.060 2 Spherical 15000 15000 150 0.058 0 0 0 1 0.010 Figure 34. Indicator variogram and fitted model computed for alteration in the CHn model units. (a) Stratigraphically vertical and (b) stratigraphically horizontal directions. Separation Distance, h, in feet 0 5000 10000 15000 Variogram 0.00 0.02 0.04 0.06 0.08 0.10 Data N-S Data E-W a priori Variance Model N-S Model E-W CHnHorizAltVario.wmf (b) Separation Distance, in feet 0 50 100 150 200 250 300 350 Variogram 0.00 0.02 0.04 0.06 0.08 0.10 Data a priori Variance Model CHnVertAltVario.wmf (a) Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 85 of 150 Tcp Model Unit–The variogram model for hydrous-phase mineral alteration indicators in the Tcp model unit is presented in Figure 35. The parameters of the fitted variogram model are given in Table 23. Details of the variogram modeling exercise may be found in the scientific notebook (Rautman and McKenna 1998, pp. 772 to 786). 6.4.8 Postprocessing of Simulated Models 6.4.8.1In corporation of S pecific Attributes into Si mulated M odels The actual rocks at Yucca Mountain are the composite result of numerous geologic processes that overlap in both space and time. Consequently, rock properties modeling involves more than the simple generation of a set of porosity values. This is particularly true for the models of derivative material properties that have been generated by coregionalization with porosity. This section presents the techniques used to incorporate two specific types of secondary geologic attributes into raw simulated property models. Vitrophyres–The widely variable hydraulic conductivity values associates with densely welded vitrophyric core samples of uniformly low porosity from the Topopah Spring Tuff (lithostratigraphic units Tptrv1 and Tptpv3) have been described by Flint (1998, p. 38) as resulting from microfractures present within these glassy, brittle rocks. Additional consideration of these samples (Flint 1998, Figure 12) suggests that vitrophyric samples may independently be identified by their uniformly very low porosity (less than approximately 0.05; see Section 6.4.5.1, Figure 22). Accordingly, the simulated models of saturated hydraulic conductivity for the TSw model unit (only) were post-processed (software routine VITROPHYRE, Table 2) such that if the corresponding porosity value was less than 0.05 (simulations paired by run number), the coregionalized hydraulic conductivity value was discarded. Under the assumption that such low-porosity grid nodes represent vitrophyre or other essentially nonporous brittle materials, a value of Ks was generated by random sampling from a uniform population with a range of 10-14 to Figure 35. Indicator variogram and fitted model computed for alteration in the Tcp model unit. (a) Stratigraphically vertical and (b) both stratigraphically horizontal directions (isotropic). Separation Distance, in feet 0 100 200 300 Variogram 0.00 0.05 0.10 0.15 TcpVertAltVario.wmf (a) Separation Distance, in feet 0 5000 10000 15000 Variogram 0.00 0.05 0.10 0.15 TcpHoirzAltVario.wmf (b) Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 86 of 150 10-6 m/sec as derived from the histogram of Figure 22. A conceptual representation of the logic underlying the modeling of vitrophyric rocks is presented in Figure 36. Hydrous-Phase Mineral Alteration–Hydrous-phase mineral alteration is inferred to represent a secondary alteration process that affected vitric tuffaceous materials at some particular time after formation of the original rock mass, although generally before tectonic faulting and tilting. Consequently, there appears to be little or no direct correlation of saturated hydraulic conductivity with matrix porosity for altered (zeolitized) samples, see Section 6.4.5.1, Figure 21). Recall also from Figure 21 that the histogram of altered hydraulic conductivity values appears virtually indistinguishable from a gaussian population, given the relatively small sample size. This modeling philosophy has been implemented for the rock properties model by postprocessing (software routine ZEOLITE5, Table 2) the initial coregionalized hydraulic conductivity models (for both the CHn and Tcp model units) grid node by grid node together with a corresponding indicator kriging model (Section 6.4.7.2) indicating the probability of significant hydrous-phase mineral alteration. If the grid node under consideration was considered unaltered (palt < 0.50), the coregionalized Ks value was retained, and the processing moved to the next grid node. If the node was considered altered (palt > 0.50), then the coregionalized Ks value was discarded in favor of a Figure 36. Logic diagram for post-processing porosity and hydraulic conductivity simulations to recognize vitrophyre rock type. See text for discussion. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 87 of 150 normally distributed random value sampled from a population with the appropriate mean and variance. A schematic diagram of this postprocessing procedure is presented in Figure 37. 6.4.8.2 Uncertainty Modeling Geostatistical generation of the material property models is essentially complete at this point. As part of the current modeling exercise, 50 replicate, statistically indistinguishable models of porosity for each model unit (one set each for matrix and lithophysal porosity in the TSw model unit) and 50 replicate models for each one of the derivative properties (bulk density, matrix saturated hydraulic conductivity, and—for the TSw model unit—thermal conductivity as well) have been generated. Each of the replicate simulations honors the measured porosity data at sample locations (subject to the discretization limits), exhibits the full range of variability captured by the histogram of the relevant property, exhibits the appropriate range of spatial correlation (variogram), and (for the derivative properties) exhibits the appropriate correlation coefficient with porosity. In effect, there is nothing objective about any one simulated (or coregionalized) model to prefer it over any other model of that suite. Indeed, the only meaningful distinguishing feature within a suite of replicate models is the arbitrarily selected random number seed which initiated the simulation process. Figure 37. Logic diagram for post-processing porosity and alteration indicator simulations to recognize hydraulic conductivity dependence on alteration state. See text for discussion. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 88 of 150 Because there are few, if any, objective differences to distinguish the members of each suite of simulated property models, it thus logically follows that the variability among members of a suite represents an empirical estimate of the geologic uncertainty associated with each material property. “Geologic uncertainty” is defined in this context as that which results from less-thanexhaustive sampling or other measurement. The difficulty arises, however, as how best to represent this space of uncertainty in a simple and concise manner. An “uncertainty model” has been generated (software routine ETYPE; Table 2) for each material property-modeling unit combination by computing the node-by-node standard deviations of each set of 50 replicate simulated models. This process produces uncertainty models that are themselves spatially heterogeneous. By theory and in practice, variability among simulations— and uncertainty as defined by the standard deviation—is small in close proximity to measured sample values. Variability among simulations and uncertainty is high at great distances from measured data, or in the vicinity of conflicting measured values. Note that there is no particular reason other than general reader familiarity for selecting the standard deviation as “the” measure of uncertainty. Values such as the total range of the modeled property or the interquartile range could easily be computed during this post-processing step. With respect to alternative uncertainty models, it is also important to remember that the “best” measure of geologic uncertainty for the potential nuclear-waste repository at Yucca Mountain is the impact of that uncertainty on some relevant measure of repository system performance. Potential examples of such global performance measures might be particle-tracking travel times or cumulative radionuclide release rates, as suggested by the consequence analysis step shown in the conceptual diagram of Figure 3. Development of these types of comprehensive uncertainty assessments is beyond the scope of the rock properties modeling effort. 6.4.8.3 “Expected-Value” Modeling The post-processing step used to generate the standard-deviation uncertainty model for each rock property-modeling unit combination described in Section 6.4.8.2 is easily adapted to produce some type of summary statistical measure associated with the “central-tendency” of the simulation process. In effect, the geostatistical simulation process develops a 3-D array of spatially correlated probability density functions, one distribution for each of the discretizing grid nodes within a model unit. Continuing with the probability density function analogy, certain rock property values are more likely to exist at a given location than other values. Note, however, that “more likely to exist” can be defined in more than one manner. A set of summary “expected-value” models has been generated (software routine ETYPE; Table 2) for each suite of simulated models by calculating the arithmetic mean of the 50 replicate simulated values generated at each grid node. Journel (1983, pp. 459 to 461; see also Deutsch and Journel 1992, pp. 76 and 225) defined this conditional expectation representation of central tendency as the “E-type estimate.” Note that computing the median (“M-type estimate,” Journel 1983, p. 460) of the 50 replicate simulated values would be an equally legitimate measure of “expectation,” although one that is somewhat more computationally intensive. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 89 of 150 Because of the logistical difficulty of presenting the full simulated results for 50 models times 19 unique material-property/model-unit combinations, the results of this geostatistical modeling exercise described in Section 6.5 that follows are presented in terms of the E-type estimate. It is important to remember, however, that the three characteristic properties of simulated models described as desirable in the first paragraph of Section 6.4.6.2 no longer necessarily apply to these summary models. Compare, for example, the statistical nature of the simulated and summary models from ISM2, as discussed by Rautman and McKenna (1997, e.g., pp. 99 to 104), and in abbreviated form for ISM3 as presented in the discussion of model validation in Section 6.7, below. The first characteristic, that involving reproduction of the measured (porosity) values at the locations of actual measurement, is maintained. However, the ensemble of modeled E-type values no longer represents the full range of univariate variability of the measurement ensemble (the second characteristic). Additionally, the two-point spatial correlation character (variogram) of the E-type model no longer reproduces that of the underlying measurements (the third characteristic). Specifically, because of averaging across the replicate simulations, the E-type model typically grades relatively smoothly and continuously from one (exactly reproduced) measured value to the next (in three dimensions). Thus the apparent spatial continuity of the Etype model typically is much greater than that observed for the data themselves. This is the socalled smoothing effect that is typical of virtually all interpolation (in contrast to simulation) algorithms, including kriging, nearest-neighbor estimation, and inverse-distance-to-a-power weighting. 6.5 RESULTS AND DISCUSSION The “real” results of the rock properties modeling exercise consist of the numerical model files containing the replicate Monte Carlo simulated models and the associated summary E-type and uncertainty model files for each of the four modeling units and each of the modeled properties. Because of the impracticality of presenting and discussing each of these simulated and summary models individually, the results presented in this section are “illustrative” and focus on gross heterogeneity features as revealed by the summary E-type models. Each model unit and each modeled material property within that model unit is discussed in turn. Discussion of the replicate individual simulated models is generally restricted to the section on Model Validation. Again, because of the impracticality of “validating” each and every replicate simulation, the discussion in this section will be illustrative. Because of technical limitations involving plotting the individual drill holes on some model views, refer to Figures 4 though 7 for information regarding the spatial distribution of input measurements. Note that the distribution of available drill hole information is unit-specific. Also note that although the color-scale legend on each figure generally represents the full range of the modeled property present in the underlying numerical file, it is only fortuitous if the highest (or lowest) value is actually shown in the particular model view. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 90 of 150 6.5.1 PTn Model Unit Heterogeneity of matrix porosity within the PTn model unit is shown in Figure 38 in both stratigraphic and real-world coordinates. As indicated by the projection arrows connecting the two halves of the figure, the vertically exaggerated rectangular volume in stratigraphic coordinates is back transformed to real-world coordinates, such that the material property values assume their correct relative positions within the tilted and faulted strata of Yucca Mountain, shown with no vertical exaggeration. Cross-sectional views of porosity heterogeneity are presented in Figure 39. Porosity values within the PTn model unit are generally high, varying principally from about 30 to more than 60 percent (0.30 to 0.60). Porosity values appear relatively continuous over distances of 5,000 to 10,000 feet, as expected from the input range of spatial continuity. Porosity trends are prominently anisotropic from northwest to southeast, also as expected from the variogram model (top surface of the block diagram of Figure 38). Bulk density (Figure 40) varies spatially approximately inversely with porosity, as expected from the strong (r = –0.912) negative correlation coefficient. Bulk density values vary from less than 1.0 g/cm3 to nearly 2.0 g/cm3. Densities are generally low across the modeled area. Prominent regions of higher density are associated with drill hole H-6, shown on the southernmost east-west cross section, and particularly with drill hole G-2 near the northern boundary of the modeled region at the intersection of the north-south and northernmost east-west cross sections. At the G-2 locality, density exceeds 1.9 g/cm3 at two horizons, presumably corresponding to the Pah Canyon and Yucca Mountain Tuffs, which are described as moderately welded in this part of the model area. Heterogeneity in matrix saturated hydraulic conductivity values is shown in Figure 41. Hydraulic conductivity values are generally between 10–5 and 10–7 m/sec (note that in this and following similar figures, hydraulic conductivity values are shown in log10 units; i.e., 10-7 m/sec = –7.000). Lower conductivities on the order of 10–8 m/sec are modeled in the vicinities of drill hole H-6 (not visible in figure) and G-2, coincident with the lower matrix porosities in these regions from which the hydraulic conductivity values are coregionalized. 6.5.2 TSw Model Unit Heterogeneity of material properties within the TSw model unit is presented in Figures 42 to 48. Note that matrix porosity, as indicated in Figure 42 and 43, is very low, mostly less than 10 to 15 percent, and relatively constant in magnitude across the entire modeled region. This minimal variability is consistent with the definition of this unit as densely welded tuff. However, the front face of the block diagram, shown in stratigraphic coordinates in Figure 42, and the cross-section views of Figure 43 indicate that increased matrix porosity is associated with the major lithophysae-bearing intervals within the unit, and most particularly with what in the GFM would be the vapor-phase corroded crystal-rich nonlithophysal interval near the top of the unit. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 91 of 150 Figure 38. Perspective diagrams showing E-type model matrix porosity in the PTn model unit in both stratigraphic and real-world coordinates. Light-grey objects in real-world-coordinate view are drill holes and workings of the ESF. Easting and northing values are Nevada state plane coordinates in feet. ESF Nevada State Plane Coordinates Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 92 of 150 Figure 39. Cross-sectional views showing E-type heterogeneity in matrix porosity in the PTn model unit in stratigraphic coordinates. Vertical exaggeration 10x. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 93 of 150 Figure 40. Cross-sectional view showing E-type heterogeneity of bulk density in the PTn model unit in stratigraphic coordinates. Vertical exaggeration 10x. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 94 of 150 Figure 41. Cross-sectional views showing E-type heterogeneity in matrix saturated hydraulic conductivity in the PTn model unit in stratigraphic coordinates. Vertical exaggeration 10x. Approximate Location of G-2 Approximate Location of H-6 Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 95 of 150 ESF Figure 42. Perspective diagrams showing E-type model matrix porosity in the TSw model unit in both stratigraphic and real-world coordinates. Light-grey objects in real-world-coordinate view are drill holes and workings of the ESF. Easting and northing values are Nevada state plane coordinates in feet. Note compressed porosity scale compared with Figure 44. Nevada State Plane Coordinates Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 96 of 150 Figure 43. Cross-sectional views showing E-type heterogeneity of matrix porosity in the TSw model unit in stratigraphic coordinates. Vertical exaggeration 10x. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 97 of 150 Figure 44. Perspective diagrams showing E-type model lithophysal porosity in the TSw model unit in both stratigraphic and real-world coordinates. Light-grey objects in real-world-coordinate view are drill holes and workings of the ESF. Easting and northing values are Nevada state plane coordinates in feet. Note expanded porosity scale compared with Figure 42. ESF Nevada State Plane Coordinates Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 98 of 150 Figure 45. Cross-sectional views showing E-type heterogeneity of lithophysal porosity in the TSw model unit in stratigraphic coordinates. Vertical exaggeration 5x. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 99 of 150 Figure 46. Cross-sectional views showing E-type heterogeneity of bulk density in the TSw model unit in stratigraphic coordinates. Vertical exaggeration 5x. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 100 of 150 Figure 47. Cross-sectional views showing E-type heterogeneity of thermal conductivity in the TSw model unit in stratigraphic coordinates. Vertical exaggeration 5x. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 101 of 150 Figure 48. Cross-sectional views showing E-type heterogeneity of matrix saturated hydraulic conductivity in the TSw model unit in stratigraphic coordinates. Vertical exaggeration 5x. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 102 of 150 The heterogeneity of lithophysal porosity, in contrast to the matrix porosity, is shown in Figure 44 as perspective views, and in Figure 45 in cross-sectional view. Although even lithophysal porosity is generally low compared to the porosity of a nonwelded welded tuff, such as the PTn model unit, maximum porosity values within the lithophysal intervals locally exceed 30 to 35 percent (0.30 to 0.35). Figure 45 clearly indicates two such intervals of high porosity, corresponding approximately to the upper and lower lithostratigraphic units (Tptpul and Tptpll), separated by a low-porosity (on the order of 10 percent) interval equivalent to the middle nonlithophysal lithostratigraphic unit (Tptpmn). Note however, that there is a fairly large amount of lateral heterogeneity within each of the elevated porosity intervals. It is the measured porosity data themselves, as propagated away from drill hole locations by the spatial continuity model, that produce both the apparent layering and the variations within those layers. There are no detailed lithostratigraphic (or other) subunits explicitly modeled within the TSw model unit. All property heterogeneity in the rock properties model are functions strictly of the measured material properties. Bulk density heterogeneity, which is coregionalized from lithophysal porosity, is illustrated in Figure 46. Density values are typically above 2.0 g/cm3 throughout most of the relatively lithophysae-free region. Bulk density is particularly high (approaching 2.5 g/cm3) in the lower parts of the TSw model unit, as indicated by the red colors on the figure. A prominent high density interval is associated with the lower vitrophyre in the central part of the modeled region (approximately corresponding to lithostratigraphic unit Tptpv3). However, bulk-rock density values associated with the upper lithophysal horizon in particular may be as low as 1.5 to 1.8 g/ cm3. The alternation of lithophysal and nonlithophysal intervals are particularly clearly visualized through the bulk density model Thermal conductivity is also coregionalized from lithophysal porosity, in an effort to predict the thermal conductivity of volumes of rock influenced by the presence of lithophysal cavities one decimeter or larger in diameter. Figure 47 presents the E-type model of spatial heterogeneity in thermal conductivity. High values of thermal conductivity (greater than approximately 1.3 to 1.4 W/m-K and shown in yellow and orange tones) are associated with the lower portion of the TSw model unit, and with the presumed low-lithophysae “middle nonlithophysal” lithostratigraphic unit (Tptpmn). Particularly high thermal conductivity values, approaching 1.5 W/m-K, are present at the very base of the TSw unit, presumably associated with the densely welded vitric lithostratigraphic unit. In contrast, the most lithophysal portions of the unit, which contain lithophysal cavities up to a meter in diameter, appear to be characterized by bulk-rock thermal conductivities less than 1.0 W/m-K (blue and blue-green colors in Figure 47). Heterogeneity in matrix saturated hydraulic conductivity is presented in Figure 48, and has been coregionalized from the matrix porosity model shown in Figures 42 and 43. Matrix conductivities are less than 10–11 m/sec through much of the lower part of the model unit. As expected from the higher matrix porosity values associated with the lithophysae-bearing portions of the TSw model unit, matrix hydraulic conductivity values are markedly higher, 10–9 to 10–10 m/sec (yellow to green tones), in these vapor-phase altered portions of the unit. Some of these higher conductivity values may also be associated with vapor-phase corrosion of the welded tuff within the crystal-rich nonlithophysal unit. Note, however, that “matrix” hydraulic conductivity Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 103 of 150 does not include conductivity attributable to flow through lithophysal cavities or fractures under true saturated conditions. 6.5.3 CHn Model Unit Variations in the matrix porosity of the CHn model unit are presented in Figures 49 and 50. Porosity values are generally high at 20 to 40 percent (0.20 to 0.40) throughout the unit, particularly in contrast to the low porosity values typical of the overlying TSw model unit (0.10 to 0.15). Expanding the porosity color scale indicates that a mass of particularly high porosity occupies the central portion of the modeled volume. Porosity values locally approach 50 percent (0.50) within this region. Variations in bulk density in the CHn model unit are presented in Figure 51. Density values vary from more than 2.0 g/cm3 to less than 1.3 g/cm3, depending upon location, although the majority of this nonwelded unit exhibits limited variation in bulk density at 1.5 to 1.75 g/cm3. Generally speaking, bulk density varies inversely with porosity, as anticipated from the coregionalization relationship. Heterogeneity in matrix saturated hydraulic conductivity for the CHn model unit is illustrated in Figure 52. Although hydraulic conductivity is derived by coregionalization with matrix porosity, the relationship between the two material properties is not precisely straightforward because of the presence of hydrous-phase mineral alteration (predominantly zeolitic) within the unit. Matrix conductivities are typically 10–6 to 10–7 m/sec (greens to reds) within the unaltered portion of the CHn, and typically less than 10–11 m/sec elsewhere (blue). The block diagram in the upper part of Figure 52 clearly indicates that vitric (to potentially devitrified) materials are limited to the upper portion of the model unit, and more particularly to the southwestern portion of the modeled volume. The overall impression is of a wedge of vitric (unaltered) material tapering to a feather edge toward the northeast. Zeolitic (altered) rocks are shown in tones of blue underlying and replacing the green- through red-colored volume to the north. Recall that the saturated hydraulic conductivity of altered samples are essentially uncorrelated with matrix porosity. The lower portion of Figure 52 presents cross-sectional views of the hydraulic conductivity field within the model unit, and attempts to illustrate some of the complex interfingering relationships of altered and unaltered rock types. 6.5.4 Tcp Model Unit Figures 53 and 54 present the spatial heterogeneity of matrix porosity within the Tcp model unit. Because the Prow Pass Tuff is mostly nonwelded, the porosity values are typically high, varying from 20 to nearly 40 percent across large volumes of the model. Lower porosity values, typically less than 15 percent (0.15), are present along the northern boundary of the modeled volume, and as a poorly defined lobate mass (indicated by the 0.15 porosity isoshell), which may correspond to the “moderately welded” portion of this unit, generally low within the central-eastern part of the region. Variations in bulk density, shown in Figure 55, substantiate the variations in porosity described in the preceding paragraph. Densities well in excess of 2.1 g/cm3 are prominently displayed along Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 104 of 150 ESF Cross-Block Drift Nevada State Plane Coordinates Figure 49. Perspective diagrams showing E-type model matrix porosity in the CHn model unit in both stratigraphic and real-world-coordinates. Light-grey objects in real-world-coordinate view are drill holes and workings of the ESF. Easting and northing values are Nevada state plane coordinates in feet. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 105 of 150 Figure 50. Cross-sectional views showing E-type heterogeneity of matrix porosity in the CHn model unit in stratigraphic coordinates. Vertical exaggeration 10x. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 106 of 150 Figure 51. Cross-sectional views showing E-type heterogeneity of bulk density in the CHn model unit in stratigraphic coordinates. Vertical exaggeration 10x. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 107 of 150 Figure 52. Block and cross-sectional views showing E-type heterogeneity of matrix saturated hydraulic conductivity in the CHn model unit in stratigraphic coordinates. Vertical exaggeration 10x. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 108 of 150 ESF Cross-Block Drift Nevada State Plane Coordinates Figure 53. Perspective diagrams showing E-type model matrix porosity in the Tcp model unit in both stratigraphic and real-world coordinates. Light-grey objects in real-world-coordinate view are drill holes and workings of the ESF. Easting and northing values are Nevada state plane coordinates in feet. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 109 of 150 Figure 54. Cross-sectional views showing E-type heterogeneity of matrix porosity in the Tcp model unit in stratigraphic coordinates. Vertical exaggeration 10x. Note adjusted porosity scale in lower figure. Objects extending away from cross sections in lower figure are isoshells enclosing all regions of lowest porosity (less than approximately 0.15). Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 110 of 150 Figure 55. Cross-sectional views showing E-type heterogeneity of bulk density in the Tcp model unit in stratigraphic coordinates. Vertical exaggeration 10x. Approximate Location of G-2 Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 111 of 150 the northern boundary of the model in the vicinity of drill hole G-2. Higher densities on the order of 2.0 g/cm3 are also visible in the east-central portion of the block, corresponding to the low porosity lobe. Elsewhere across the modeled volume, bulk densities are more typically between 1.75 and 2.0 g/cm3 and are shown in green colors. Heterogeneity in matrix saturated hydraulic conductivity is presented in Figure 56. In a manner similar to the overlying CHn model unit, a bimodal distribution of conductivity values, corresponding to altered and unaltered rock types, is quite prominent. Lower hydraulic conductivity values, typically less than 10–10 m/sec are associated with regions affected by hydrous-phase mineral alteration. Markedly higher values of hydraulic conductivity, varying from 10–10 to 10–7 m/sec, are associated with the vitric-to-devitrified continuum of matrix porosity values in regions unaffected by alteration. The block diagram in the upper portion of Figure 56 indicates that the upper and lower margins of the Prow Pass Tuff are essentially completely altered. Reference to the cross-sectional views in the lower half of the figure indicates that the unaltered portions of the unit correspond to the devitrified interior core of the ash flow. 6.6 UNCERTAINTIES AND LIMITATIONS “Uncertainty,” in the context of a stochastic modeling analysis, assumes a specific meaning as a descriptor of one’s state of knowledge, in that a quantitative assessment of uncertainty is typically one of the explicit objectives of the analysis. This meaning is in contrast to another possible definition, which is expressed more precisely as involving limitation or doubt as to the accuracy or relevance of a result. This distinction is maintained in this report. This section first describes a number of limitations of both methodology and data that detract from the exactness or accuracy of the models generated by this analysis. The results of a stochastic uncertainty analysis is then presented, which attempts to quantify rigorously the space of geologic uncertainty that results from less-than-completely “exhaustive” site characterization. 6.6.1 Limitations There are a number of factors affecting this analysis that may best be described as limitations of the data or of the modeling process itself. These limitations include errors and biases in the sample data used in the analysis, the methodological use of porosity as a surrogate for other material properties, the combination of numerous lithostratigraphic units into the four major modeling units, and the effect of geologic departures from the assumptions inherent in the use of the stratigraphic coordinate system. 6.6.1.1 Errors and Biases in Sample Data Stochastic simulation is a statistical-probabilistic methodology, with reproduction of various target statistical measures an important part. As such, errors (“uncertainties”) incorporated into the statistical description of the rock mass unquestionably will be propagated through the simulation process into the output models. These errors are of two principal types: measurement error and systematic bias. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 112 of 150 Figure 56. Block diagram and cross-sectional views showing E-type heterogeneity of matrix saturated hydraulic conductivity in the Tcp model unit in stratigraphic coordinates. Vertical exaggeration 10x. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 113 of 150 Measurement error or analytical uncertainty in the input data is not addressed in these rock property models. Philosophically, the underlying assumption is that measurement uncertainty, such as would be captured by replicate measurements of a given material property on the same physical sample, is small compared to the lateral heterogeneity of the actual rock mass. It is also assumed that these measurement types of errors are essentially unbiased, such that some may be high whereas others may be correspondingly low. Statistically, this assumption may be stated that the errors are of mean zero with a small relative variance. Systematic biases, on the other hand, represent a more severe challenge for rock properties modeling using a statistical-probabilistic approach. Two examples of systematic bias within the data sets used for this modeling activity are known, and have been compensated for to some extent during creation of the individual simulated models. These biases are discussed in the following paragraphs. Measurement Sensitivity Limits–Measurement of the matrix saturated hydraulic conductivity values for core samples appears to have had a lower “detection” limit of roughly 10–12 to 10–11 m/sec. Samples whose hydraulic conductivity are lower than this lower limit are reported in the data set as “no flow.” Omitting these samples entirely would lead to an unrealistically high set of modeled hydraulic conductivity values. On the other hand, substituting the no-flow samples with an arbitrary low value prior to the simulation process would tend to skew the results towards that arbitrary low conductivity. Therefore, as described in Section 6.4.7, the effect of these no-flow samples has been simulated explicitly during post-processing by setting the appropriate fraction of values equal to the arbitrary value of 10–14 m/sec at randomly selected grid nodes within each model unit. This “fix” to approximate the non-negligible number of non-flowing laboratory samples assumes that there is no particular spatial correlation among these samples. Some of these samples, or additional immediately adjoining samples, more recently have been retested in a permeameter with a lower detection limit somewhat less that 10–13 m/sec, but these revised values for the non-flowing samples have not been incorporated into the current generation of rock properties models. Preferential Sampling Bias–Another, more pronounced example of preferential bias in the target sample population involves the whole-rock thermal conductivity modeling. As described originally by Rautman and McKenna (1997, pp. 40 to 43), the laboratory measurements of thermal conductivity are systematically biased by preferential sampling of more coherent, generally lower-porosity/higher thermal conductivity, core specimens. This bias was identified by differences in the histograms of porosity for those laboratory thermal-test specimens (Figure 27) versus the overall histogram of porosity for the Topopah Spring welded unit. An additional confounding influence with respect to thermal conductivity is the effect of larger-than-core-size lithophysal cavities on the thermal conductivity of the rock mass as a whole. An attempt was made to reduce the impact of this identified sampling bias for thermal conductivity by constructing an “unbiased” reference distribution (see Rautman and McKenna 1997, pp. 40 to 43) of thermal conductivity values for simulation using a porosity-weighted distribution of estimated thermal conductivities (Section 6.4.5.3), where the porosity values were obtained by relatively rigorous systematic sampling of the entire TSw model unit. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 114 of 150 Similar sampling bias also affects matrix hydraulic conductivity determinations, in that laboratory testing is skewed slightly toward measurements of the more conductive samples (low permeability samples take longer to run) in certain drill holes. This bias was not addressed explicitly in this modeling work as the total number of laboratory hydraulic conductivity determinations is quite large (400-plus) compared to thermal conductivity (~50 total; 35 for the TSw model unit), and several drill holes (especially UZ-16, SD-9) were sampled on a quite systematic basis for the hydraulic property (compare, for example, the locations and measured values of hydraulic conductivity specimens sampled from SD-9 (Rautman and Engstrom 1996b, Figure 11) with the corresponding information for SD-7 (Rautman and Engstrom 1996a, Figure 9)). 6.6.1.2 Porosity as a Surrogate A fundamental limitation of the rock properties modeling effort clearly is the use of porosity as a surrogate for derivative properties of more general interest to design and performance assessment analysts (Section 5). These derivative properties, such as matrix saturated hydraulic conductivity, are related to the principal modeled property only through a correlation coefficient, generally indicated as r. To the extent that the absolute value of r is less than one (recall that r = –1.0 implies a perfect inverse relationship), this modeling of a surrogate thus increases the uncertainty in the secondary properties. Nevertheless, use of a non-zero correlation coefficient combined with incorporation of spatial correlation in the modeling of those secondary properties works to decrease uncertainty in those secondary properties. This is in contrast to modeling methodologies that discount either or both of these observable statistical characteristics. Knowing that a particular region exhibits high porosity values most likely translates to higherthan- average matrix permeability in the same region (a positive r-value). However, actual measured values of derivative properties are not reproduced within the simulated (coregionalized) models in the same manner that measured porosity values located at a grid node are reproduced by construction. Again, the real issue is whether the modeled uncertainty in material properties translates to unacceptable uncertainty in an objective performance measure when evaluated over a number of statistically indistinguishable simulated models (Figure 3). Another limitation induced by the modeling process into the rock properties models through the porosity-as-a-surrogate mechanism involves those hydraulic conductivity values that are not correlated with porosity. Flint (1998, Figures 12(a) and (b)) describes a group of low-porosity samples that exhibit apparently random permeabilities with respect to their uniformly low porosity values, and Flint interprets these samples as exhibiting behavior consistent with the existence of microfractures, largely in vitrophyric rocks (lithostratigraphic units Tptrv1 and Tptpv3). The implication is that the permeability measured is not truly a matrix property, even though microfracture-related flow is measurable at the core scale. These erratic, out-of-porositycharacter permeability values have been modeled as a random overprint imposed only on extremely low-porosity (< 0.05) grid nodes (interpreted as representing vitrophyre). Because the entire sampled population of microfractured “vitrophyre-like” samples consists of a mere dozen or so individuals, it is simply impossible to determine if these values are spatially correlated in their own right. However, to the extent that the measured data truly represent vitrophyre as a rock type (geologically restricted to the upper (Tptrv1) and lower (Tptpv3) margins of the TSw unit), it is thus possible to generate microfractured permeability values at inappropriate spatial Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 115 of 150 locations within the TSw model unit. This limitation is presumed to be relatively minor, as the restriction of generating these values to extremely low-porosity grid nodes (typically between five and seven percent of the TSw grid) suggests that such rocks might be susceptible to microfracturing even though they would not belong to the vitrophyre-type small-scale lithostratigraphic units (see also Flint 1998, Figures 12(a) and (b)). Additionally, conditioning of the interior of the TSw model unit to porosity values substantially in excess of 0.05 produces models that are very unlikely to exhibit extremely low porosity values except near the margins where measured porosities of this magnitude are observed. A somewhat similar limitation affects the modeling of altered hydraulic conductivity values within the CHn and Tcp model units. Although the style of hydrous-mineral-phase alteration responsible for these reduced matrix permeability values is clearly correlated spatially (Section 6.4.7.2), it is unclear whether the permeability values themselves within those altered regions are correlated. Adequate, spatially distributed measured Ks data do not exist to provide an estimate of the spatial correlation structure of permeability itself in these materials. In any event, as there is no reliable relationship between porosity and matrix permeability for these altered specimens, it is unclear what would serve as a surrogate for modeling spatial continuity for these materials. Accordingly, the spatial distribution of altered permeability has been treated as random within the spatial envelope of altered rocks, and the affected grid nodes merely assigned values sampled from a normal population with the appropriate mean and variance. To the extent that altered permeabilities are, in fact, spatially correlated, this modeling approach increases uncertainty. Across the full suite of simulated models, however, the variability of rock properties (and presumably of process modeling results as well) is likely to be greater than had the properties been spatially correlated, thus allowing a quantitative evaluation of the consequences of that increased space of uncertainty. 6.6.1.3 Underestimation of Porosity and Hydraulic Conductivity in the TSw Model Unit Another limitation related somewhat to the use of porosity as a surrogate for hydraulic conductivity affects modeling of parts of the TSw model unit. Recall from Section 6.4.4.2 that for non-cored drill holes in this model unit, the available petrophysical data were able to provide only an estimate of lithophysal porosity, because the density logging tool is sensitive to the total amount of void space in the rock mass, including the influence of large lithophysal cavities. Accordingly, a “surrogate for porosity-as-a-surrogate” was adopted whereby the water-filled porosity data from the computed VWC trace were inserted into the matrix porosity data files for the named lithophysal zones only. Because matrix saturations, particularly in the crystal-rich lithophysal zone (Tptrl) and crystal-poor upper lithophysal zone (Tptpul) are less than one where these units are present above the static water level, substitution of the VWC data for matrix porosity underestimates the true matrix porosity (such as would be obtained from core) of the rock (see Figure 12). Coregionalization of matrix saturated hydraulic conductivity from porosity models conditioned to these lowered matrix porosity data produces models of Ks that are systematically somewhat low in localized regions. The impact of this limitation on the overall material property models is probably somewhat limited. First, the effect is limited almost exclusively to the upper lithophysal intervals (Tptrl and Tptpul), as matrix saturations within the crystal-poor lower lithophysal zone (Tptpll) are Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 116 of 150 typically sufficiently high that there is little mismatch between measured core porosity values and the water-filled porosity (VWC) measurements (see Figure 12). Second, the matter becomes an issue only for areas populated by non-cored drill holes. The immediate Yucca Mountain-ESF area is dominated by drill holes for which core samples were obtained (Figure 5). Close within the range of influence of these cored holes, the porosity (and hydraulic conductivity) models are strongly conditioned by the presence of laboratory measured matrix porosity data. Furthermore, examination of Figure 5 indicates that the principal holes lacking core, the WT- series, are located primarily in regions near the periphery or surrounding the volume modeled by the rock properties model. It is presumed that these regions are of less direct interest in the evaluation of the Yucca Mountain site. Third, the additional uncertainty caused by this substitution of VWC data for measured matrix porosity values has already been incorporated into the simulated models. For regions near the main part of the ESF where non-cored drill holes (such as H-5, Figure 5) compete with cored holes (e.g., SD-9), any discordance between the laboratory measurements and the VWC substitute will result in the simulation algorithm generating a wider range of simulated values (in the appropriate stratigraphic interval) than would be the case in the absence of that discordance. The space of uncertainty, as measured across the suite of simulations (Figure 3), has been increased, which is a realistic reflection of the state of knowledge associated with the limitations imposed by the use of non-core drilling techniques: the matrix porosity is unknown. 6.6.1.4 Use of Major Stratigraphic Units as Modeling Units A quite different, but potentially significant limitation of the approach used in development of the rock properties model is the use of composite major stratigraphic intervals and an internal stratigraphic coordinate system as the geometric basis for modeling. There are only four such model units; compare these to the virtual plethora of different lithostratigraphic units tabulated in Table 10. To the extent that the stratigraphic coordinate transformation for each of these major modeling units does not reposition equivalent parts of the model unit at the same stratigraphic position, the model will attempt to simulate continuity (improperly) between rocks formed at significantly different pressure-temperature conditions. These attempts will create an increase of uncertainty across the suite of replicate simulations as the simulation algorithm tries to resolve any inconsistency of the measured material properties in between drill hole locations. This limitation is probably of minimal effect within the major ash-flow units at Yucca Mountain, particularly for the Topopah Spring welded unit, which was effectively an instantaneous deposit of massive proportions. Although the rock unit thins southward away from its source, the same physical/chemical conditions responsible for the ultimate physical properties of the rock almost certainly varied with relative vertical position within the cooling rock mass. The same logical argument applies to a large extent to the Prow Pass Tuff modeling unit, and (to a lesser degree) to the multiple-cooling-unit Calico Hills nonwelded interval. The justification for treating the PTn modeling interval as a single rock properties unit is somewhat weaker, as this modeling unit contains two distinctly different, significant if only locally, developed pyroclastic-flow deposits (the Pah Canyon and Yucca Mountain Tuffs), separated by intervals of unrelated and even reworked volcanic materials. The decision to model a single PTn entity is based on two factors: (1) In general, the properties of the rocks within the Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 117 of 150 PTn unit are quite similar (almost all nonwelded tuffaceous materials, particularly within the region transected by the main part of the ESF), and especially in comparison with over- and underlying materials. (2) The individual genetic units are typically very thin, leading not only to a vastly increased bookkeeping task for selecting and tracking sample data but—more importantly—to markedly reduced statistical mass relevant to any one unit. Ultimately the choice to represent the PTn model unit as a whole was a pragmatic determination, presumably suitable for modeling at the site scale, although perhaps not appropriate for more detailed use. 6.6.1.5 Faulting, Erosion, and the Stratigraphic Coordinate System Another limitation related to the use of a stratigraphic coordinate system is that the presence of erosional unconformities or within-unit faulting will work to confound the petrologic and material-property equivalence of rocks assigned the same stratigraphic (vertical) coordinates. Faults are known to affect several of the drill holes at Yucca Mountain, specifically WT-1, WT-11, ONC-1, and UZ-7a. Drill hole p#1 is affected by erosion, and Moyer and Geslin (1995, pp. 8 to 31) report progressive lateral truncation of inferred depositional units within the Calico Hills Formation and Prow Pass Tuff across the model area. The effects of fault displacements have been included in the computation of stratigraphic coordinates based on the best available information, and in all cases, this compensation involves the same separations and uncertainties as the GFM (GFM3). The locations and extents of erosional complexities are probably less well constrained than the effects of faulting. However, in all cases, there should be no discontinuities between the adjustments to stratigraphic coordinates applied to the rock properties modeling effort and the representation of the GFM. A factor working to offset uncertainties related to the stratigraphic coordinate transformation is that all of the properties modeling activities were conducted within that conceptual and mathematical framework. Specifically, the quantitative description of spatial correlation behavior (variograms) was conducted after conversion to stratigraphic coordinates. If undetected faulting or erosion worked to juxtapose samples of differing rock properties, the observed range of spatial correlation should be reduced, and this higher lateral variability would be reflected in the simulated property models. For a given set of conditioning data, a lesser degree of spatial correlation will also translate into more variability across the suite of realizations, and thus more uncertainty is (accurately) reflected in the suite of simulated models. 6.6.2 Stochastic Uncertainty Assessment The entire rock properties modeling effort has been designed as one method for quantifying the geologic uncertainty—that which results from less-than-exhaustive site characterization—in the material properties used in downstream design and performance assessment analyses. The stochastic simulation process is intended to generate an arbitrary number of individual rock property models, each one of which reproduces the measured site characterization data and exhibits the full range of heterogeneity and spatial continuity observed in those data. This uncertainty in material properties then must be propagated into uncertainty in some particular performance measure, through some relevant consequence analysis, as depicted in the conceptual representation of a Monte Carlo modeling process presented in Figure 3. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 118 of 150 Because the emphasis in a geostatistical modeling analysis such as this one is on the joint characteristics of heterogeneity and spatial continuity, it logically follows that estimates of uncertainty are spatially variable as well. By theory and intuition, uncertainty as a statement of confidence is low in the immediate vicinity of observations, where the effect of conditioning the simulated models to measured property values leads to the construction of well-constrained probability density functions. At great distances from measured values, the probability density functions are less well constrained, which leads to the generation of more disparate values when considered across a suite of simulations. The space of uncertainty is greater. However, even at these locations, there is the constraint imposed by the statistical character of the data ensemble as a whole. At a minimum the unconditional probability density function is effectively equal to the histogram of all relevant measured values, in contrast to the uniform distribution of all physically possible values. The influence of spatial correlation in the modeling process further works to reduce uncertainty, regardless of the presence or absence of observed information in the vicinity. Because the simulated models are built sequentially, inclusion of previously simulated grid nodes in the local search neighborhood used in constructing the probability density function means that low values are likely to be generated in the vicinity of other low values and high values are likely to be generated near other highs. Note, however, that this relationship is not absolute. Because the value at each simulated grid node is derived by random sampling from the relevant probability density function, there is a finite chance within any one realization at each node of generating an “unusual” value from the tails of the distribution. Because the properties at any unsampled location are, in fact, not known, it is possible that the true property value might be quite different from the expectation. Across multiple simulations, it is therefore possible to account quantitatively for the performance or design consequences of that uncertainty by propagating variations in input properties through numerical process models to describe uncertainty in performance measures (Figure 3). 6.6.2.1 Expectations versus Individual Outcomes These comments regarding the representation of spatial uncertainty related to less-thanexhaustive site characterization apply to the individual stochastically simulated rock property models. In addition to these individual simulations or outcomes, summary type models have been created that correspond to the mathematical expectation of a probability distribution. These summary models were generated by computing the node-by-node arithmetic mean value of the suite of individual simulations (50 in the present case), and these are typically referred to as E-type models. Note, however, that simply because a value is mathematically the “most likely” does not mean that value is the appropriate value to use for a particular purpose. In fact, the “expected” value may be a very unlikely value for any actual realization or outcome. The mathematical expectation of one-half where zero is “heads” and one is “tails” does not mean that the result of tossing a coin is “most likely” to be standing on edge. In similar manner, an expected porosity value of 0.20 at a particular grid location does not necessarily imply that 0.20 is the value of porosity that should be entered into a numerical flow-and-transport model at this location. It means that over a (large) number of individual simulated outcomes, the average simulated Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 119 of 150 porosity at the location, individual values of which may have varied from, say, 5 to 45 percent, was 0.20. What those simulated values were, however, is discarded in the summary process. The limitations of E-type models, which are in common with the shortcomings of most interpolation-type modeling algorithms in representing the real world, are fairly well known. For example, the range of variability is always reduced with respect to variability of the input data (the so-called smoothing effect). Modeled values outside the limits of the measured values generally are impossible, even though it is rather unlikely that physical sampling has actually encountered the absolute highest or lowest value present in the real world. Also, the apparent strength of spatial correlation is typically much greater (a consequence of smoothing). And finally, the consequences of applying a numerical physical-process algorithm to a smoothed input property field may be very different from the results of applying that same algorithm to one (or more) of the underlying more realistically varying individual “outcome” models. This effect may be particularly important in performance assessment analyses, for which the “expected” behavior of the physical system is of less interest than high-consequence/low probability events representing the tails of a probability density function. Also, the flow consequences of smoothly varying (i.e., continuous) material properties may be quite different than the consequences of a more heterogeneous property representation. 6.6.2.2 Summary Uncertainty Models Uncertainty in rock material properties should be evaluated rigorously in terms of the consequences of that variability on a particular computed performance measure. However, the post-processing step that produces the E-type models is easily adapted to generate models of spatially distributed variability across the suite of simulations. Such a spatially varying representation may be thought of as a first-order “uncertainty model.” This type of summary model has been generated for this modeling activity as the node-by-node standard deviations of the 50 individual simulated rock property models. There is no particular reason to prefer the standard deviation approach other than general user familiarity. It would also be possible to represent in one summary model the variability of individual stochastic simulations through use of the inter-quartile range or the tenth-to-ninetieth-percentile difference. In almost any such uncertainty model, the magnitude of the spatially varying uncertainty will decrease effectively to zero at sample locations and increase to some maximum value at great distances from conditioning samples or in close proximity to samples exhibiting conflicting values. In all cases, it is important to separate conceptually the difference between spatial heterogeneity and uncertainty. PTn Model Unit–A block view of the PTn model unit in stratigraphic coordinates is presented in Figure 57 showing the uncertainty model of porosity for this unit. Uncertainty is generally low (shown in blue colors) in the immediate vicinity of drill holes containing conditioning measurements (Figure 4), and increases to higher values away from these locations. Uncertainty is spatially heterogeneous within the model as well. However, because most drill holes penetrate most of each unit, the general pattern of heterogeneity in uncertainty will be that exhibited on the top surface of the model unit. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 120 of 150 Figure 57. Uncertainty model showing E-type standard deviation of matrix porosity in the PTn model unit. Color scale is in porosity units as a fraction. Vertical exaggeration 10x. Legibility of the coordinates can be deduced from Figure 42 on page 95 of this document. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 121 of 150 TSw Model Unit–Block models presenting the uncertainty models of matrix and lithophysal porosity in the TSw model unit are presented in Figures 58 and 59, respectively. Uncertainty in the figures is calculated as the node-by-node standard deviation of the replicate simulated material property models. As anticipated, uncertainty is lowest in the vicinity of drill holes containing conditioning measurements (Figure 5), and it increases away from those drill hole locations. Note that there are slight differences in the uncertainty models for matrix and lithophysal porosity, even though drill hole coverage at this stratigraphic level is essentially identical. As is typical in many real-world data sets, the variance (standard deviation) is a function of the magnitude of the variable under consideration. CHn Model Unit–The uncertainty model for matrix porosity in the CHn model unit is presented in Figure 60, again as the node-by-node standard deviation of the 50 replicate simulated models of this property. Low values of uncertainty are indicated by the shades of blue, and these regions are associated with the drill holes containing conditioning data that penetrate this model unit (Figure 6). Spatial heterogeneity of uncertainty can be seen on the front face of the block model. However, the dominant pattern of uncertainty will be vertically downward associated with the vertical drill holes. Tcp Model Unit–A block view of uncertainty in matrix porosity for the Tcp model unit, as computed as the standard deviation of the replicate simulated models, is presented in Figure 61. Uncertainty is lowest in the immediate vicinity of the drill holes penetrating this unit (Figure 7), and increases away from these locations of conditioning data. It is interesting to note the marked increase in uncertainty within the Tcp model unit, in comparison with the uncertainty model for the CHn model unit presented previously in Figure 60. This increase is most noticeable in the northeast corner of the modeled volume (note also the increase in the maximum value of the standard deviation in these two figures, colored red on the color scale). The cause of this marked change in uncertainty is the loss of drill hole WT-16 from the Tcp data set (compare Figures 6 and 7, showing the drill hole locations). Additional increases in modeled uncertainty, particularly in the eastern portion of the area, in the Tcp model unit result from the loss of drill holes WT-14 and WT-15 (outside the modeled volume) and drill hole ONC-1 (within the volume). 6.7 MODEL VALIDATION A fundamental premise of the Monte Carlo simulation approach is that each individual realization is a plausible model of the unknown real world and that variation among the different stochastic realizations represents a probabilistic distribution of outcomes consistent with all that is known. Presumably, the only meaningful difference between realization 1 and realization N is that a different random number seed was used to initiate the simulation process (definition of random path and initialization of random number generators). Recalling that conditional simulations theoretically possess the attributes of data reproduction (including reproduction of ensemble statistical character) described in Section 6.4.6.2, it should be possible to test the validity of the individual simulated models in terms of statistical similarity to the data, by examining the three relevant criteria for simulated models specified in Section 6.4.6.2: 1. Reproduction of the known data values at the same locations within the model as represented by the real-world samples Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 122 of 150 Figure 58. Uncertainty model showing E-type standard deviation of matrix porosity in the TSw model unit. Color scale is in porosity units as a fraction. Vertical exaggeration 5x. Legibility of the coordinates can be deduced from Figure 42 on page 95 of this document. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 123 of 150 Figure 59. Uncertainty model showing E-type standard deviation of lithophysal porosity in the TSw model unit. Color scale is in porosity units as a fraction. Vertical exaggeration 10x. Legibility of the coordinates can be deduced from Figure 42 on page 95 of this document. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 124 of 150 Figure 60. Uncertainty model showing E-type standard deviation of matrix porosity in the CHn model unit. Color scale is in porosity units as a fraction. Vertical exaggeration 10x. Legibility of the coordinates can be deduced from Figure 42 on page 95 of this document. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 125 of 150 Legibility of the coordinates can be deduced from Figure 42 on page 95 of this document. Figure 61. Uncertainty model showing E-type standard deviation of matrix porosity in the Tcp model unit. Color scale is in porosity units as a fraction. Vertical exaggeration 10x. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 126 of 150 2. Reproduction of the full range of measurement variability, as represented by the univariate descriptive statistics of the known data values (the histogram) 3. Reproduction of the bivariate statistics, or two-point spatial correlation structure, of the input data (the variogram). Because of the large number of suites of simulated models generated for the rock properties model, it is impractical to present validation statistics for every material property-model unit combination. Instead, “illustrative” examples validating selected simulated suites will be presented here. Similar validation exercises for all the material property models constituting RPM3 may be found in the scientific notebook covering this modeling activity (Rautman 1999). 6.7.1 Validation of TSw Lithophysal Porosity and Whole-Rock Thermal Conductivity Criterion number 1 (above) for validating a descriptive model consists of reproducing the known data values at the same locations within the model as represented by the real-world samples. Satisfaction of this criterion can be demonstrated by extracting a vertical profile through the modeled volume at a drill hole location and comparing that property profile to the measured property values that were used to condition the simulation process. Although this comparison is, perhaps, somewhat simple-minded, a numerical property model that does not reproduce measured values at the locations of those values can be demonstrated to be a distortion of the known reality, and can potentially be considered invalidated. In practice, reproduction of measured drill hole data by a material property model is only approximate, principally because the model can only represents rock properties at discrete grid locations (656-ft/200-m horizontal spacing in RPM3). Depending upon the distance from the drill hole to the nearest set of grid nodes, a certain amount of completely expectable variation from the actual measured values may be introduced into the comparison. Recall also that a simulated model is generated by sampling at random from the locally conditioned probability density function appropriate for that grid node. Nevertheless, the model should reproduce fairly accurately the major features of the measured values, particularly when several individual simulations are considered. An “odd” property value sampled by chance from a tail of the conditional probability density function, even a well constrained one located near a conditioning datum, in one stochastic simulation is unlikely to be repeated in multiple simulations. Therefore, the summary E-type model, which summarizes the entire suite of simulated models, should reproduce the measured values fairly closely, even accounting for a modest mismatch between the locations of the drill hole and the vertical grid nodes comprising the profile. Figure 62 presents three drill hole profiles extracted from five different simulations of lithophysal porosity for the TSw model unit. For each profile, the detailed (nominal 3-ft spacing) measured porosity data are shown, as are the discretized model (~16-ft vertical spacing) porosities for each different simulation. Also shown for comparison is the E-type porosity model (also discretized, but shown as a continuous curve for clarity on the figure). The distance from each drill hole to the nearest set of grid nodes is shown on the figure for reference. Note that an argument could be made (for the five simulations shown) that the inter-simulation variability exhibited at drill hole G-4 (66 ft distant) appears somewhat less than the variation present at drill hole SD-7 (332 ft distant). Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 127 of 150 Figure 62. Lithophysal porosity profiles extracted from five simulated models of the TSw model unit, showing comparison of results with the measured lithophysal porosity values and with the results for the E-type model. (a) SD-7, (b) G-4, (c) H-5. Distance to Grid: 332 ft Porosity, as a fraction 0.0 0.2 0.4 0.6 Stratigraphic Depth 0 200 400 600 800 1000 (a) Distance to Grid: 259 ft Porosity, as a fraction 0.0 0.2 0.4 0.6 0 200 400 600 800 1000 Data Sim. 03 Sim. 17 Sim. 23 Sim. 34 Sim. 47 E-type (c) G-4: TSw Distance to Grid: 66 ft Porosity, as a fraction 0.0 0.2 0.4 0.6 0 200 400 600 800 1000 (b) Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 128 of 150 With respect to validation criterion 2, which deals with reproduction of the ensemble statistical character of the data, Figure 63 presents histograms of four different simulations of lithophysal porosity (parts (a) through (d)). Also shown for comparison are the histograms of the input conditioning data values (part (e) of the figure), and of the summary E-type model (part (f)). The figure clearly indicates quite exact reproduction of the input histogram by the several simulated models. This is particularly so given that the input histogram is, in fact, only a sampling of all the potential measurements of lithophysal porosity that could be obtained from the Yucca Mountain site. In contrast, the histogram of the E-type summary model is distinctly more normal in character, and the tails of the distribution (both high and low values) are noticeable compressed in comparison with the measured data. In fact, reduction of the high-valued tail is quite marked. This is despite the fact that the mean porosity of both simulations/data and the E-type model is virtually identical at 0.12. Figure 64 presents variograms summarizing the spatial continuity pattern of the simulated models of lithophysal porosity for the TSw model unit. Reproduction of the bivariate statistical character or spatial correlation structure of the measured data is criterion 3 proposed for validation of a descriptive model. These figures present the average variogram across all 50 replicate simulations, together with the plus/minus one standard deviation values. According to the figure, reproduction of the input model variogram (Table 15) is quite excellent in the vertical and north-south directions. Reproduction of the input model is somewhat less in the east-west direction of minimum horizontal connectivity, with the overall model indicating a longer range of correlation in this direction, and therefore less horizontal anisotropy than had been modeled. However, note that the average variogram, in fact, approximates the experimental variogram derived directly from the data at the longer lag separations, especially those around 12,000 to 13,000 ft. Presumably there is something about the conditioning data values themselves that (1) causes the experimental variogram to decrease from its higher values at slightly shorter lag distances and (2) constrains the simulation algorithm to produce a more continuous continuity structure in the east-west direction. Note that in each part of Figure 64, the computed variogram for the E-type summary model indicates a markedly longer range of correlation than either the experimental variogram or the variograms of the simulated models. This behavior is typical of E-type models, as the smoothing effect increases the apparent continuity of the values. 6.7.2 Validation of TSw Thermal Conductivity Models Validation of the simulated thermal conductivity models for the TSw model unit is complicated somewhat by the very nature of these coregionalized representations. Thermal conductivity is severely undersampled at Yucca Mountain (the TSw model unit is represented by only 35 samples (Table 19). Because of the way the coregionalization process has been implemented in this modeling exercise, there are no conditioning thermal conductivity data against which one can judge the validity of these models via criterion 1. Also, because of the mechanics of the coregionalization process itself, the spatial continuity patterns or variogram of the thermal conductivity values as a whole will be essentially identical to that of the lithophysal porosity models from which the coregionalization was performed. Criterion 3 is thus functionally not applicable as well. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 129 of 150 Figure 63. Histograms and summary statistics of four randomly selected individual simulations of lithophysal porosity in the TSw model unit (a)–(d), compared to the original porosity measurements (e) and the E-type summary (f). (a) (b) (c) (d) (e) (f) Frequency Porosity, as a fraction 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.000 0.040 0.080 0.120 TSw Lithophysal Porosity; Sim. 18 Number of Data89182 mean 0.1437 std. dev. 0.0675 coef. of var 0.4698 maximum0.5593 upper quartile0.1840 median 0.1390 lower quartile 0.0940 minimum0.0002 Frequency Porosity, as a fraction 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.000 0.040 0.080 0.120 TSw Lithophysal Porosity; Sim. 23 Number of Data89182 mean 0.1432 std. dev. 0.0695 coef. of var 0.4856 maximum 0.5491 upper quartile0.1850 median 0.1380 lower quartile 0.0920 minimum0.0006 Frequency Porosity, as a fraction 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.000 0.040 0.080 0.120 TSw Lithophysal Porosity; Sim. 34 Number of Data89182 mean 0.1425 std. dev. 0.0686 coef. of var 0.4814 maximum0.5502 upper quartile0.1840 median 0.1380 lower quartile 0.0920 minimum.00001 Frequency Porosity, as a fraction 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.000 0.040 0.080 0.120 TSw Lithophysal Porosity; Sim. 45 Number of Data89182 mean 0.1438 std. dev. 0.0685 coef. of var 0.4762 maximum0.5472 upper quartile0.1850 median 0.1390 lower quartile0.0940 minimum0.0010 Frequency Porosity, as a fraction 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.000 0.020 0.040 0.060 0.080 0.100 TSw Lithoph Poros (POREF) Number of Data11796 number trimmed235 mean 0.1462 std. dev. 0.0757 coef. of var 0.5180 maximum0.5510 upper quartile0.1900 median 0.1380 lower quartile0.0900 minimum0.0010 Frequency Porosity, as a fraction 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.000 0.040 0.080 0.120 0.160 TSw Lithoph. Porosity E-type Model Number of Data89182 mean 0.1423 std. dev. 0.0471 coef. of var 0.3311 maximum0.4214 upper quartile0.1725 median 0.1437 lower quartile0.1100 minimum0.0035 TRI-6115-965-0 TSwLithoporHists.wmf Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 130 of 150 Figure 65 addresses criterion 2, and the figure presents histograms and summary statistics for four individual coregionalized simulations, plus the target histogram reproduced from Figure 28. Again, similar to the case with lithophysal porosity, Figure 65 also presents the histogram and statistics of the E-type summary model. As was the case above for the underlying lithophysal porosity values, reproduction of the input measurements by the various simulated models is excellent. The normalization of the histogram of the E-type summary model (part (f) of the figure) is even more pronounced than was the case for the lithophysal porosity. One other criterion bears on validation of coregionalized models, such as these for thermal conductivity: the cross-variable correlation with the underlying primary property model. Figure 66 presents a crossplot of thermal conductivity versus lithophysal porosity for Figure 64. Reproduction of input variograms for simulated lithophysal porosity models for the TSw model unit: (a) stratigraphic vertical; (b) stratigraphic horizontal, azimuth = 0°; (c) stratigraphic horizontal, azimuth = 90°. Line with error bars is average variogram with plus/minus one standard deviation. Heavy solid line without error bars is input model. E-type model normalized by variance. Vertical Separation, in feet 0 100 200 300 400 Gamma 0.0 0.5 1.0 1.5 Mean Simulation Input Model Data E-type Model (a) TSwPorValVario1.wmf Separation, in feet 0 5000 10000 15000 20000 Gamma 0.0 0.5 1.0 1.5 Mean Simulation Input Model Data E-type Model (b) TSwPorValVario2.wmf Separation, in feet 0 5000 10000 15000 20000 Gamma 0.0 0.5 1.0 1.5 Mean Simulation Input Model Data E-type Model (c) TSwPorValVario3.wmf Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 131 of 150 Figure 65. Histograms and summary statistics of four randomly selected individual simulations of thermal conductivity in the TSw model unit (a)–(d), compared to the original thermal conductivity measurements (e) and the E-type summary model (f). Note: the identical statistics of parts (a) through (d) are correct. (a) (b) (c) (d) (e) (f) Frequency Thermal Conductivity, W/m-K 0.00 1.00 2.00 3.00 0.000 0.040 0.080 0.120 Simulated TSw Thermal K; Sim. 18 Number of Data89182 mean 1.1831 std. dev. 0.1822 coef. of var 0.1540 maximum2.1147 upper quartile1.3285 median 1.2034 lower quartile1.0550 minimum0.0000 Frequency Thermal Conductivity, W/m-K 0.00 1.00 2.00 3.00 0.000 0.040 0.080 0.120 Simulated TSw Thermal K; Sim. 23 Number of Data89182 mean 1.1831 std. dev. 0.1822 coef. of var 0.1540 maximum2.1147 upper quartile1.3285 median 1.2034 lower quartile 1.0550 minimum0.0000 Frequency Thermal Conductivity, W/m-K 0.00 1.00 2.00 3.00 0.000 0.040 0.080 0.120 Simulated TSw Thermal K; Sim. 34 Number of Data89182 mean 1.1831 std. dev. 0.1822 coef. of var 0.1540 maximum 2.1147 upper quartile1.3285 median 1.2034 lower quartile 1.0550 minimum 0.0000 Frequency Thermal Conductivity, W/m-K 0.00 1.00 2.00 3.00 0.000 0.040 0.080 0.120 Simulated TSw Thermal K; Sim. 45 Number of Data89182 mean 1.1831 std. dev. 0.1822 coef. of var 0.1540 maximum2.1147 upper quartile1.3285 median 1.2034 lower quartile1.0550 minimum0.0000 Frequency Thermal Conductivity, W/m-K 0.00 1.00 2.00 3.00 0.000 0.040 0.080 0.120 Simulated TSw Thermal K Number of Data6063 mean 1.1831 std. dev. 0.1819 coef. of var 0.1538 maximum1.5502 upper quartile1.3285 median 1.2032 lower quartile1.0550 minimum0 .6761 Frequency Thermal Conductivity, W/m-K 0.00 1.00 2.00 3.00 0.000 0.050 0.100 0.150 0.200 TSw Thermal Conductivity E-type Model Number of Data89182 mean 1.1831 std. dev. 0.0976 coef. of var 0.0825 maximum1.5098 upper quartile1.2500 median 1.1781 lower quartile1.1210 minimum0.7732 TRI-6115-967-0 TSwThKHists.wmf Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 132 of 150 one randomly selected coregionalized model. From Figure 26, the target correlation coefficient is –0.776 (based on an r2 value of 0.586 with a negative slope). As indicated on Figure 66(a), the actual correlation coefficient is –0.746, which appears quite good, considering that the target value is derived from a correlation of laboratory-sized core specimens whereas the simulated model is based on an implied large-scale correlation of lithophysal porosity with whole-rock thermal conductivity. Figure 66(b) presents the same cross-variable correlation plot for the E-type models of lithophysal porosity and thermal conductivity. For this comparison, the actual correlation coefficient is –0.983, markedly in excess of the target value of –0.746. The crosssimulation averaging process involved in constructing the E-type models has induced a nearly one-to-one correspondence between the two material properties. 6.7.3 Validation of Matrix Porosity in the CHn Model Unit Validation profiles for matrix porosity in the CHn model unit are presented in Figure 67, in order to allow evaluation of criterion 1 for this part of the rock properties model. Again, the format is to present vertical profiles extracted from 5 different simulated models at drill hole locations, and to compare those simulated values with both the measured RH porosity data and with the values extracted from the E-type summary model. Although there are fewer distinctive features in the porosity data profiles within the CHn model unit than there were for the TSw unit, reproduction of both the general magnitude of porosity and of what features there are is quite good. Note, particularly, reproduction of the relatively large Figure 66. Crossplot showing (a) cross-variable correlation between coregionalized thermal conductivity and the underlying lithophysal porosity simulation (5-percent subsample) and (b) cross-variable correlation for the equivalent E-type models. Compare Figure 26(b). Porosity, as a fraction 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Thermal Conductivity, in W/m-K 0 1 2 3 r ² = 0.557 r = -0.746 ThKlithoPorValid.wmf (a) Porosity, as a fraction 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Thermal Conductivity, in W/m-K 0 1 2 3 r2 = 0.966 r = -0.983 ThKlithoPorEtypeValid.wmf (b) Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 133 of 150 Figure 67. Matrix porosity profiles extracted from five simulated models of the CHn model unit, showing comparison of results with the measured lithophysal porosity values and with the results for the E-type model. (a) SD-7, (b) SD-9, (c) H-5. Distance to Grid: 332 ft Porosity, as a fraction 0.0 0.2 0.4 0.6 Stratigraphic Depth 0 100 200 300 400 (a) Distance to Grid: 259 ft Porosity, as a fraction 0.0 0.2 0.4 0.6 0 100 200 300 400 (b) Distance to Grid: 259 ft Porosity, as a fraction 0.0 0.2 0.4 0.6 0 100 200 300 400 Data Sim. 06 Sim. 13 Sim. 29 Sim. 31 Sim. 42 E-type (c) Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 134 of 150 wavelength porosity features in drill hole SD-7 (part (a) of Figure 67). Reproduction of shorter wavelength features deeper in that same drill hole are not particularly well reproduced. However, in this latter case, the wavelength of the porosity excursions is approximately the same as the grid discretization (5 m/16 ft). The missing data values (offscale) just below a stratigraphic depth of 300 feet in drill hole H-5 (part (c)) provide an interesting, if somewhat subtle, example of simulation in the vicinity of a data absence. It appears as though the variability among the 5 simulations shown right at this stratigraphic elevation is somewhat greater than it is several grid nodes higher on the profile. Although it is unclear that any firm conclusions in this regard should be drawn from examination of only 5 simulations, the proposed increase in variance at this level is consistent with simulation theory. Validation criterion 2 is addressed in Figure 68. Examination of the histograms associated with the four randomly selected simulated models indicates a very close resemblance to the histogram of the measured RH porosity data. The reduction of the tails of the statistical distribution for the E-type summary model is quite marked, as is the normalization of the histogram itself. Both of these features are consistent with the statistical distortions observed above with respect to lithophysal porosity in the TSw model unit. Variogram reproduction (validation criterion 3) for the simulated models of matrix porosity in the CHn model unit is presented in Figure 69. As was the case with the variogram comparison for the TSw model unit, the comparison is with the mean variogram computed over all 50 replicate stochastic realizations. Examination of the figure suggests that the input model has been reproduced quite well in all three principal directions of continuity. A small mismatch may be observed in part (c) of the figure for the minimum horizontal direction of correlation, although the degree of mismatch is less than was observed for lithophysal porosity. Again, consideration of the experimental (data) variogram with respect to the input model suggests that the simulated models may have been conditioned to produce a somewhat longer range of correlation than indicated by the mathematical modeling of the experimental data. Also, as was the case with the variograms for lithophysal porosity in the TSw model unit (Figure 64), the range of correlation exhibited by the E-type model is, in general, markedly longer than that expressed by the data or the simulated models. 6.7.4 Validation of Matrix Saturated Hydraulic Conductivity in the CHn Model Unit Validation of the derivative models of matrix saturated hydraulic conductivity in the CHn model unit is one step more complicated than the validation of the thermal conductivity models discussed in Section 6.7.2. The cause of this additional complication is that the nonwelded tuffs of the Calico Hills interval have been partially altered to hydrous-phase minerals, principally zeolites. Accordingly, the material properties of this interval partially reflect the original neardepositional processes, while in part reflecting the operation of the later-stage diagenetic processes of alteration. According to the conceptual model used for modeling this unit (and the underlying Prow Pass Tuff), porosity is largely unaffected by the alteration process (see Section 6.4.4.3 and 6.4.4.4). Presumably then, the statistical character of the matrix porosity models is essentially that of the measured data: a situation that was validated by the histograms and variograms shown in Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 135 of 150 Figure 68. Histograms and summary statistics of four randomly selected individual simulations of matrix porosity in the CHn model unit (a)–(d), compared to the original porosity measurements (e) and the E-type summary model (f). (a) (b) (c) (d) (e) (f) Frequency Porosity, as a fraction 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.000 0.040 0.080 0.120 CHn Porosity; Sim 2 Number of Data35088 mean 0.2612 std. dev. 0.0724 coef. of var 0.2772 maximum0.5241 upper quartile0.3070 median 0.2580 lower quartile0.2200 minimum0.0007 Frequency Porosity, as a fraction 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.000 0.040 0.080 0.120 CHn Porosity; Sim 13 Number of Data35088 mean 0.2599 std. dev. 0.0696 coef. of var 0.2676 maximum0.5400 upper quartile0.3030 median 0.2560 lower quartile0.2200 minimum0.0010 Frequency Porosity, as a fraction 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.000 0.040 0.080 0.120 CHn Porosity; Sim 29 Number of Data35088 mean 0.2624 std. dev. 0.0725 coef. of var 0.2762 maximum0.5120 upper quartile0.3070 median 0.2580 lower quartile0.2200 minimum0.0010 Frequency Porosity, as a fraction 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.000 0.040 0.080 0.120 CHn Porosity; Sim 41 Number of Data35088 mean 0.2609 std. dev. 0.0705 coef. of var 0.2703 maximum0.5365 upper quartile0.3060 median 0.2570 lower quartile 0.2200 minimum0.0002 Frequency Porosity, as a fraction 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.000 0.100 0.200 0.300 CHn Porosity E-type Model Number of Data35088 mean 0.2569 std. dev. 0.0368 coef. of var 0.1434 maximum0.4724 upper quartile0.2747 median 0.2529 lower quartile0.2378 minimum0.0448 Frequency Porosity, as a fraction 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.000 0.040 0.080 0.120 CHn POREF and RH data for FY99 Number of Data2711 number trimmed857 mean 0.2539 std. dev. 0.0753 coef. of var 0.2967 maximum0.5110 upper quartile0.3000 median 0.2520 lower quartile 0.2130 minimum0.0010 TRI-6115-956-0 CHnPorosHists.wmf Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 136 of 150 Figures 68 and 69. However, the statistical character of the derivative property, hydraulic conductivity should show a bimodal character (altered vs. unaltered) reflecting its two-stage genetic history. Figure 70 presents the now-familiar comparison of histograms representing four randomly selected coregionalized simulations of matrix saturated hydraulic conductivity for the CHn model unit. Also shown on the figure are the equivalent statistical distributions of the measured Ks values and the values generated for the E-type summary model. Examination of the figure suggests that the essential features of the measured data—a broad spread of high conductivity vitric conductivities corresponding to the full range of porosities in the CHn, a cluster of lower “zeolitic” conductivities, and a spike of arbitrarily valued non-flowing grid nodes representing the “no-flow” laboratory samples—have been well reproduced in the simulated models. The Figure 69. Reproduction of input variograms for simulated matrix porosity models for the CHn model unit: (a) stratigraphic vertical; (b) stratigraphic horizontal, azimuth = 0°; (c) stratigraphic horizontal, azimuth = 90°. Line with error bars is average variogram with plus/minus one standard deviation. Heavy solid line without error bars is input model. E-type model normalized by variance. Vertical Separation, in feet 0 100 200 300 400 Gamma 0.0 0.5 1.0 1.5 Mean Simulation Input Model Data E-type Model (a) CHnPorValVario1.wmf Separation, in feet 0 5000 10000 15000 20000 Gamma 0.0 0.5 1.0 1.5 Mean Simulation Input Model Data E-type Model (b) CHnPorValVario2.wmf Separation, in feet 0 5000 10000 15000 20000 Gamma 0.0 0.5 1.0 1.5 Mean Simulation Input Model Data E-type Model (c) CHnPorValVario3.wmf Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 137 of 150 Figure 70. Histograms and summary statistics of four randomly selected individual simulations of matrix saturated hydraulic conductivity in the CHn model unit (a)–(d), compared to the original Ks measurements (e) and the E-type summary model (f). (a) (b) (c) (d) (e) (f) Frequency log10 Saturated Hydraulic Conductivity, m/sec -15.0 -11.0 -7.0 -3.0 0.000 0.050 0.100 0.150 0.200 0.250 CHn Zeolitic Ksat; Sim. 2 Number of Data35088 mean -10.7578 std. dev. 2.5090 coef. of var undefined maximum -4.0000 upper quartile -9.9418 median -10.4019 lower quartile -14.0000 minimum -14.0000 Frequency log10 Saturated Hydraulic Conductivity, m/sec -15.0 -11.0 -7.0 -3.0 0.000 0.050 0.100 0.150 0.200 0.250 CHn Zeolitic Ksat; Sim. 13 Number of Data35088 mean -10.7425 std. dev. 2.4847 coef. of var undefined maximum -4.0001 upper quartile -9.9259 median -10.4009 lower quartile -14.0000 minimum -14.0000 Frequency log10 Saturated Hydraulic Conductivity, m/sec -15.0 -11.0 -7.0 -3.0 0.000 0.050 0.100 0.150 0.200 0.250 CHn Zeolitic Ksat; Sim. 29 Number of Data35088 mean -10.7647 std. dev. 2.4580 coef. of var undefined maximum -4.0011 upper quartile -9.9386 median -10.3998 lower quartile -14.0000 minimum -14.0000 Frequency log10 Saturated Hydraulic Conductivity, m/sec -15.0 -11.0 -7.0 -3.0 0.000 0.050 0.100 0.150 0.200 0.250 CHn Zeolitic Ksat; Sim. 41 Number of Data35088 mean -10.8020 std. dev. 2.3923 coef. of var undefined maximum -4.0003 upper quartile -9.9335 median -10.4034 lower quartile -14.0000 minimum -14.0000 Frequency log10 Saturated Hydraulic Conductivity, m/sec log10 Saturated Hydraulic Conductivity, m/sec -15.0 -11.0 -7.0 -3.0 0.000 0.100 0.200 0.300 CHn Ksat E-type Model Number of Data35088 mean -10.7911 std. dev. 1.7067 coef. of var undefined maximum -4.2669 upper quartile -11.1917 median -11.4231 lower quartile -11.6144 minimum -12.4858 Frequency -15.0 -11.0 -7.0 -3.0 0.000 0.050 0.100 0.150 0.200 0.250 CHn Measured Ksat Data Number of Data116 mean -11.0457 std. dev. 2.2273 coef. of var undefined maximum -4.7959 upper quartile -10.1600 median -10.5058 lower quartile -14.0000 minimum -14.0000 TRI-6115-955-0 CHnKsatHists.wmf Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 138 of 150 statistical character of the E-type summary model reflects the bimodal alteration state, although the two populations are more clearly separated and more normalized in shape. Also, the spike of “no-flow” grid nodes has been averaged out in the E-type post-processing step, as these nonflowing nodes represent a random process superimposed on each model of altered conductivity. Because they are neither spatially correlated nor conditioned to measured values, non-flowing nodes in one realization rarely are replicated at the same location in other simulations. The cross-variable correlation of hydraulic conductivity with matrix porosity in the CHn is presented for one representative simulated model in Figure 71(a). Part (b) of the figure represents the subset of grid nodes that are associated with unaltered materials, whereas part (c) of the figure represents those grid nodes that have been modeled as altered. What can be observed clearly in Figure 71(a) is the alteration-related segregation of conductivity values into the two populations (three, counting the non-flowing nodes arbitrarily assigned a Ks of 10–14 m/sec in part (c) of the figure). Dependence of unaltered hydraulic conductivity on porosity is clearly indicated in part (b); compare with Figure 20(b). The target r2 value was 0.604, which is somewhat higher than the modeled r2 value of 0.419. However, the target value is based on the correlation of measured Ks values with porosity for all units, not just for the CHn. The inclusion of other units in the target calculation involves a much wider range of porosity values than exhibited by data for the CHn model unit. The result of including correlated hydraulic conductivity values beyond the range of porosities observed for the CHn is to strengthen the apparent correlation. Note that there are only 12 unaltered Ks sample values for the CHn model unit, far too few data from which to develop a meaningful unit-specific target correlation. Figure 71(d) presents a similar cross-variable comparison of the E-type models of matrix porosity and saturated hydraulic conductivity for the CHn model unit. Although the resemblance to Figure 71(a) is unmistakable, it is also clear that the strength of the correlation between the two properties for the unaltered rock type (part (b) of the figure) is markedly greater and that the range of variability in porosity across the E-type model has been compressed. Variability of hydraulic conductivity has also been compressed for the altered lithology, and the “typical” or average altered conductivity has been decreased at least one log unit (compare the vertical position of the “horizontal” cloud of points with that shown in part (c)). Additionally, the cluster of values corresponding to the simulated non-flowing grid nodes in parts (a) and (c) of the figure has been eliminated by the averaging process. 6.7.5 Summary Statement Regarding Validation Based on the information presented in Figures 62 through 71, as well as on additional similar information presented in the scientific notebook associated with this modeling activity (Rautman 1999), it appears reasonable to conclude that the individual simulated rock properties models (including the coregionalized models of derivative properties) do meet the criteria that have been set out above with respect to “model validation,” and that the models as entire entities do closely resemble the input data used to construct these descriptive models. Specifically, the primary porosity models reproduce the input measurements used to condition the simulations within the limits imposed by the discretization of the models. Additionally, both the primary simulated models and the derivative coregionalized models exhibit the full range of material property variability exhibited by the conditioning data ensemble. Furthermore, the simulated Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 139 of 150 Porosity, as a fraction 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Matrix Saturated Hydraulic Conductivity, in log10 m/sec -14 -12 -10 -8 -6 -4 -2 (d) ZKsEtypeValPorXplot.wmf Figure 71. Crossplots showing cross-variable correlation between coregionalized saturated hydraulic conductivity and the underlying matrix porosity simulation for (a) entire modeled unit regardless of lithology and for (b) the unaltered and (c) altered facies. Simulation 29, 10- percent subsample. Compare Figures 20(b) and 21(b). (d) Equivalent crossplot for E-type models of matrix porosity and full-field saturated hydraulic conductivity. Porosity, as a fraction 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Matrix Saturated Hydraulic Conductivity, in log10 m/sec -14 -12 -10 -8 -6 -4 -2 r ² = 0.419 r = +0.647 ZKs029Unalt_PorXplot.wmf (b) Porosity, as a fraction 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Matrix Saturated Hydraulic Conductivity, in log10 m/sec -14 -12 -10 -8 -6 -4 -2 ZKs029Alt_PorXplot.wmf (c) Porosity, as a fraction 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Matrix Saturated Hydraulic Conductivity, in log10 m/sec -14 -12 -10 -8 -6 -4 -2 (a) ZKs029ValPorXplot.wmf Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 140 of 150 models reproduce the bivariate spatial correlation structure (variogram) observed in the data ensemble; the coregionalized models will reproduce this structure by construction. For the coregionalized models, the correlation coefficient between the derivative property and the underlying porosity simulation is reproduced within reasonable limits, given the fact that the target correlations are typically based on global correlations (without regard for model unit). The summary E-type models also reproduce the input conditioning measurements, but the statistical character of the models as a whole departs from the ensemble statistics of the underlying data. Univariate variability (histograms) is reduced, with the extreme tails of the distribution of values being reduced and the form of the distribution normalized to a greater or lesser extent. The spatial correlation structure (variograms) of the E-type models is similarly distorted. Continuity is observed to be somewhat greater in the summarized models. Additionally, cross-variable correlations appear to be strengthened by the averaging process implicit in the E-type models, approaching in some instances a one-to-one correlation. All of these three latter characteristics are completely expectable from the mechanics of the summary process. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 141 of 150 7. CONCLUSIONS The Rock Properties Model (RPM) is one component of the Integrated Site Model, which also includes the Geologic Framework Model and the Mineralogic Model. The RPM provides exhaustive, three-dimensional, discretized, numerical representations of several important hydrologic and thermal rock properties (porosity, bulk density, matrix saturated hydraulic conductivity, and thermal conductivity) that are intended for further use in numerical design and performance assessment analyses. The four composite modeling units defined for this analysis encompass the majority of the rocks within the unsaturated zone throughout the immediate vicinity of the potential repository at Yucca Mountain. 7.1 DESCRIPTIVE RESULTS Pre-modeling geostatistical analysis has produced a statistical and spatial continuity description of porosity for each of the four modeling units. Spatial correlation patterns are unit-specific and are relatively complex; nested variogram models were required to fit the observed continuity patterns. Spatial correlation is anisotropic, both vertically and in the stratigraphically horizontal plane. For porosity, a significant fraction of the total variability is reached within a distance of 2,000 to 6,000 feet, although another significant fraction exhibits a maximum range of between 10,000 and 38,000 feet. The spatial correlation of hydrous-phase mineral alteration was also evaluated for the lower two modeling units. Variogram patterns again are unit-specific, and the range of maximum correlation is observed to vary between 7,000 and 15,000 feet. Cross-variable correlations between porosity and the several secondary material properties have been analyzed on a global basis. The individual model components have been constructed using geostatistical simulation methods, conditioned to drill hole-based measurements of porosity, with the result that the individual stochastic models are heterogeneous, spatially correlated, and essentially indistinguishable statistically from the set of measured values. These simulated models of porosity have been provided as input to a linear coregionalization algorithm, which has been used to generate simulated models of the derivative material properties, such as bulk density, matrix hydraulic conductivity, and thermal conductivity. These derivative property models are also spatially correlated and heterogeneous and are close statistical replicas of the set of measured secondary property data. Cross-variable correlations exist, and the strength of these correlations is approximately that described by the input sample correlation coefficients. Unitspecific post-processing steps have been used to impart additional desired property attributes such as microfractured hydraulic conductivity associated with very low-porosity (vitrophyric) welded tuffs and different hydraulic conductivities associated with altered and unaltered nonwelded tuffs. Sets of 50 simulations for each material property in each modeling unit have been post-processed to provide spatially varying models of the node-by-node standard deviation of the replicate individual models. These representations of the spatially heterogeneous inter-simulation variability serve as a first-pass uncertainty model for each material property. The same post-processing run has also been used to generate a summary expected-value type (E-type) model, defined as the arithmetic mean value of the 50 input replicate simulated models. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 142 of 150 Although the averaging process distorts the statistical character of the conditioning data and simulated models, these E-type models are useful for visualizing large-scale heterogeneities in material properties across the site area. Visual examination of block diagrams and cross-sectional views through the various model components indicate that there is substantial, complex internal heterogeneity in all four modeling units, although the E-type models reveal material properties broadly compatible with the conceptual model of layered volcanic stratigraphy at Yucca Mountain. Features such as highporosity “lithophysal” intervals that correlate approximately with formally named lithophysal zones have been generated based solely on measured porosity values. However, it can be demonstrated that such high porosity intervals exist outside the formally named zones. Matrix porosity is also higher associated with these intervals of high, “lithophysal” porosity. Major heterogeneity of matrix saturated hydraulic conductivity, which reflects principally the distribution of hydrous-phase mineral alteration in the CHn and Prow Pass Tuff model units, is internally complex. Some interfingering of altered and unaltered material is indicated. 7.2 MODEL VALIDATION A “model validation” exercise has been conducted in which these exhaustive, descriptive models of material properties are compared with the sparse, spatially distributed input properties used to generate the models. Emphasis has been placed on examining the following three questions: • Do the models reproduce the input data at the locations of those data? • Do the models reproduce the full range of variability exhibited by the ensemble of measured values? • Do the models reproduce the spatial correlation structure exhibited by the data ensemble? A fourth question is relevant to the coregionalized models of derivative properties: • Do the models reproduce the observed correlation coefficient with porosity? The individual simulated models appear to pass these validation criteria, and are thus concluded to be good statistical models of the various material properties. The E-type summary models, however, although they still reproduce the input data at the locations of those data, are demonstrated to reproduce neither the univariate variability nor the observed spatial correlation structure of the data. Thus, although they may be useful for some purposes, the E-type summary models are deemed statistical distortions of the underlying data used in their construction. 7.3 UNCERTAINTY AND RESTRICTIONS The concept of “uncertainty” in terms of limitation on the “accuracy” of the RPM does not particularly apply to rock properties models in the usual sense, because the entire intent of the rock properties modeling activities is not to make location-specific predictions, but rather to provide an entire suite of equally plausible, realistic models that collectively describes the total space of uncertainty associated with the relevant material properties. No representation is made that the individual simulated models are “accurate” (except at data locations). However, these Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 143 of 150 same individual models are represented as geologically reasonable and consistent with all known information, and collectively they have been demonstrated (through validation) to be nearperfect replicas of the (geo)statistical character of the data ensemble used in their construction. It is fully intended that the consequences of the geologic uncertainty represented by the suites of models be evaluated through Monte Carlo-style flow-and-transport calculations or other suitable analyses. Methodological uncertainties involved in the RPM include those inherent in any statistical study, and these relate principally to uncertainties (error or bias) present in the input measurements used to condition the modeling algorithms. Measurement errors are presumed to be small in comparison with spatial heterogeneity. Two instances of measurement bias have been identified, but both have been compensated for to the extent possible. Some increased uncertainty is induced in the models of derivative material properties through the use of coregionalization (the concept of porosity-as-a-surrogate, based solely on the correlation coefficient) as the modeling algorithm (in contrast to conditional cosimulation). However, because of the overall modeling focus on defining spaces of uncertainty rather than on locally “accurate” predictions, the variability across a given suite of coregionalized models will be larger than it would otherwise, and hence that increased uncertainty is already accounted for. The use of stratigraphic coordinates for modeling purposes can act to increase uncertainty in the presence of significant erosion within a modeling unit or in the presence of faulted sections. However, instances of these confounding geologic features are known to be relatively few in the specific rocks being modeled, and the increased precision gained by correlating materials formed under similar temperature-pressure conditions afforded by the stratigraphic-coordinate conversion overall appear to be well worth the local distortions caused by erosion and/or faulting. Restrictions on subsequent use of the rock properties models are minimal and most are completely consistent with the Monte Carlo modeling methodology adopted as the conceptual framework for this activity. It is the individual simulated models (matched by run number across different properties to preserve the joint spatial character) that are intended for subsequent use in various relevant consequence analyses. The summary E-type models are not intended for general use in downstream analyses, except perhaps for individually justified special cases, because of their known distortion of the (geo)statistical character of the measured data. Evidently, the models are only useful for the stratigraphic intervals considered. Rocks belonging to the welded portion of the Tiva Canyon Tuff, as well as overlying units, were not modeled. Neither were deeper rocks, occurring principally in the saturated zone at Yucca Mountain, beginning at the stratigraphic level of the Bullfrog Tuff. The models are not valid beyond the defined lateral limits of the RPM, although the (geo)statistical descriptions and relevant measured values could be used to create other models not subject to these lateral limits. However, such extended modeling should be restricted to the outflow facies (in contrast to the caldera-margin or intra-caldera facies) of the several ash-flow tuff units involved in the four modeling units. A final restriction involves the saturated hydraulic conductivity models of the Topopah Spring welded modeling unit (only). These hydraulic conductivity models represent matrix permeability only, and they do not consider the effect of large lithophysal cavities that may be connected by Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 144 of 150 either matrix porosity or by fractures. Thus, the matrix hydraulic conductivity models are intended more as the upper permeability limit of the unsaturated hydraulic conductivity function at near-full matrix saturation, rather than as representing the permeability of the TSw model unit were it present within the saturated zone. As the TSw unit is almost invariably high within the unsaturated zone in the immediate vicinity of the potential repository (the modeled area), this restriction is of little practical effect. Note that the substitution of water-filled porosity (volumetric water content) within the upper lithophysal zone (in particular) for drill holes from which no core measurements are available leads to a slight underestimation of the true (fully saturated) matrix hydraulic conductivity in those regions. Again, however, the modeled values are reasonable estimates of the upper permeability limit of the unsaturated hydraulic conductivity function at the ambient matrix saturation. This document and its conclusions may be affected by technical product input information that requires confirmation. Any changes to the document or its conclusions that may occur as a result of completing the confirmation activities will be reflected in subsequent revisions. The status of the input information quality may be confirmed by review of the Document Input Reference System database. Several key inputs to the Rock Properties Model have To Be Verified (TBV) status; however, excluding them would prohibit construction of the model. Accordingly, outputs from the Rock Properties Model are also TBV. Those TBV data on which the RPM is based currently are being evaluated to verify their quality-assurance status. Data identified to be unqualified through this verification activity will be subject to an independent qualification process, in accordance with AP-2III.2Q, to ensure that those data on which the RPM is based are fully qualified. Because it is anticipated that all data on which the RPM is based ultimately will be qualified, there is no need at this time to develop criteria, pursuant to AP-3.10Q, by which to assess the impact and appropriateness of the use of unqualified data on the applicability or validity of the RPM. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 145 of 150 8. ATTACHMENTS ATTACHMENT I. Document Input Reference System (DIRS) - Removed. See electronic DIRS database. ATTACHMENT II. Sigma Plot Transform “STRATC4” ATTACHMENT III. GSLIB Routine “TRANS” ATTACHMENT IV. Validation of UNCERT Program “Vario” ATTACHMENT V. Validation of UNCERT Program “Variofit” Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 146 of 150 INTENTIONALLY LEFT BLANK Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 147 of 150 9. REFERENCES 9.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. Buesch, D.C.; Spengler, R.W.; Moyer, T.C.; and Geslin, J.K. 1996. Proposed Stratigraphic Nomenclature and Macroscopic Identification of Lithostratigraphic Units of the Paintbrush Group Exposed at Yucca Mountain, Nevada. Open-File Report 94-469. Denver, Colorado: U.S. Geological Survey. ACC: MOL.19970205.0061. CRWMS M&O (Civilian Radioactive Waste Management System Management and Operating Contractor) 1996. Determining Porosity, Balanced Water Content and Other Rock Properties from Geophysical Logs for the Modern Borehole Data Set at Yucca Mountain, Nevada, June 1996. BAAA00000-01717-0200-00013 REV 00. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19970210.0171. 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. Analysis and Modeling Report (AMR) Development Plan (DP) for Rock Properties Model Version 3.1. TDP-NBS-GS-000023. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19990917.0187. CRWMS M&O 1999d. Combined Porosity from Geophysical Logs. BAA000000-01717-0200- 00003 REV 00. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19991208.0437. Deutsch, C.V. and Journel, A.G. 1992. GSLIB Geostatistical Software Library and User's Guide. New York, New York: Oxford University Press. TIC: 224174. Deutsch, C.V. and Journel, A.G. 1998. GSLIB Geostatistical Software Library and User's Guide. Second 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 J. R. Dyer (DOE) 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. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 148 of 150 Flint, L.E. 1998. Characterization of Hydrogeologic Units Using Matrix Properties, Yucca Mountain, Nevada. Water-Resources Investigations Report 97-4243. Denver, Colorado: U.S. Geological Survey. ACC: MOL.19980429.0512. Journel, A.G. 1983. “Nonparametric Estimation of Spatial Distributions.” Mathematical Geology, 15, (3), 445-468. New York, New York: Plenum Publishing Corporation. TIC: 224736. Journel, A.G. and Huijbregts, C.J. 1978. Mining Geostatistics. New York, New York: Academic Press. TIC: 210128. Luster, G.R. 1986. Raw Materials for Portland Cement: Applications of Conditional Simulation of Coregionalization. Ph.D. dissertation. Ann Arbor, Michigan: University Microfilms International. TIC: 226080. Moyer, T.C. and Geslin, J.K. 1995. Lithostratigraphy of the Calico Hills Formation and Prow Pass Tuff (Crater Flat Group) at Yucca Mountain, Nevada. Open-File Report 94-460. Denver, Colorado: U.S. Geological Survey. ACC: MOL.19941208.0003. Nelson, P.H. 1996. Computation of Porosity and Water Content from Geophysical Logs, Yucca Mountain, Nevada. Open-File Report 96-078. Denver, Colorado: U.S. Geological Survey. ACC: MOL.19980529.0444. Rautman, C. 1999. Scientific Notebook, Geostatistical Modeling of Porosity and Derivative Properties (FY 1999). Scientific Notebook SNL-SCI-011. ACC: MOL.19990623.0025. Rautman, C.A and Engstrom, D.A. 1996a. Geology of the USW SD-7 Drill Hole Yucca Mountain, Nevada. SAND96-1474. Albuquerque, New Mexico: Sandia National Laboratories. ACC: MOL.19971218.0442. Rautman, C.A. and Engstrom, D.A. 1996b. Geology of the USW SD-12 Drill Hole Yucca Mountain, Nevada. SAND96-1368. Albuquerque, New Mexico: Sandia National Laboratories. ACC: MOL.19970613.0101. Rautman, C.A. and McKenna, S.A. 1997. Three-Dimensional Hydrological and Thermal Property Models of Yucca Mountain, Nevada. SAND97-1730. Albuquerque, New Mexico: Sandia National Laboratories. ACC: MOL.19980311.0317. Rautman, C. and McKenna, S. 1998. Geostatistical Modeling of Porosity and Derivative Properties for Fiscal Year 1998. Scientific Notebook SNL-SCI-006. ACC: MOL.19981027.0187. Sawyer, D.A.; Fleck, R. J.; Lanphere, M.A.; Warren, R.G.; Broxton, D.E.; and Hudson, M.R. 1994. “Episodic Caldera Volcanism in the Miocene Southwestern Nevada Volcanic Field: Revised Stratigraphic Framework, 40Ar/39Ar Geochronology, and Implications for Magmatism and Extension.” Geological Society of American Bulletin, 106, (10), 1304-1318. Boulder, Colorado: Geological Society of America. TIC: 222523. SNL (Sandia National Laboratories) 1997. 3-D Rock Properties Modeling for FY 1998. Work Agreement WA-0344, Rev. 00. Albuquerque, New Mexico: Sandia National Laboratories. ACC: MOL.19990823.0276. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 149 of 150 SNL 1999. 3-D Rock Properties Modeling for FY 1999. Work Agreement WA-0358, Rev. 00. Albuquerque, New Mexico: Sandia National Laboratories. ACC: MOL.19991001.0158. 9.2 PROCEDURES AP-3.10Q, Rev. 0, ICN 0. Analyses and Models. Washington, D.C.: U.S. Department of Energy, Office of Civilian Radioactive Waste Management. ACC: MOL.19990225.0335. AP-3.10Q, Rev. 1, ICN 0. Analyses and Models. Washington, D.C.: U.S. Department of Energy, Office of Civilian Radioactive Waste Management. ACC: MOL.19990702.0314. AP-SI.1Q, Rev. 0, ICN 0. Software Configuration Management. Washington, D.C.: U.S. Department of Energy, Office of Civilian Radioactive Waste Management. ACC: MOL.19990505.0125. AP-SI.1Q, Rev. 1, ICN 0. Software Configuration Management. Washington, D.C.: U.S. Department of Energy, Office of Civilian Radioactive Waste Management. ACC: MOL.19990630.0395. AP-SIII.2Q, Rev. 0, ICN 0. Qualification of Unqualified Data and the Documentation of Rationale for Accepted Data. Washington, D.C.: U.S. Department of Energy, Office of Civilian Radioactive Waste Management. ACC: MOL.19990702.0308. QAP-2-0, Rev. 4. Conduct of Activities. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19971031.0060. QAP-2-0, Rev. 5. Conduct of Activities. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19980826.0209. QAP-SIII-3, Rev. 2. Scientific Notebooks. Las Vegas, Nevada: CRWMS M&O. ACC: MOL.19971105.0040. SNL QAIP 19-1, Rev. 04. Software Quality Assurance Requirements. Albuquerque, New Mexico: Sandia National Laboratories. ACC: MOL.19981103.0547. SNL QAIP 20-2, Rev. 01. Scientific Notebooks. Albuquerque, New Mexico: Sandia National Laboratories. ACC: MOL.19990224.0438. SNL QAIP 20-2, Rev. 02. Scientific Notebooks. Albuquerque, New Mexico: Sandia National Laboratories. ACC: MOL.19990224.0508. 9.3 SOURCE DATA, LISTED BY DATA TRACKING NUMBER GS960708312132.002. Porosity, Water Content, Mineralogy and Other Data Derived from Geophysical Logs and Cores for 26 Boreholes. Submittal Date: 07/09/1996. GS980708312242.010. Physical Properties of Borehole Core Samples, and Water Potential Measurements Using the Filter Paper Technique, for Borehole Samples from USW WT-24. Submittal Date: 07/27/1998. Title: Rock Properties Model (RPM3.1) Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: 150 of 150 GS980808312242.014. Physical Properties of Borehole Core Samples and Water Potential Measurements Using the Filter Paper Technique for Borehole Samples from USW SD-6. Submittal Date: 08/11/1998. LA9910DB831321.001. Mineralogic Variation in Drill Holes. Submittal date: 10/29/1999. MO9510RIB00002.004. RIB ITEM: Stratigraphic Characteristics: Geologic/Lithologic Stratigraphy. Submittal Date: 06/26/1996. MO9804MWDGFM03.000. A 3-D Geological Framework and Integrated Site Sub-Model (GFM3.0) of Yucca Mountain, Nevada. Submittal date: 03/31/1998. MO9811MWDGFM03.000. Input Data To Geologic Framework Model GFM3.0. Submittal date: 11/30/1998. MO9901MWDGFM31.000. Geologic Framework Model. Submittal date: 01/06/1999. MO9910POROCALC.000. Combined Porosity from Geophysical Logs. Submittal date: 10/05/1999. MO9911INPUTRPM.000. Core Porosities, Bulk Densities, Particle Densities (RH & OD Values) Plus Saturated Hydraulic Conductivities and Associated Porosities. Submittal date: 10/12/1999. SNF40060198001.001. Unsaturated Zone Lithostratigraphic Contacts in Borehole USW WT-24. Submittal date: 10/15/1998. SNF40060298001.001. Unsaturated Zone Lithostratigraphic Contacts in Borehole USW SD-6. Submittal Date: 10/15/1998. SNL01A05059301.005. Laboratory Thermal Conductivity Data for Boreholes UE25 NRG-4, NRG-5; USW NRG-6 and NRG-7/7A. Submittal Date: 02/07/1996. SNL01A05059301.007. Calculated Porosities for Thermal Conductivity Specimens from Boreholes UE25 NRG-4, UE25 NRG-5, USW NRG-6, and USW NRG-7/7a. Submittal Date: 10/14/1998. 9.4 BASELINED SOFTWARE GSLIB V 1.4 SGSIM V 1.4. STN: 10110-1.4MSGSIMV1.40-00. GSLIB V 2.0 IK3D V 2.0. STN: 10122-2.0MIK3DV2.0-00. GSLIB V 1.4 NSCORE V 1.201. STN: 10109-1.4MNSCOREV1.201-01. EARTHVISION V 4.0. STN: 30035 V4.0. MATCHUP Version 12/5/98. STN: 10196-12598-00. MATCH Version 12/5/98. STN: 10195-12598-00. Title: Rock Properties Model (RPM3.1) Attachment I Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: I-1 of I-2 ATTACHMENT I. DOCUMENT INPUT REFERENCE SYSTEM (DIRS) REMOVED See electronic DIRS database Title: Rock Properties Model (RPM3.1) Attachment I Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: I-2 of I-2 Title: Rock Properties Model (RPM3.1) Attachment II Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: II-1 of II-10 ATTACHMENT II. SIGMA PLOT TRANSFORM “STRATC4” Title: Rock Properties Model (RPM3.1) Attachment II Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: II-2 of II-10 Title: Rock Properties Model (RPM3.1) Attachment II Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: II-3 of II-10 SIGMA PLOT TRANSFORM “STRATC4” 1. Software Routine Identification Software Name and Version Number: “stratc4”, ver. 4 Note: version number was not assigned originally, but was added for documentation to be included in the Rock Properties Model Analysis Model Report. However, the digit “4” in the routine name was intended to indicate that this is the fourth iteration of a software routine to perform a similar function. SRR Document Identification Number: none SRR Media Number: none 2. Description and Testing Documentation that the software routine or macro provides correct results for a specified range of input parameters • Description and Equations of Mathematical Models, Algorithms, and Numerical Solution Techniques, as Applicable: The equation that is implemented in software routine “stratc4” to convert drillhole measured depths to “stratigraphic coordinates,” more specifically to stratigraphic depths (StratDepth), is as follows (equation 2, p. 55, of this report): , where SampleDepth is the depth of the relevant sample measurement in a given drill hole, UnitTop and UnitBottom are the measured or projected top and bottom contacts of re relevant model unit at the drill hole location used to determine the thickness of the unit, and NominalThickness is the appropriate unit-specific scaling constant (Table 12 of this report). The conversion produces a “stratigraphic depth” such that the very top of the geologic unit in question is at a “depth” of zero, the very bottom of the unit is at a “depth” equal to the nominal stratigraphic thickness of the unit, and all depths in between are scaled proportionately to the nominal thickness. Stratigraphic coordinates are dimensionless. • Description of Software Routine Including the Execution Environment Software routine “stratc4” is a “transform” implemented within the industry-standard, commercial statistical and graphics software package, Sigma Plot (SPSS, Inc.; 233 South Wacker Drive, 11th Floor; Chicago, Illinois 60606-6307). The routine is essentially a “macro” (however, see next paragraph). Sigma Plot transform syntax does not appear to StratDepth SampleDepth UnitTop – UnitBottom UnitTop – --------------------------------------------------------------- NominalThickness · = d d Title: Rock Properties Model (RPM3.1) Attachment II Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: II-4 of II-10 be version specific, but this routine has been used for the current application within Sigma Plot, version 5.00, running under Windows NT, version 4. Use of the transform is intimately associated with the spreadsheet mechanics of Sigma Plot. Note that a “transform” in Sigma Plot parlance is different than a “macro.” A Sigma Plot macro consists of a recording of a sequence of mouse-based actions for creating or modifying a graph. In contrast, a transform performs one or more mathematical operations on a set of input data and provides essentially what amounts to an internal programming language for performing those operations. The Sigma Plot spreadsheets designed for use with software routine “stratc4” are intended to contain various sets of core-based and geophysics-based rock property measurements, and to present those data values graphically as drillhole profiles displaying the selected measurements as a function of depth. Accordingly, the true, measured depths associated with each set of property values are arranged in a columnar manner, with four blank columns to the right of the true-depth column (to be filled in with stratigraphic depths); property values are typically arranged to the right of the four blank columns as a logical grouping. Multiple sets of such depth and property data may be present in a single spreadsheet. Additionally, a column of top and bottom depth values corresponding to the top and bottom of each geologic unit of interest is provided elsewhere in the spreadsheet. The numeric identifiers of the true-depth and “tops” columns are required as manual input to the software routine (involves editing the transform listing within Sigma Plot). The transform is invoked according to Sigma Plot conventions, whereupon the appropriate unit tops and bottoms are read from the spreadsheet, and the coordinate-conversion equation is computed repeatedly for each depth value present in the specified true-depth column, with the results being inserted into the appropriate output columns (presumably the four blank colums set up for this purpose). • Description of Test Cases Any arbitrary set of “depth” values can be used to test the stratigraphic depth calculations performed by software routine “stratc4”. The only requirement is that some number of “true-depth” values be contained within the “top” and “bottom” depths of the “geologic unit” being considered. As an alternative, the following data set, which consists of the depths associated with selected laboratory core property measurements taken from drillhole USW SD-7 (the property data themselves are irrelevant and are not modified by the transform), may be used to validate the software routine. Note that within Sigma Plot, these data entries are associated with specific cells identified by row and column numbers. Title: Rock Properties Model (RPM3.1) Attachment II Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: II-5 of II-10 The column headed “unit designation” is not used by the transform, and it is included only for reader reference. “Tops” Depth Unit Designation “True-Depth” 316.0000 top PTn 302.8000 316.0000 306.1000 308.6000 384.3000 base PTn 312.1000 384.3000 top TSw 320.9000 323.7000 1308.0000 base TSw 330.1000 1308.0000 top CH 332.9000 338.9000 1621.5000 base CH 341.5000 1621.5000 top PP 368.0000 369.0000 2183.9000 base Prow 371.9000 2183.9000 top Bullfrog 386.3000 389.7000 392.7000 395.8000 399.1000 401.4000 404.6000 • Description of Test Results The following lines are copied from the Sigma Plot spreadsheet in which the stratigraphic coordinate conversion was computed. In the original, each entry is associated with a particular cell, identified by row and column numbers. True Depth PTn StratD TSw StratD CHn StratD Tcp StratD 302.8000 -38.6530 -88.2321 -1282.5518 -937.9090 306.1000 -28.9898 -84.6595 -1278.3413 -935.5619 308.6000 -21.6691 -81.9530 -1275.1515 -933.7838 312.1000 -11.4202 -78.1639 -1270.6858 -931.2945 320.9000 14.3485 -68.6370 -1259.4577 -925.0356 323.7000 22.5476 -65.6057 -1255.8852 -923.0441 330.1000 41.2884 -58.6771 -1247.7193 -918.4922 332.9000 49.4876 -55.6458 -1244.1467 -916.5007 338.9000 67.0571 -49.1502 -1236.4912 -912.2333 341.5000 74.6706 -46.3354 -1233.1738 -910.3841 368.0000 152.2694 -17.6464 -1199.3620 -891.5363 369.0000 155.1977 -16.5638 -1198.0861 -890.8250 371.9000 163.6896 -13.4243 -1194.3860 -888.7624 386.3000 205.8565 2.1652 -1176.0128 -878.5206 389.7000 215.8126 5.8461 -1171.6746 -876.1024 392.7000 224.5974 9.0939 -1167.8469 -873.9687 Title: Rock Properties Model (RPM3.1) Attachment II Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: II-6 of II-10 Several items are of note. First, the input data (first listing, above) represent only a very small portion of the relevant drillhole data set. Specifically, the “true depths” include only those associated with the PTn geologic unit, plus a few extra depths both above and below for illustrative purposes. As expected, true depths higher in the hole than the top of the PTn unit (given above as 316 ft), are converted to meaningless stratigraphic depths (“PTn StratD”) less than zero. Second, depths below the bottom of the PTn unit (given above as 384.3 ft) are converted to equally meaningless stratigraphic depths in excess of the nominal stratigraphic thickness of the PTn unit, which is 200.0 (dimensionless). Third, using similar logic, almost all stratigraphic depths for the TSw geologic unit (“TSw StratD”) are negative (those above the top of the TSw unit at 384.3 ft), as are all of the indicated stratigraphic depths for the CHn unit and Tcp unit (“CHn StratD” and “Tcp StratD”), as none of the measured depths of with this illustrative data set are associated with those deeper units. And finally, note as logically consistent, that the lowermost valid PTn stratigraphic depth (at a true depth of 371.9 ft) is immediately succeeded by the first valid TSw stratigraphic depth (at a true depth of 386.3 ft). Shaded values in the listing indicate numbers that are meaningless (though arithmetically correct) for one or more of the reasons discusssed above. A manual calculation, conducted by calculator (a Hewlett-Packard HP-11C; Hewlett- Packard Co.; 3000 Hanover Street; Palo Alto, CA 94304), for the first valid PTn “true depth” value (at 320.9 ft) is as follows: , which has been rounded up to 14.3485 by the display-format specifier in the results presented above (set within Sigma Plot). • Range of Input Parameter Values for Which the Results Were Verified There are no “input parameter values” specified by the user, other than the column identifiers specifying where to look for the unit tops and bottoms and for the “true depths” that are to be converted. In terms of input data values, any real number recognized by Sigma Plot is a mathematically valid entry for both the unit contacts (top and bottom) and for the “true depths.” Whether or not specific values are legitimate values for a particular drillhole or sample depth is another issue (discussion above). The user is responsible for providing geologically meaningful input values. 395.8000 233.6750 12.4499 -1163.8915 -871.7639 399.1000 243.3382 16.0225 -1159.6810 -869.4168 401.4000 250.0732 18.5125 -1156.7464 -867.7809 404.6000 259.4436 21.9768 -1152.6635 -865.5050 True Depth PTn StratD TSw StratD CHn StratD Tcp StratD StratDepth 320.9 316.0 – 384.3 316.0 – ------------------------------ 200.0 14.34846266 = · = Title: Rock Properties Model (RPM3.1) Attachment II Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: II-7 of II-10 • Identification of Any Limitations on Software Routine Applications or Validity There are no particular limitations to the software routine, per se. Any valid real number recognized by Sigma Plot is acceptable as input. However, the routine does not limit the calculation of the stratigraphic depth to the geologic unit of interest. Accordingly, the routine may generate statigrahic “depths” that are outside the legitimate range of zero through NominalThickness. Such values, although representing arithmetically correct results, are meaningless. The user is responsible for identifying and dealing properly with such values in later applications. 3. Supporting Information • Directory Listing of Executable and Data Files Because software routine “stratc4” is a macro designed to operate inside another software package, there is no “executable” per se. The Sigma Plot tranform capability operates in an “interpreter” mode, compiling and executing each line of the transform at runtime. There are also no “data files” per se. All input data is contained within the composite Sigma Plot spreadsheet and graphics file. • Computer Listing of Source Code, if Available The following listing of the transform code has been annotated (in a different font) to illustrate some salient features and required user-input values. jsv5D ;Macro to compute stratigraphic coordinates ; version 4 ; note: uses contacts of t/m units in spreadsheet ; bottom of one unit can be NONcoincident with top of next ctop = 1 ; column with tops <--User Input dp = 4 ; true depth column <--User Input nc1 = dp + 1 ; column in which to put StratZ for PTn nc2 = dp + 2 ; column in which to put StratZ for TSw nc3 = dp + 3 ; column in which to put StratZ for CHn nc4 = dp + 4 ; column in which to put StratZ for PP top1 = cell( ctop,2 ) ; PTn parameters bot1 = cell( ctop,4 ) nt1 = 200.0 <-- NominalThickness for PTn top2 = cell( ctop,5 ) ; TSw parameters bot2 = cell( ctop,7 ) nt2 = 1000.0 <-- NominalThickness for TSw top3 = cell( ctop,8 ) ; CHn parameters Title: Rock Properties Model (RPM3.1) Attachment II Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: II-8 of II-10 bot3 = cell( ctop,10 ) nt3 = 400.0 <-- NominalThickness for CHn top4 = cell( ctop,11) ;PP parameters bot4 = cell( ctop,13) nt4 = 400 <-- NominalThickness for Tcp thk1 = bot1 - top1 col(nc1 ) = ( col( dp ) - top1 ) / thk1 * nt1 thk2 = bot2 - top2 col(nc2 ) = ( col( dp ) - top2 ) / thk2 * nt2 thk3 = bot3 - top3 col(nc3 ) = ( col( dp ) - top3 ) / thk3 * nt3 thk4 = bot4 - top4 col(nc4 ) = (col( dp ) - top4 ) / thk4 * nt4 • Computer Listing of Test Data Input and Output, Identifying Software Routine Name and Version Number A printout (taken from an Acrobat portable-document-format [.pdf] file) of the relevant parts of the test Sigma Plot spreadsheet is reproduced below in figure . No “software names” or “version numbers” are associated with the internal appearance of the spreadsheet. Title: Rock Properties Model (RPM3.1) Attachment II Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: II-9 of II-10 Test Case for "stratc4" 3dm units FY97 3DM 3 Core Depth PTn StratC TSw StratC CHn StratC Tcp StratC 1 316.0000 top PTn 302.8000 -38.6530 -88.2321 -1282.5518 -937.9090 2 316.0000 306.1000 -28.9898 -84.6595 -1278.3413 -935.5619 3 308.6000 -21.6691 -81.9530 -1275.1515 -933.7838 4 384.3000 base PTn 312.1000 -11.4202 -78.1639 -1270.6858 -931.2945 5 384.3000 top TSw 320.9000 14.3485 -68.6370 -1259.4577 -925.0356 6 323.7000 22.5476 -65.6057 -1255.8852 -923.0441 7 1308.0000 base TSw 330.1000 41.2884 -58.6771 -1247.7193 -918.4922 8 1308.0000 top CH 332.9000 49.4876 -55.6458 -1244.1467 -916.5007 9 338.9000 67.0571 -49.1502 -1236.4912 -912.2333 10 1621.5000 base CH 341.5000 74.6706 -46.3354 -1233.1738 -910.3841 11 1621.5000 top PP 368.0000 152.2694 -17.6464 -1199.3620 -891.5363 12 369.0000 155.1977 -16.5638 -1198.0861 -890.8250 13 2183.9000 base Prow 371.9000 163.6896 -13.4243 -1194.3860 -888.7624 14 2183.9000 top Bullfrog 386.3000 205.8565 2.1652 -1176.0128 -878.5206 15 389.7000 215.8126 5.8461 -1171.6746 -876.1024 16 392.7000 224.5974 9.0939 -1167.8469 -873.9687 17 395.8000 233.6750 12.4499 -1163.8915 -871.7639 18 399.1000 243.3382 16.0225 -1159.6810 -869.4168 19 401.4000 250.0732 18.5125 -1156.7464 -867.7809 20 404.6000 259.4436 21.9768 -1152.6635 -865.5050 Figure II-1. Printed spreadsheet cells from Sigma Plot file containing the test case data and outp (reproduced from a .pdf file). Columns are numbered sequentially from the left (see col. “3” but alphanumeric titles may be entered for user convenience. Title: Rock Properties Model (RPM3.1) Attachment II Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: II-10 of II-10 INTENTIONALLY LEFT BLANK Title: Rock Properties Model (RPM3.1) Attachment III Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: III-1 of 18 ATTACHMENT III. GSLIB ROUTINE “TRANS” Title: Rock Properties Model (RPM3.1) Attachment III Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: III-2 of III-18 Title: Rock Properties Model (RPM3.1) Attachment III Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: III-3 of 18 Software routine TRANS, version 1.3, is used directly from the geostatistical software library, GSLIB (Deutsch and Journel 1992, pp. 213–214). The routine is intended to convert any particular univariate distribution of values to match that of a different “reference” distribution. For the rock properties modeling exercise, TRANS was used to transform coregionalized, standard-normal (m=0, s2=1) models to match the target statistics of the desired secondary property in order to create the full-field derivative-property models. Documentation of TRANS was inadvertently omitted from the scientific notebooks covering this work, and hence this material is included as this attachment. Source Code Listing program trans C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%% C % C Copyright (C) 1992 Stanford Center for Reservoir Forecasting. All % C rights reserved. Distributed with: C.V. Deutsch and A.G. Journel. % C ‘‘GSLIB: Geostatistical Software Library and User’s Guide,’’ Oxford % C University Press, New York, 1992. % C % C The programs in GSLIB are distributed in the hope that they will be % C useful, but WITHOUT ANY WARRANTY. No author or distributor accepts % C responsibility to anyone for the consequences of using them or for % C whether they serve any particular purpose or work at all, unless he % C says so in writing. Everyone is granted permission to copy, modify % C and redistribute the programs in GSLIB, but only under the condition % C that this notice and the above copyright notice remain intact. % C % C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%% c----------------------------------------------------------------------- c c Univariate Transformation c ************************* c c Reads in a reference distribution and a number of other distributions c and then transforms the values in each of the second distributions c such that their histograms match that of the reference distribution. c c c INPUT/OUTPUT Parameters: c c distin the input file with the reference distribution c ivr,iwt parameters to read in reference distribution c tmin,tmax trimming limits c datafl the input file with the data to be transformed c ivrd,iwtd parameters to read in reference distribution c tmin,tmax trimming limits c nxyz,nsim number of points to transform at a time and the c number of times to expect nxyz in the file c outfl the output file for transformed data c zmin,zmax minimum and maximum data values Title: Rock Properties Model (RPM3.1) Attachment III Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: III-4 of III-18 c ltail,ltpar option to handle values in lower tail c utail,utpar option to handle values in upper tail c c c c Original: C.V. Deutsch Date: May 1991 c----------------------------------------------------------------------- parameter(MAXREF=25000, MAXDAT=100000, MV=20, EPSLON=1.0e-12, + VERSION=1.300) implicit real*8 (a-h,o-z) character distin*40,datafl*40,outfl*40,str*40 real*8 rcdf(MAXREF),rvr(MAXREF),dcdf(MAXDAT),dvr(MAXDAT), + indx(MAXDAT),var(MV),ltpar,utpar integer ltail,utail logical testfl data lin/1/,lout/2/ c c Note VERSION number before anything else: c write(*,9999) VERSION 9999 format(/’ TRANS Version: ‘,f5.3/) c c Get the name of the parameter file - try the default name if no input: c write(*,*) ‘Which parameter file do you want to use?’ read (*,’(a40)’) str if(str(1:1).eq.’ ‘)str=’trans.par ‘ inquire(file=str,exist=testfl) if(.not.testfl) then write(*,*) ‘ERROR - the parameter file does not exist,’ write(*,*) ‘ check for the file and try again ‘ stop endif open(lin,file=str,status=’OLD’) c c Find Start of Parameters: c 1 read(lin,’(a4)’,end=97) str(1:4) if(str(1:4).ne.’STAR’) go to 1 c c Read Input Parameters: c read(lin,’(a40)’,err=97) distin write(*,*) ‘ reference distribution: ‘,distin read(lin,*,err=97) ivr,iwt write(*,*) ‘ columns: ‘,ivr,iwt read(lin,*,err=97) tmin,tmax write(*,*) ‘ trimming limits: ‘,tmin,tmax read(lin,’(a40)’,err=97) datafl write(*,*) ‘ data file: ‘,datafl read(lin,*,err=97) ivrd,iwtd write(*,*) ‘ columns: ‘,ivrd,iwtd read(lin,*,err=97) tmind,tmaxd Title: Rock Properties Model (RPM3.1) Attachment III Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: III-5 of 18 write(*,*) ‘ trimming limits: ‘,tmind,tmaxd read(lin,*,err=97) nxyz,nsim write(*,*) ‘ number to transform: ‘,nxyz,nsim read(lin,’(a40)’,err=97) outfl write(*,*) ‘ output file: ‘,outfl read(lin,*,err=97) zmin,zmax write(*,*) ‘ data limits: ‘,zmin,zmax read(lin,*,err=97) ltail,ltpar write(*,*) ‘ lower tail option: ‘,ltail,ltpar read(lin,*,err=97) utail,utpar write(*,*) ‘ upper tail option: ‘,utail,utpar close(lin) c c Check for error situation: c if(ltail.ne.1.and.ltail.ne.2) then write(*,*) ‘ERROR invalid lower tail option ‘,ltail write(*,*) ‘ only allow 1 or 2 - see manual ‘ stop endif if(utail.ne.1.and.utail.ne.2.and.utail.ne.4) then write(*,*) ‘ERROR invalid upper tail option ‘,ltail write(*,*) ‘ only allow 1,2 or 4 - see manual ‘ stop endif if(utail.eq.4.and.utpar.lt.1.0) then write(*,*) ‘ERROR invalid power for hyperbolic tail’,utpar write(*,*) ‘ must be greater than 1.0!’ stop endif c c Read in the reference distribution: c inquire(file=distin,exist=testfl) if(.not.testfl) then write(*,*) ‘ERROR: No reference distribution file’ stop endif open(lin,file=distin) c c Proceed with reading in distribution: c read(lin,’(a)’,err=98) str read(lin,*,err=98) nvari do 2 i=1,nvari 2 read(lin,’()’,err=98) c c Read as much data as possible: c ncut = 0 nt = 0 tcdf = 0 3 read(lin,*,end=4,err=98) (var(j),j=1,nvari) c Title: Rock Properties Model (RPM3.1) Attachment III Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: III-6 of III-18 c Trim this data? c if(var(ivr).lt.tmin.or.var(ivr).ge.tmax) then nt = nt + 1 go to 3 endif if(iwt.ge.1) then if(var(iwt).le.EPSLON) then nt = nt + 1 go to 3 endif endif c c Accept this data: c ncut = ncut + 1 if(ncut.gt.MAXREF) then write(*,*) ‘ERROR: exceeded available storage for reference’ write(*,*) ‘ have ‘,MAXREF,’ available’ stop endif rvr(ncut) = var(ivr) if(iwt.ge.1) then rcdf(ncut) = var(iwt) else rcdf(ncut) = 1.0 endif tcdf = tcdf + rcdf(ncut) c c Go back for another data: c go to 3 4 close(lin) write(*,*) write(*,*) ‘ TRANS: there were ‘,ncut,’ reference values’ write(*,*) ‘ trimmed ‘,nt,’ values’ write(*,*) ‘ the total weight ‘,tcdf c c Sort the Reference Distribution and Check for error situation: c call sortem(1,ncut,rvr,1,rcdf,c,d,e,f,g,h) if(ncut.le.1.or.tcdf.le.EPSLON) then write(*,*) ‘ERROR: too few data or too low weight’ stop endif if(utail.eq.4.and.rvr(ncut).le.0.0) then write(*,*) ‘ERROR can not use hyperbolic tail with ‘ write(*,*) ‘ negative values! - see manual ‘ stop endif if(zmin.gt.rvr(1)) then write(*,*) ‘ERROR zmin should be no larger than smallest’ write(*,*) ‘ refernce value’ write(*,*) ‘ zmin = ‘,zmin,’ vr1 ‘,rvr(1) Title: Rock Properties Model (RPM3.1) Attachment III Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: III-7 of 18 stop endif if(zmax.lt.rvr(ncut)) then write(*,*) ‘ERROR zmax should be no less than the largest’ write(*,*) ‘ reference value’ write(*,*) ‘ zmax = ‘,zmax,’ vrnt ‘,rvr(ncut) stop endif c c Turn the (possibly weighted) distribution into a cdf that is useful: c oldcp = 0.0 cp = 0.0 tcdf = 1.0 / tcdf do 5 i=1,ncut cp = cp + rcdf(i) * tcdf rcdf(i) =(cp + oldcp) * 0.5 oldcp = cp 5 continue c c Get median and write some info to the screen: c half = 0.5 call locate(rcdf,ncut,half,j) gmedian = xlinint(rcdf(j),rcdf(j+1),rvr(j),rvr(j+1),half) write(*,*) ‘ the median ‘,gmedian write(*,*) c c Get the output and distribution files ready: c inquire(file=datafl,exist=testfl) if(.not.testfl) then write(*,*) ‘ERROR ‘,datafl,’ does not exist!’ write(*,*) ‘ you need some distributions to work with’ stop endif open(lin,file=datafl) open(lout,file=outfl) read(lin,’(a40)’,err=96) str read(lin,*,err=96) nvari do 6 i=1,nvari 6 read(lin,’(a40)’,err=96) str write(lout,110) 110 format(‘Trans output’,/,’1’,/,’transformed variable’) c c MAIN Loop over all the increments to transform: c do 10 isim=1,nsim c c Read in the data values: c tcdf = 0.0 num = 0 do 11 i=1,nxyz Title: Rock Properties Model (RPM3.1) Attachment III Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: III-8 of III-18 read(lin,*,end=14,err=98) (var(j),j=1,nvari) num = num + 1 if(num.gt.MAXDAT) then write(*,*) ‘ERROR: exceeded storage for data’ write(*,*) ‘ have ‘,MAXDAT,’ available’ stop endif dvr(num) = var(ivrd) indx(num) = dble(real(num)) if(dvr(num).ge.tmind.and.dvr(num).lt.tmaxd) then dcdf(num) = 1.0 if(iwtd.ge.1) dcdf(num) = var(iwtd) tcdf = tcdf + dcdf(num) endif 11 continue 14 if(tcdf.le.EPSLON) then write(*,*) ‘ERROR: no data’ stop endif c c Turn the (possibly weighted) data distribution into a useful cdf: c call sortem(1,num,dvr,2,dcdf,indx,d,e,f,g,h) oldcp = 0.0 cp = 0.0 tcdf = 1.0 / tcdf do 12 i=1,num if(dvr(i).ge.tmind.and.dvr(i).lt.tmaxd) then cp = cp + dcdf(i) * tcdf dcdf(i) =(cp + oldcp) * 0.5 oldcp = cp endif 12 continue c c Now, get the right order back: c call sortem(1,num,indx,2,dcdf,dvr,d,e,f,g,h) c c Go through all the data back transforming them to the reference CDF: c do 13 i=1,num if(dvr(i).ge.tmind.and.dvr(i).lt.tmaxd) then zval = getz(dcdf(i),ncut,rcdf,rvr,zmin,zmax, + ltail,ltpar,utail,utpar) else zval = -999.0 endif write(lout,101) zval 101 format(f9.4) 13 continue c c END Main loop: c 10 continue Title: Rock Properties Model (RPM3.1) Attachment III Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: III-9 of 18 stop 96 stop ‘ERROR in distribution file!’ 97 stop ‘ERROR in parameter file!’ 98 stop ‘ERROR in global data file!’ end double precision function getz(cdfval,ncut,rcdf,rvr,zmin,zmax, + ltail,ltpar,utail,utpar) c----------------------------------------------------------------------- c c Get a Z value for a give CDF c **************************** c c Performs specified interpolation in the lower and upper tail. c c c c c c c Original: C.V. Deutsch May 1991 c----------------------------------------------------------------------- parameter(EPSLON=1.0e-12,UNEST=-1.0) implicit real*8 (a-h,o-z) real*8 rcdf(ncut),rvr(ncut),ltpar,utpar,lambda integer ltail,utail,cclow,cchigh zero = 0.0 one = 1.0 c c Figure out part of distribution 0=lower tail, 1=middle, 2=upper tail: c ipart = 1 if(cdfval.le.rcdf(1)) ipart = 0 if(cdfval.ge.rcdf(ncut)) ipart = 2 c c ARE WE IN THE LOWER TAIL? c if(ipart.eq.0) then getz = rvr(1) if(ltail.eq.1) then c c Linear interpolation to lower limit “zmin”: c getz = powint(zero,rcdf(1),zmin,rvr(1),cdfval,one) else if(ltail.eq.2) then c c Power model interpolation to lower limit “zmin”: c cpow = 1.0/ltpar getz = powint(zero,rcdf(1),zmin,rvr(1),cdfval,cpow) endif c Title: Rock Properties Model (RPM3.1) Attachment III Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: III-10 of III-18 c OR, ARE WE IN THE MIDDLE? c else if(ipart.eq.1) then call locate(rcdf,ncut,cdfval,cclow) cchigh = cclow + 1 getz = xlinint(rcdf(cclow),rcdf(cchigh), + rvr(cclow),rvr(cchigh),cdfval) c c OR, ARE WE IN THE UPPER TAIL? c else if(ipart.eq.2) then c c Linear interpolation to upper limit zmax? c if(utail.eq.1) then getz = xlinint(rcdf(ncut),one,rvr(ncut), + zmax,cdfval) else if(utail.eq.2) then cpow = 1.0 / utpar getz = powint(rcdf(ncut),one,rvr(ncut), + zmax,cdfval,cpow) else if(utail.eq.4) then c c Fit a Hyperbolic Distribution? Figure out “lambda”: c lambda = (rvr(ncut)**utpar)*(1.0-rcdf(ncut)) getz = (lambda/(1.0-cdfval))**(1.0/utpar) endif endif c c Check for bounds and return: c if(getz.lt.zmin) getz = zmin if(getz.gt.zmax) getz = zmax return end double precision function xlinint(xlow,xhigh,ylow,yhigh,xval) c----------------------------------------------------------------------- c c Linearly interpolates the value of y between (xlow,ylow) and c (xhigh,yhigh) for a value of x. c c----------------------------------------------------------------------- parameter(EPSLON=1.0e-12) real*8 xlow,xhigh,ylow,yhigh,xval if((xhigh-xlow).le.EPSLON) then xlinint = xlow else xlinint = ylow + (yhigh-ylow)*(xval -xlow) / (xhigh-xlow) end if return Title: Rock Properties Model (RPM3.1) Attachment III Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: III-11 of 18 end double precision function powint(xlow,xhigh,ylow,yhigh,xval,pow) c----------------------------------------------------------------------- c c Power interpolate the value of y between (xlow,ylow) and (xhigh,yhigh) c for a value of x and a power pow. c c----------------------------------------------------------------------- parameter(EPSLON=1.0e-20) implicit real*8 (a-h,o-z) powint = ylow + (yhigh-ylow)* + (((xval-xlow)/amax1(EPSLON,(xhigh-xlow)))**pow) return end subroutine locate(xx,n,x,j) c----------------------------------------------------------------------- c c Given an array “xx” of length “n”, and given a value “x”, this routine c returns a value “j” such that “x” is between xx(j) and xx(j+1). xx c must be monotonic, either increasing or decreasing. j=0 or j=n is c returned to indicate that x is out of range. c c Bisection Concept From “Numerical Recipes”, Press et. al. 1986 pp 90. c----------------------------------------------------------------------- real*8 xx(n),x c c Initialize lower and upper methods: c jl = 0 ju = n c c If we are not done then compute a midpoint: c 10 if(ju-jl.gt.1) then jm = (ju+jl)/2 c c Replace the lower or upper limit with the midpoint: c if((xx(n).gt.xx(1)).eqv.(x.gt.xx(jm))) then jl = jm else ju = jm endif go to 10 endif c c Return with the array index: c Title: Rock Properties Model (RPM3.1) Attachment III Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: III-12 of III-18 j = jl return end subroutine sortem(ib,ie,a,iperm,b,c,d,e,f,g,h) c----------------------------------------------------------------------- c c Quickersort Subroutine c ********************** c c This is a subroutine for sorting a real array in ascending order. This c is a Fortran translation of algorithm 271, quickersort, by R.S. Scowen c in collected algorithms of the ACM. c c The method used is that of continually splitting the array into parts c such that all elements of one part are less than all elements of the c other, with a third part in the middle consisting of one element. An c element with value t is chosen arbitrarily (here we choose the middle c element). i and j give the lower and upper limits of the segment being c split. After the split a value q will have been found such that c a(q)=t and a(l)<=t<=a(m) for all i<=l7 no other array is permuted. c c b,c,d,e,f,g,h arrays to be permuted according to array a. c c OUTPUT PARAMETERS: c c a = the array, a portion of which has been sorted. c c b,c,d,e,f,g,h =arrays permuted according to array a (see iperm) c c NO EXTERNAL ROUTINES REQUIRED: Title: Rock Properties Model (RPM3.1) Attachment III Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: III-13 of 18 c c----------------------------------------------------------------------- implicit real*8 (a-h,o-z) dimension a(*),b(*),c(*),d(*),e(*),f(*),g(*),h(*) c c The dimensions for lt and ut have to be at least log (base 2) n c integer lt(64),ut(64),i,j,k,m,p,q c c Initialize: c j = ie m = 1 i = ib iring = iperm+1 if (iperm.gt.7) iring=1 c c If this segment has more than two elements we split it c 10 if (j-i-1) 100,90,15 c c p is the position of an arbitrary element in the segment we choose the c middle element. Under certain circumstances it may be advantageous c to choose p at random. c 15 p = (j+i)/2 ta = a(p) a(p) = a(i) go to (21,19,18,17,16,161,162,163),iring 163 th = h(p) h(p) = h(i) 162 tg = g(p) g(p) = g(i) 161 tf = f(p) f(p) = f(i) 16 te = e(p) e(p) = e(i) 17 td = d(p) d(p) = d(i) 18 tc = c(p) c(p) = c(i) 19 tb = b(p) b(p) = b(i) 21 continue c c Start at the beginning of the segment, search for k such that a(k)>t c q = j k = i 20 k = k+1 if(k.gt.q) go to 60 if(a(k).le.ta) go to 20 c c Such an element has now been found now search for a q such that a(q) ¬ î ï í ï ì = Title: Rock Properties Model (RPM3.1) Attachment V Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: V-4 of V-6 Figure V-2 presents the replotted variogram model from VARIOFIT plus the independently computed model using the same parameters, generated using the graphics/spreadsheet package, SIGMAPLOT. Note that the two models appear to overlie one another precisely. Figure V-1. Reproduced screen image of an arbitrary nested variogram model generated using the UNCERT program, VARIOFIT Model Parameters Nest Range C Model 2 5000 25 Spherical 1 1000 10 Spherical Nugget 2 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 0.0 5.0 10.0 15.0 20.0 25.0 30.0 35. 0 40.0 45.0 50.0 Variofi t Test Case X Y Figure V-2. Comparison of model variograms generated by program VARIOFIT and by independent computation using SIGMAPLOT Separation, h 0 1000 2000 3000 4000 5000 6000 Gamma 0 10 20 30 40 50 "variofit" Output Sigma Plot Calculation g = 2 + 10 Sph(1000) + 25 Sph(5000) Title: Rock Properties Model (RPM3.1) Attachment V Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: V-5 of V-6 Conclusion Based on the nearly exact correspondence between the variogram models computed by the two different methods (the maximum difference in gamma is 0.0005 compared to a total sill value of 37.0), we conclude that VARIOFIT is validated for use with the spherical model equation. A similar exercise could be conducted for other variogram model types, if necessary. Listing of SigmaPlot Transform ;transform to compute nested spherical variograms c0 = 2.0 c1 = 10. a1 = 1000 c2 = 25. a2 = 5000 n = 8 ; column in which to start variogram col(n) = data( 0, 7000, 35) ;lag distances to compute variogram for ; intermediate computational results: first nugget, then nest 1, then nest 2 col(n+2) = c0 col(n+3) = if( col(n) < a1 , c1*( 1.5*col(n) / a1 - 0.5*(col(n)^3) / (a1^3) ) , c1 ) col(n+4) = if( col(n) < a2 , c2*( 1.5*col(n) / a2 - 0.5*(col(n)^3) / (a2^3) ) , c2 ) ; output model and preservation of input parameters Output is sum of 3 cols col(n+1) = col(n+2) + col(n+3) + col(n+4) cell(18,1) = c0 cell(18,3) = c1 cell(18,4) = a1 cell(18,6) = c2 cell(18,7) = a2 ; compute difference between two methods col(7) = col(6) - col(9) ; preserve maximum delta gamma cell(18,9) = max(col(7) ) Title: Rock Properties Model (RPM3.1) Attachment V Document Identifier: MDL-NBS-GS-000004 REV 00 ICN 01 Page: V-6 of V-6 INTENTIONALLY LEFT BLANK