CRA-2009 Main | References | CFR Index | Search CRA-2009 | About | View Appendix PA As PDF

 

Title 40 CFR Part 191
Subparts B and C
Compliance Recertification
Application
for the
Waste Isolation Pilot Plant

Appendix PA-2009
Performance Assessment

United States Department of Energy
Waste Isolation Pilot Plant

Carlsbad Field Office
Carlsbad, New Mexico

 


Appendix PA-2009
Performance Assessment


Table of Contents

  PA-1.0  Introduction

  PA-2.0  Overview and Conceptual Structure of the PA

      PA-2.1  Overview of Performance Assessment

          PA-2.1.1  Changes in the CRA-2009 PA

          PA-2.1.2  Conceptual Basis for PA

          PA-2.1.3  Undisturbed Repository Performance

          PA-2.1.4  Disturbed Repository Performance

              PA-2.1.4.1  Cuttings and Cavings

              PA-2.1.4.2  Spallings

              PA-2.1.4.3  Direct Brine Flow

              PA-2.1.4.4  Mobilization of Actinides in Repository Brine

              PA-2.1.4.5  Long-Term Brine Flow up an Intrusion Borehole

              PA-2.1.4.6  Groundwater Flow in the Culebra

              PA-2.1.4.7  Actinide Transport in the Culebra

              PA-2.1.4.8  Intrusion Scenarios

          PA-2.1.5  Compliance Demonstration Method

          PA-2.1.6  Results of the PA

      PA-2.2  Conceptual Structure of the PA

          PA-2.2.1  Regulatory Requirements

          PA-2.2.2  Probabilistic Characterization of Different Futures

          PA-2.2.3  Estimation of Releases

          PA-2.2.4  Probabilistic Characterization of Parameter Uncertainty

      PA-2.3  PA Methodology

          PA-2.3.1  Identification and Screening of FEPs

          PA-2.3.2  Scenario Development and Selection

              PA-2.3.2.1  Undisturbed Repository Performance

              PA-2.3.2.2  Disturbed Repository Performance

                   PA-2.3.2.2.1  Disturbed Repository M Scenario

                  PA-2.3.2.2.2  Disturbed Repository E Scenario

                  PA-2.3.2.2.3  The E2 Scenario

                   PA-2.3.2.2.4  The E1 Scenario

                  PA-2.3.2.2.5  The E1E2 Scenario

              PA-2.3.2.3  Disturbed Repository ME Scenario

              PA-2.3.2.4  Scenarios Retained for Consequence Analysis

          PA-2.3.3  Calculation of Scenario Consequences

  PA-3.0  Probabilistic Characterization of Futures

      PA-3.1  Probability Space

      PA-3.2  AICs and PICs

      PA-3.3  Drilling Intrusion

      PA-3.4  Penetration of Excavated/Nonexcavated Area

      PA-3.5  Drilling Location

      PA-3.6  Penetration of Pressurized Brine

      PA-3.7  Plugging Pattern

      PA-3.8  Activity Level

      PA-3.9  Mining Time

        PA-3.10  Scenarios and Scenario Probabilities

        PA-3.11  CCDF Construction

  PA-4.0  Estimation of Releases

      PA-4.1  Results for Specific Futures

      PA-4.2  Two-Phase Flow:  BRAGFLO

          PA-4.2.1  Mathematical Description

          PA-4.2.2  Initial Conditions

          PA-4.2.3  Creep Closure of Repository

          PA-4.2.4  Fracturing of MBs and DRZ

          PA-4.2.5  Gas Generation

          PA-4.2.6  Capillary Action in the Waste

          PA-4.2.7  Shaft Treatment

          PA-4.2.8  Option D Panel Closures

              PA-4.2.8.1  Panel Closure Concrete

              PA-4.2.8.2  Panel Closure Abutment with MBs

              PA-4.2.8.3  DRZ Above the Panel Closure

              PA-4.2.8.4  Empty Drift and Explosion Wall Materials

          PA-4.2.9  Borehole Model

            PA-4.2.10  Castile Brine Reservoir

            PA-4.2.11  Numerical Solution

            PA-4.2.12  Gas and Brine Flow across Specified Boundaries

            PA-4.2.13  Additional Information

      PA-4.3  Radionuclide Transport in the Salado: NUTS

          PA-4.3.1  Mathematical Description

          PA-4.3.2  Calculation of Maximum Concentration ST( Br, Ox, El)

          PA-4.3.3  Radionuclides Transported

          PA-4.3.4  NUTS Tracer Calculations

          PA-4.3.5  NUTS Transport Calculations

          PA-4.3.6  Numerical Solution

          PA-4.3.7  Additional Information

      PA-4.4  Radionuclide Transport in the Salado:  PANEL

          PA-4.4.1  Mathematical Description

          PA-4.4.2  Numerical Solution

          PA-4.4.3  Implementation in PA

          PA-4.4.4  Additional Information

      PA-4.5  Cuttings and Cavings to Surface:  CUTTINGS_S

          PA-4.5.1  Cuttings

          PA-4.5.2  Cavings

              PA-4.5.2.1  Laminar Flow Model

              PA-4.5.2.2  Turbulent Flow Model

              PA-4.5.2.3  Calculation of Rf

          PA-4.5.3  Additional Information

      PA-4.6  Spallings to Surface: DRSPALL and CUTTINGS_S

          PA-4.6.1  Summary of Assumptions

          PA-4.6.2  Conceptual Model

              PA-4.6.2.1  Wellbore Flow Model

                  PA-4.6.2.1.1  Wellbore Initial Conditions

                  PA-4.6.2.1.2  Wellbore Boundary Conditions

              PA-4.6.2.2  Repository Flow Model

              PA-4.6.2.3  Wellbore to Repository Coupling

                  PA-4.6.2.3.1  Flow Prior to Penetration

                   PA-4.6.2.3.2  Flow After Penetration

                  PA-4.6.2.3.3  Cavity Volume After Penetration

                  PA-4.6.2.3.4  Waste Failure

                   PA-4.6.2.3.5  Waste Fluidization

           PA-4.6.3  Numerical Model

               PA-4.6.3.1  Numerical Method—Wellbore

              PA-4.6.3.2  Numerical Method—Repository

              PA-4.6.3.3  Numerical Method—Wellbore to Repository Coupling

          PA-4.6.4  Implementation

          PA-4.6.5  Additional Information

      PA-4.7  DBR to Surface:  BRAGFLO

          PA-4.7.1  Overview of Conceptual Model

          PA-4.7.2  Linkage to Two-Phase Flow Calculation

          PA-4.7.3  Conceptual Representation for Flow Rate rDBR(t)

          PA-4.7.4  Determination of Productivity Index Jp

          PA-4.7.5  Determination of Waste Panel Pressure pw( t) and DBR

          PA-4.7.6  Boundary Value Pressure pwf

          PA-4.7.7  Boundary Value Pressure pwE1

              PA-4.7.7.1  Solution for Open Borehole

              PA-4.7.7.2  Solution for Sand-Filled Borehole

          PA-4.7.8  End of DBR

          PA-4.7.9  Numerical Solution

            PA-4.7.10  Additional Information

      PA-4.8  Groundwater Flow in the Culebra Dolomite

          PA-4.8.1  Mathematical Description

          PA-4.8.2  Implementation

          PA-4.8.3  Computational Grids and Boundary Value Conditions

          PA-4.8.4  Numerical Solution

          PA-4.8.5  Additional Information

      PA-4.9  Radionuclide Transport in the Culebra Dolomite

          PA-4.9.1  Mathematical Description

              PA-4.9.1.1  Advective Transport in Fractures

              PA-4.9.1.2  Diffusive Transport in the Matrix

              PA-4.9.1.3Coupling Between Fracture and Matrix Equations

              PA-4.9.1.4  Source Term

              PA-4.9.1.5  Cumulative Releases

          PA-4.9.2  Numerical Solution

              PA-4.9.2.1  Discretization of Fracture Domain

              PA-4.9.2.2  Discretization of Matrix Equation

              PA-4.9.2.3  Fracture-Matrix Coupling

              PA-4.9.2.4  Cumulative Releases

          PA-4.9.3  Additional Information

  PA-5.0  Probabilistic Characterization of Subjective Uncertainty

      PA-5.1  Probability Space

      PA-5.2  Variables Included for Subjective Uncertainty

      PA-5.3  Variable Distributions

      PA-5.4  Correlations and Dependencies

      PA-5.5  Separation of Aleatory and Epistemic Uncertainty

  PA-6.0  Computational Procedures

      PA-6.1  Sampling Procedures

      PA-6.2  Sample Size for Incorporation of Subjective Uncertainty

      PA-6.3  Statistical Confidence on Mean CCDF

      PA-6.4 Generation of Latin Hypercube Samples

      PA-6.5  Generation of Individual Futures

      PA-6.6  Construction of CCDFs

      PA-6.7  Mechanistic Calculations

          PA-6.7.1  BRAGFLO Calculations

          PA-6.7.2  NUTS Calculations

          PA-6.7.3  PANEL Calculations

          PA-6.7.4  DRSPALL Calculations

          PA-6.7.5  CUTTINGS_S Calculations

          PA-6.7.6  BRAGFLO Calculations for DBR Volumes

          PA-6.7.7  MODFLOW Calculations

          PA-6.7.8  SECOTP2D Calculations

      PA-6.8Computation of Releases

          PA-6.8.1  Undisturbed Releases

          PA-6.8.2  Direct Releases

              PA-6.8.2.1  Construction of Cuttings and Cavings Releases

              PA-6.8.2.2  Construction of Spallings Releases

              PA-6.8.2.3  Construction of DBRs

          PA-6.8.3  Radionuclide Transport Through the Culebra

          PA-6.8.4  Determining Initial Conditions for Direct and Transport Releases

              PA-6.8.4.1  Determining Repository and Panel Conditions

              PA-6.8.4.2  Determining Distance from Previous Intrusions

          PA-6.8.5  CCDF Construction

      PA-6.9  Sensitivity Analysis

          PA-6.9.1  Scatterplots

          PA-6.9.2  Regression Analysis

          PA-6.9.3  Stepwise Regression Analysis

  PA-7.0  Results for the Undisturbed Repository

      PA-7.1  Salado Flow

          PA-7.1.1  Pressure in the Repository

          PA-7.1.2  Brine Saturation in the Waste

          PA-7.1.3  Brine Flow Out of the Repository

      PA-7.2  Radionuclide Transport

          PA-7.2.1  Radionuclide Transport to the Culebra

          PA-7.2.2  Radionuclide Transport to the Land Withdrawal Boundary

  PA-8.0  Results for a Disturbed Repository

      PA-8.1  Drilling Scenarios

      PA-8.2  Mining Scenarios

      PA-8.3  Salado Flow

          PA-8.3.1  Pressure in the Repository

          PA-8.3.2  Brine Saturation

          PA-8.3.3  Brine Flow out of the Repository

      PA-8.4  Radionuclide Transport

          PA-8.4.1  Radionuclide Source Term

          PA-8.4.2  Transport through MBs and Shaft

          PA-8.4.3  Transport to the Culebra

              PA-8.4.3.1  Single Intrusion Scenarios

              PA-8.4.3.2  Multiple Intrusion Scenario

          PA-8.4.4  Transport Through the Culebra

              PA-8.4.4.1  Partial Mining Results

              PA-8.4.4.2  Full Mining Results

      PA-8.5  Direct Releases

          PA-8.5.1  Cuttings and Cavings

          PA-8.5.2  Spallings

              PA-8.5.2.1  DRSPALL Results

              PA-8.5.2.2  CUTTINGS_S Results

          PA-8.5.3  DBRs

  PA-9.0  Normalized Releases

      PA-9.1  Cuttings and Cavings

      PA-9.2  Spallings

      PA-9.3  Direct Brine

      PA-9.4  Groundwater Transport

      PA-9.5  Total Normalized Releases

    PA-10.0  References

List of Figures

Figure PA-1.      Overall Mean CCDFs for Total Normalized Releases:  CRA-2009 PA and CRA-2004 PABC

Figure PA-2.  Computational Models Used in PA

Figure PA-3.  Construction of the CCDF Specified in 40 CFR Part 191 Subpart B

Figure PA-4.      Distribution of CCDFs Resulting from Possible Values for the Sampled Parameters

Figure PA-5.      Example Mean Probability of Release and the Confidence Limits on the Mean from CRA-2009 PA

Figure PA-6.  Example CCDF Distribution From CRA-2009 PA (Replicate 1)

Figure PA-7.  Logic Diagram for Scenario Analysis

Figure PA-8.   Conceptual Release Pathways for the UP Scenario

Figure PA-9.  Conceptual Release Pathways for the Disturbed Repository M Scenario

Figure PA-10.    Conceptual Release Pathways for the Disturbed Repository Deep Drilling E2 Scenario

Figure PA-11.    Conceptual Release Pathways for the Disturbed Repository Deep Drilling E1 Scenario

Figure PA-12.    Conceptual Release Pathways for the Disturbed Repository Deep Drilling E1E2 Scenario

Figure PA-13.  CDF for Time Between Drilling Intrusions

Figure PA-14.  Discretized Locations for Drilling Intrusions

Figure PA-15.  Computational Grid Used in BRAGFLO for PA

Figure PA-16.  Definition of Element Depth in BRAGFLO Grid

Figure PA-17.  Identification of Individual Cells in BRAGFLO Grid

Figure PA-18.    Schematic View of the Simplified Shaft Model (numbers on right indicate dimensions in meters)

Figure PA-19.  Schematic Side View of Option D Panel Closure

Figure PA-20.  Representation of Option D Panel Closures in the BRAGFLO Grid

Figure PA-21. Selecting Radionuclides for the Release Pathways Conceptualized by PA

Figure PA-22.  Detail of Rotary Drill String Adjacent to Drill Bit

Figure PA-23.  Schematic Diagram of the Flow Geometry Prior to Repository Penetration

Figure PA-24.  Schematic Diagram of the Flow Geometry After Repository Penetration

Figure PA-25.  Effective Wellbore Flow Geometry Before Bit Penetration

Figure PA-26.  Effective Wellbore Flow Geometry After Bit Penetration

Figure PA-27.  Finite-Difference Zoning for Wellbore

Figure PA-28.  DBR Grid Used in PA

Figure PA-29.  Assignment of Initial Conditions for DBR Calculation

Figure PA-30.  Borehole Representation Used for Poettmann-Carpenter Correlation

Figure PA-31.  Areas of Potash Mining in the McNutt Potash Zone

Figure PA-32.    Modeling Domain for Groundwater Flow (MODFLOW) and Radionuclide Transport (SECOTP2D) in the Culebra

Figure PA-33.  Boundary Conditions Used for Simulations of Brine Flow in the Culebra

Figure PA-34.    Finite-Difference Grid Showing Cell Index Numbering Convention Used by MODFLOW

Figure PA-35.  Parallel-Plate, Dual-Porosity Conceptualization

Figure PA-36.    Schematic of Finite-Volume Staggered Mesh Showing Internal and Ghost Cells

Figure PA-37.  Illustration of Stretched Grid Used for Matrix Domain Discretization

Figure PA-38.    Correlation Between S_HALITE:PRMX_LOG and S_HALITE:COMP_RCK

Figure PA-39.    Correlation between CASTILLER:PRMX_LOG and CASTILLER:COMP_RCK

Figure PA-40.    Dependency between WAS_AREA:GRATMICI and WAS_AREA:GRATMICH

Figure PA-41.  Logic Diagram for Determining the Intrusion Type

Figure PA-42.  Processing of Input Data to Produce CCDFs

Figure PA-43.    Pressure in the Waste Panel Region, Replicate R1, Scenario S1, CRA-2009 PA

Figure PA-44.    Primary Correlations of Pressure in the Waste Panel Region with Uncertain Parameters, Replicate R1, Scenario S1, CRA-2009 PA

Figure PA-45.    Brine Saturation in the Waste Panel Region, Replicate R1, Scenario S1, CRA-2009 PA

Figure PA-46.    Primary Correlations of Brine Saturation in the Waste Panel Region with Uncertain Parameters, Replicate R1, Scenario S1, CRA-2009 PA

Figure PA-47.    Brine Flow Away From the Repository, Replicate R1, Scenario S1, CRA-2009 PA

Figure PA-48.    Brine Flow via All MBs Across the LWB, Replicate R1, Scenario S1, CRA-2009 PA

Figure PA-49.    Primary Correlations of Total Cumulative Brine Flow Away from the Repository with Uncertain Parameters, Replicate R1, Scenario S1, CRA-2009 PA

Figure PA-50.    Pressure in the Waste Panel Region for Replicate R1, Scenario S2, CRA-2009 PA

Figure PA-51.    Pressure in the Waste Panel Region for Replicate R1, Scenario S4, CRA-2009 PA

Figure PA-52.    Primary Correlations of Pressure in the Waste Panel Region with Uncertain Parameters, Replicate R1, Scenario S2, CRA-2009 PA

Figure PA-53.    Primary Correlations of Pressure in the Waste Panel Region with Uncertain Parameters, Replicate R1, Scenario S4, CRA-2009 PA

Figure PA-54.    Brine Saturation in the Waste Panel Region for Replicate R1, Scenario S2, CRA-2009 PA

Figure PA-55.    Brine Saturation in the Waste Panel Region for Replicate R1, Scenario S4, CRA-2009 PA

Figure PA-56.    Primary Correlations of Brine Saturation in the Waste Panel Region with Uncertain Parameters, Replicate R1, Scenario S2, CRA-2009 PA

Figure PA-57.    Primary Correlations of Brine Saturation in the Waste Panel Region with Uncertain Parameters, Replicate R1, Scenario S4, CRA-2009 PA

Figure PA-58.    Total Cumulative Brine Outflow in Replicate R1, Scenario S2, CRA-2009 PA

Figure PA-59.    Brine Flow via All MBs Across the LWB in Replicate R1, Scenario S2, CRA-2009 PA

Figure PA-60.    Total Cumulative Brine Outflow in Replicate R1, Scenario S4, CRA-2009 PA

Figure PA-61.    Brine Flow via All MBs Across the LWB in Replicate R1, Scenario S4, CRA-2009 PA

Figure PA-62.    Primary Correlations of Cumulative Brine Flow Away from the Repository with Uncertain Parameters, Replicate R1, Scenario S2, CRA-2009 PA

Figure PA-63.    Primary Correlations of Cumulative Brine Flow Away from the Repository with Uncertain Parameters, Replicate R1, Scenario S4, CRA-2009 PA

Figure PA-64.    Total Mobilized Concentrations in Salado Brine, Replicate R1, CRA-2009 PA

Figure PA-65.    Total Mobilized Concentrations in Castile Brine, Replicate R1, CRA-2009 PA

Figure PA-66.  Cumulative Normalized Release to the Culebra, Scenario S2, CRA-2009 PA

Figure PA-67.  Cumulative Normalized Release to the Culebra, Scenario S3, CRA-2009 PA

Figure PA-68.  Cumulative Normalized Release to the Culebra, Scenario S4, CRA-2009 PA

Figure PA-69.  Cumulative Normalized Release to the Culebra, Scenario S5, CRA-2009 PA

Figure PA-70.    Cumulative Normalized Release to the Culebra, Replicate R1, Scenario S6, CRA-2009 PA

Figure PA-71.    Scatterplot of Waste Particle Diameter Versus Spallings Volume, CRA-2009 PA

Figure PA-72. Scatterplot of Waste Permeability Versus Spallings Volume, CRA-2009 PA

Figure PA-73.    Sensitivity of DBR Volumes to Pressure and Mobile Brine Saturation, Replicate R1, Scenario S2, Lower Panel, CRA-2009 PA.  (Symbols Indicate the Range of Mobile Brine Saturation Given in the Legend.)

Figure PA-74.     Overall Mean CCDFs for Cuttings and Cavings Releases: CRA-2009 PA and CRA-2004 PABC

Figure PA-75.     Overall Mean CCDFs for Spallings Releases:  CRA-2009 PA and CRA 2004 PABC

Figure PA-76.  Overall Mean CCDFs for DBRs: CRA-2009 PA and CRA-2004 PABC

Figure PA-77.    Mean CCDFs for Releases from the Culebra for Replicate R2: CRA-2009 PA and CRA-2004 PABC

Figure PA-78.  The Preponderance and Distribution of Zeros Can Control the Regression

Figure PA-79.  Total Normalized Releases, Replicates R1, R2, and R3, CRA-2009 PA

Figure PA-80.    Confidence Interval on Overall Mean CCDF for Total Normalized Releases, CRA-2009 PA

Figure PA-81.    Mean CCDFs for Components of Total Normalized Releases, Replicate R1, CRA-2009 PA

Figure PA-82.    Mean CCDFs for Components of Total Normalized Releases, Replicate R2, CRA-2009 PA

Figure PA-83.    Mean CCDFs for Components of Total Normalized Releases, Replicate R3, CRA-2009 PA

Figure PA-84.    Overall Mean CCDFs for Total Normalized Releases: CRA-2009 PA and CRA-2004 PABC

List of Tables

Table PA-1.  WIPP Project Changes and Cross References

Table PA-2.       Release Limits for the Containment Requirements ( U.S. Environmental Protection Agency 1985, Appendix A, Table 1)

Table PA-3.  Parameter Values Used in Representation of Two-Phase Flow

Table PA-4.  Models for Relative Permeability and Capillary Pressure in Two-Phase Flow

Table PA-5.  Initial Conditions in the Rustler

Table PA-6.    Probabilities for Biodegradation of Different Organic Materials (WAS_AREA:PROBDEG) in the CRA-2009 PA and the CRA-2004 PA

Table PA-7.  Permeabilities for Drilling Intrusions Through the Repository

Table PA-8.  Boundary Value Conditions for Pg and Pb

Table PA-9.  Auxiliary Dirichlet Conditions for Sg and Pb

Table PA-10.     Calculated Values for Dissolved Solubility (see Appendix SOTERM-2009, Table SOTERM-16)

Table PA-11.     Scale Factor SFHum( Br, Ox, El) Used in Definition of SHum( Br, Ox, El) (see Appendix SOTERM-2009, Table SOTERM-21)

Table PA-12.     Scale Factor SFMic( Ox, El) and Upper Bound UBMic( Ox, El) (mol/L) Used in Definition of SMic( Br, Ox, El) (see Appendix SOTERM-2009, Table SOTERM-21)

Table PA-13.  Combination of Radionuclides for Transport

Table PA-14.  Initial and Boundary Conditions for Cbl( x, y, t) and Csl( x, y, t)

Table PA-15.  Uncertain Parameters in the DRSPALL Calculations

Table PA-16.  Initial Porosity in the DBR Calculation

Table PA-17.  Boundary Conditions for pb and Sg in DBR Calculations

Table PA-18.  Radionuclide Culebra Transport Diffusion Coefficients

Table PA-19.  Variables Representing Epistemic Uncertainty in the CRA-2009 PA

Table PA-20.  Sampled Parameters Added Since the CRA-2004 PA

Table PA-21.  Sampled Parameters Removed Since the CRA-2004 PA

Table PA-22.     Correlation Observed Between Variables S_HALITE:PRMX_LOG and S_HALITE:COMP_RCK in Replicate 1.  A Value of -0.99 was Specified.

Table PA-23.     Correlation Observed Between Variables CASTILER:PRMX_LOG and CASTILER:COMP_RCK in Replicate 1.  A Value of -0.75 was Specified.

Table PA-24.  Algorithm to Generate a Single Future

Table PA-25.  BRAGFLO Scenarios in the CRA-2009 PA

Table PA-26.  NUTS Release Calculations in the CRA-2009 PA

Table PA-27.  CUTTINGS_S Release Calculations in the CRA-2009 PA

Table PA-28.  MODFLOW Scenarios in the CRA-2009 PA

Table PA-29.  SECOTP2D Scenarios in the CRA-2004 PA

Table PA-30.     Number of Realizations with Radionuclide Transport to the LWB Under Partial-Mining Conditionsa,b

Table PA-31.     Number of Realizations with Radionuclide Transport to the LWB Under Full-Mining Conditionsa,b

Table PA-32.  CRA-2009 PA Cuttings and Cavings Area Statistics

Table PA-33.  CRA-2009 PA Spallings Volume Statistics

Table PA-34.  CRA-2009 PA DBR Volume Statistics

Table PA-35.     CRA-2009 PA and CRA-2004 PABCa Statistics on the Overall Mean for Total Normalized Releases in EPA Units at Probabilities of 0.1 and 0.001

 

 

 

Acronyms and Abbreviations

%                                          percent

AIC                                      active institutional control

An                                        actinide

BRAGFLO                           BRine And Gas FLOw computer code

C                                          Celsius

CCA                                     Compliance Certification Application

CCDF                                   complementary cumulative distribution function

CDF                                     cumulative distribution function

CFR                                      Code of Federal Regulations

CH-TRU                               contact-handled transuranic

Ci                                         curies

CL                                        Confidence Limit

CPR                                      cellulosic, plastic, and rubber

CRA                                     Compliance Recertification Application

DBR                                     direct brine release

DDZ                                     drilling damaged zone

DOE                                     U.S. Department of Energy

DP                                        disturbed repository performance

DRZ                                     disturbed rock zone

E                                           deep drilling scenario

EPA                                      U.S. Environmental Protection Agency

ERDA                                  U.S.Energy Research and Development Administration

FEP                                      feature, event, and process

FMT                                     Fracture-Matrix Transport

FVW                                     fraction of excavated repository volume occupied by waste

gal                                        gallon

GWB                                    Generic Weep Brine

in                                          inch

J                                           Joule

K                                          Kelvin

Kd                                         distribution coefficient

kg                                         kilogram

km                                        kilometer

km2                                       square kilometers

L                                           liter

LHS                                      Latin hypercube sampling

LWB                                     land withdrawal boundary

M                                          mining scenario

m                                          meter

m2                                         square meters

m3                                         cubic meters

MB                                       marker bed

ME                                       mining and drilling scenario

mi                                         miles

mol                                       mole

MPa                                      megapascal

MTHM                                 metric tons of heavy metal

MWd                                    megawatt-days

N                                          Newton

Pa                                         Pascal

PA                                        performance assessment

PABC                                   performance assessment baseline calculation

PAVT                                   Performance Assessment Verification Test

PCC                                      partial correlation coefficient

PCS                                      panel closure system

PDE                                      partial differential equation

PDF                                      probability distribution function

PIC                                       passive institutional control

PRCC                                   partial rank correlation coefficient

RH-TRU                               remote-handled transuranic

RKS                                     Redlich-Kwong-Soave

RoR                                      Rest of Repository

s                                           second

s2                                          seconds squared

SCF/d                                   standard cubic feet per day

SMC                                     Salado Mass Concrete

SNL                                      Sandia National Laboratories

SRC                                      standardized regression coefficient

T field                                  transmissivity field

TRU                                     transuranic

TVD                                     Total Variation Diminishing

UP                                        undisturbed repository performance

WIPP                                    Waste Isolation Pilot Plant

yr                                          year

 

Elements and Chemical Compounds

Al                                         aluminum

Am                                       americium

C                                          carbon

C6H10O5                                generic formula for CPR

Ca                                         calcium

CH4                                      methane

Cm                                        curium

CO2                                      carbon dioxide

Cr                                         chromium

Cs                                         cesium

Fe                                         iron

H2                                         hydrogen gas

H2O                                      water

H2S                                       hydrogen sulfide

I                                            iodine

Mg                                        magnesium

Mg(OH)2                                                  brucite

Mg5(CO3)4(OH)2×4H2O         hydromagnesite (5424)

MgO                                     magnesium oxide, or periclase

Mn                                        manganese

Ni                                         nickel

NO3 -                                     nitrate

Np                                        neptunium

Pb                                         lead

Pm                                        promethium

Pu                                         plutonium

Ra                                         radium

Sn                                         tin

SO4                                       sulfate

SO42-                                     sulfate ion

Sr                                         strontium

Tc                                         technetium

Th                                         thorium

U                                          uranium

V                                          vanadium

 



This appendix presents the mathematical models used to evaluate performance of the Waste Isolation Pilot Plant (WIPP) disposal system and the results of these models for the 2009 Compliance Recertification Application (CRA-2009) Performance Assessment (PA).  The term PA signifies an analysis that (1) identifies the processes and events that might affect the disposal system; (2) examines the effects of these processes and events on the performance of the disposal system; and (3) estimates the cumulative releases of radionuclides, considering the associated uncertainties, caused by all significant processes and events (40 CFR § 191.12 [U.S. Environmental Protection Agency 1993]).  PA is designed to address three primary questions about the WIPP:

Q1:        What processes and events that might affect the disposal system could take place at the WIPP site over the next 10,000 years?

Q2:        How likely are the various processes and events that might affect the disposal system to take place at the WIPP site over the next 10,000 years?

Q3:        What are the consequences of the various processes and events that might affect the disposal system that could take place at the WIPP site over the next 10,000 years?

In addition, accounting for uncertainty in the parameters of the PA models leads to a further question:

Q4:        How much confidence should be placed in answers to the first three questions?

These questions give rise to a methodology for quantifying the probability distribution of possible radionuclide releases from the WIPP repository over the next 10,000 years and characterizing the uncertainty in that distribution due to imperfect knowledge about the parameters contained in the models used to predict releases. The containment requirements of 40 CFR § 191.13 require this probabilistic methodology.

This appendix is organized as follows: Section PA-2.0 gives an overview and describes the overall conceptual structure of the CRA-2009 PA.  The WIPP PA is designed to address the requirements of section 191.13, and thus involves three basic entities:  (1) a probabilistic characterization of different futures that could occur at the WIPP site over the next 10,000 years, (2) models for both the physical processes that take place at the WIPP site and the estimation of potential radionuclide releases that may be associated with these processes, and (3) a probabilistic characterization of the uncertainty in the models and parameters that underlies the WIPP PA.  Section PA-2.0 is supplemented by Appendix SCR-2009, which documents the results of the screening process for features, events, and processes (FEPs) that are retained in the conceptual models of repository performance.

Section PA-3.0 describes the probabilistic characterization of different futures.  This characterization plays an important role in the construction of the complementary cumulative distribution function (CCDF) specified in section 191.13.  Regulatory guidance and extensive review of the WIPP site identified exploratory drilling for natural resources and the mining of potash as the only significant disruptions at the WIPP site with the potential to affect radionuclide releases to the accessible environment.  Section PA-3.0 summarizes the stochastic variables that represent future drilling and mining events in the PA.  The results of the PA for CRA-2009, as documented in Section PA-7.0, Section PA-8.0, and Section PA-9.0 of this appendix, confirm that direct releases from drilling intrusions are the major contributors to radionuclide releases to the accessible environment.

Section PA-4.0 presents the mathematical models for both the physical processes that take place at the WIPP and the estimation of potential radionuclide releases. The mathematical models implement the conceptual models as prescribed in 40 CFR § 194.23 (2004), and permit the construction of the CCDF specified in section 191.13.  Models presented in Section PA-4.0 include two-phase (i.e., gas and brine) flow in the vicinity of the repository; radionuclide transport in the Salado Formation (hereafter referred to as the Salado); releases to the surface at the time of a drilling intrusion due to cuttings, cavings, spallings, and direct brine releases (DBRs); brine flow in the Culebra Dolomite Member of the Rustler Formation (hereafter referred to as the Culebra); and radionuclide transport in the Culebra.  Section PA-4.0 is supplemented by Appendices MASS-2009, TFIELD-2009, and PORSURF-2009.  Appendix MASS-2009 discusses the modeling assumptions used in the WIPP PA.  Appendix TFIELD-2009 discusses the generation of the transmissivity fields (T fields) used to model groundwater flow in the Culebra.  Appendix PORSURF-2009 presents results from modeling the effects of excavated region closure, waste consolidation, and gas generation in the repository.

Section PA-5.0 discusses the probabilistic characterization of parameter uncertainty, and summarizes the uncertain variables incorporated into the CRA-2009 PA, the distributions assigned to these variables, and the correlations between variables.  Section PA-5.0 is supplemented by Fox (2008) and Appendix SOTERM-2009.  Fox (2008) catalogs the full set of parameters used in the CRA-2009 PA, previously referred to as the CCA Appendix PAR and the CRA-2004, Appendix PA, Attachment PAR.  Appendix SOTERM-2009 describes the actinide (An) source term for the WIPP performance calculations, including the mobile concentrations of actinides that may be released from the repository in brine.

Section PA-6.0 summarizes the computational procedures used in the CRA-2009 PA, including sampling techniques (i.e., random and Latin hypercube sampling (LHS)); sample size, statistical confidence for mean CCDF, generation of sample, generation of individual futures, construction of CCDFs, calculations performed with the models discussed in Section PA-4.0, construction of releases for each future, and the sensitivity analysis techniques in use.

Section PA-7.0 presents the results of the PA for an undisturbed repository.  Releases from the undisturbed repository are determined by radionuclide transport in brine flowing from the repository to the land withdrawal boundary (LWB) through the marker beds (MBs) or shafts.  Releases in the undisturbed scenario are used to demonstrate compliance with the individual and groundwater protection requirements in 40 CFR Part 191 (40 CFR § 194.51 and 40 CFR § 194.52).

Section PA-8.0 presents PA results for a disturbed repository.  As will be discussed in Section PA-2.3.1, the only future events and processes in the analysis of disturbed repository performance are those associated with mining and deep drilling.  Release mechanisms include direct releases at the time of the intrusion via cuttings, cavings, spallings, and DBR, and long-term releases via radionuclide transport up abandoned boreholes to the Culebra and thence to the LWB.

Section PA-9.0 presents the set of CCDFs resulting from the CRA-2009 PA.  This material supplements 40 CFR § 194.34, which demonstrates compliance with the containment requirements of section 191.13.  Section PA-9.0 presents the most significant output variables from the PA models, accompanied by sensitivity analyses to determine which subjectively uncertain parameters are most influential in the uncertainty of PA results.

The overall structure of the CRA-2009 PA does not differ from that presented in the first WIPP Compliance Certification Application (CCA) (U.S. Department of Energy 1996), the CRA-2004 (U.S. Department of Energy 2004) or the CRA-2004 Performance Assessment Baseline Calculation (PABC) (Leigh et al. 2005).  This recertification application appendix follows the approach used by Helton et al. (1998) to document the mathematical models used in the CCA PA and the results of that analysis.  Much of the content of this appendix derives from Helton et al. (1998); these authors’ contributions are gratefully acknowledged.


Because of the amount and complexity of the material presented in Appendix PA-2009, an introductory summary is provided below, followed by detailed discussions of the topics in the remainder of this section, which is organized as follows:

Section PA-2.1 – Overview of PA and the results

Section PA-2.2 – The conceptual structure of the PA used to evaluate compliance with the containment requirements

Section PA-2.3 – The overall methodology used to develop FEPs, the screening methodology applied to the FEPs, the results of the screening process, and the development of the scenarios considered in the system-level consequence analysis

The U.S. Department of Energy (DOE) continues to use the same PA methodology as in the CCA and CRA-2004 because changes that have been made since the U.S. Environmental Protection Agency (EPA) certified WIPP do not impact PA methodology.  A corresponding detailed presentation for the CRA-2004 PA methodology is provided in the CRA-2004, Chapter 6.0, Section 6.1, and a detailed presentation for the CRA-2004 PABC implementation is provided in Leigh et al. (2005).  A corresponding detailed presentation for the CCA PA methodology is provided in Helton et al. (1998, Section 2).

A demonstration of future repository performance was required by the disposal standards in Part 191.  The EPA required a PA to demonstrate that potential cumulative releases of radionuclides to the accessible environment over a 10,000-year period after disposal are less than specified limits based on the nature of the materials disposed (section 191.13).  The PA is to determine the effects of all significant processes and events that may affect the disposal system, consider the associated uncertainties of the processes and events, and estimate the probable cumulative releases of radionuclides.

A PA was included in the CCA.  This was the first demonstration of compliance by the DOE with the EPA’s disposal standards.  The EPA required a verification PA based on the CCA PA that included revised parameters and distributions.  This PA was termed the CCA Performance Assessment Verification Test (PAVT) (Trovoto 1997).  The EPA based the original certification of the WIPP on the information in the CCA and the results of the CCA PAVT (U.S. Environmental Protection Agency 1998a).  The CCA PAVT is documented in Sandia National Laboratories (SNL) 1997 and U.S. Department of Energy 1997.

The WIPP is required to be recertified every five years after first waste receipt (Public Law 02-579).  A revised PA was included in the first recertification application in 2004.  This PA is termed the CRA-2004 PA, and is documented in the CRA-2004, Chapter 6.0.  The EPA again required a verification PA using revised modeling assumptions and parameters (Cotsworth 2005).  This PA was termed the CRA-2004 PABC and was documented in Leigh et al. (2005).

As part of the five-year recertification cycle, a PA is included in the CRA-2009.  The CRA-2009 PA is a culmination of the previous PAs and is not significantly different from the CRA-2004 PABC, as the methodologies, conceptual models, and assumptions have not changed.  Updates to parameters and improvements to computer codes are incorporated in the CRA-2009.

A list of changes to PA since the CRA-2004 and citations for where they are discussed is shown in Table PA-1. In addition to the changes discussed in Table PA-1, the terminology used to describe uncertainty has been updated to reflect the usage now prevalent in the risk assessment literature. Previously, uncertainty in model parameters was referred to as “subjective uncertainty,” and that due to stochastic processes was referred to as “stochastic uncertainty.” In the years since these terms were first employed, these concepts have matured, and the term “epistemic uncertainty” is now used to describe uncertainty from lack of knowledge, while the term “aleatory uncertainty” now describes uncertainty due to natural variability, e.g. uncertainty arising from future events whose occurrence can be defined in terms of probabilities. In this text, the terms epistemic uncertainty and aleatory uncertainty are used in place of the terms subjective uncertainty and stochastic uncertainty, respectively.

Table PA-1.  WIPP Project Changes and Cross References

WIPP Project Change

Cross Reference

CRA-2004 to CRA-2004 PABC Changesa

Inventory Information

Leigh, Trone, and Fox 2005

An Solubility Values

Garner and Leigh 2005

An Solubility Uncertainty Ranges

Garner and Leigh 2005

Microbial Gas Generation Model

Nemer and Stein 2005

Culebra T Fields Mining Modification

Lowry and Kanney 2005

Anhydrite Material Parameters

Vugrin et al. 2005

SPALL Model Parameters

Vugrin 2005

CRA-2004 PABC to CRA-2009 Changesb

DBR Maximum Duration Parameter

Kirkes 2007

Conditional Relationship Between Humid and Inundated Cellulosic, Plastic, or Rubber (CPR) Degradation Rates

Kirchner 2008a

BRAGFLO Code Improvements

Nemer and Clayton 2008

Capillary Pressure and Relative Permeability Model

Nemer and Clayton 2008

Drilling Rate

Clayton 2008a

Parameter Error Corrections

·    Emplaced CPR

·    Halite/Disturbed Rock Zone (DRZ) Parameter

·    Fraction of Repository Occupied by Waste

·    NUTS and DBR Calculation Input Files

Nemer 2007a, Ismail 2007a, Dunagan 2007, Ismail 2007b, Clayton 2007

a     See Leigh et al. 2005 for additional discussion of these changes and their implementation in the CRA-2004 PABC.

b     See Clayton et al. 2008 for additional discussions of these changes and their implementation in the CRA-2009 PA.

 

From this assessment, the DOE has demonstrated that the WIPP continues to comply with the containment requirements of section 191.13.  The containment requirements are stringent and state that the DOE must demonstrate a reasonable expectation that the probabilities of cumulative radionuclide releases from the disposal system during the 10,000 years following closure will fall below specified limits.  The PA analyses supporting this determination must be quantitative and consider uncertainties caused by all significant processes and events that may affect the disposal system, including future inadvertent human intrusion into the repository.  A quantitative PA is conducted using a series of loosely coupled computer models in which epistemic parameter uncertainties are addressed by a stratified Monte Carlo sampling procedure on selected input parameters, and uncertainties related to future intrusion events are addressed using simple random sampling.

As required by regulation, results of the PA are displayed as CCDFs showing the probability that cumulative radionuclide releases from the disposal system will exceed the values calculated for scenarios considered in the analysis.  These CCDFs are calculated using reasonable and, in some cases, conservative conceptual models based on the scientific understanding of the disposal system’s behavior.  Parameters used in these models are derived from experimental data, field observations, and relevant technical literature.  Changes to the CCA and CRA-2004 parameters and models since the original certification have been incorporated into the CRA-2009 PA (Clayton 2008a).  The overall mean CCDF continues to lie entirely below the specified limits, and the WIPP therefore continues to be in compliance with the containment requirements of 40 CFR Part 191 Subpart B (see Section PA-2.1.6).  Sensitivity analysis of results shows that the location of the mean CCDF is dominated by radionuclide releases that could occur on the surface during an inadvertent penetration of the repository by a future drilling operation (Section PA-9.0).  Releases of radionuclides to the accessible environment from transport in groundwater through the shaft seal systems and the subsurface geology are negligible, with or without human intrusion, and do not significantly contribute to the mean CCDF.  No releases are predicted to occur at the ground surface in the absence of human intrusion. The natural and engineered barrier systems of the WIPP provide robust and effective containment of transuranic (TRU) waste, even if the repository is penetrated by multiple boreholes.

The foundations of PA are a thorough understanding of the disposal system and the possible future interactions of the repository, waste, and surrounding geology.  The DOE’s confidence in the results of PA is based in part on the strength of the original research done during site characterization, experimental results used to develop and confirm parameters and models, and robustness of the facility design.

The progression of compliance applications document these aspects of PA leading up to the CRA-2009 PA (i.e., the CCA, the CCA PAVT [Sandia National Laboratories 1997 and U.S. Department of Energy 1997], the CRA-2004 PA, and the CRA-2004 PABC [Leigh et al. 2005]).

The interactions of the repository and waste with the geologic system, and the response of the disposal system to possible future inadvertent human intrusion, are described in Section PA-2.1.4.

An evaluation of undisturbed repository performance, which is defined to exclude human intrusion and unlikely disruptive natural events, is required by regulation (see 40 CFR § 191.15 and 40 CFR § 191.24).  Evaluations of past and present natural geologic processes in the region indicate that none has the potential to breach the repository within 10,000 years (see the CCA, Appendix SCR, Section SCR.1).  Disposal system behavior is dominated by the coupled processes of rock deformation surrounding the excavation, fluid flow, and waste degradation.  Each of these processes can be described independently, but the extent to which they occur is affected by the others.

Rock deformation immediately around the repository begins as soon as excavation creates a disturbance in the stress field.  Stress relief results in some degree of brittle fracturing and the formation of a DRZ, which surrounds excavations in all deep mines including the WIPP repository.  For the WIPP, the DRZ is characterized by an increase in permeability and porosity, and it may ultimately extend a few meters (m) from the excavated region.  Salt will also deform by creep processes resulting from deviatoric stress, causing the salt to move inward and fill voids.  Salt creep will continue until deviatoric stress is dissipated and the system is once again at stress equilibrium (see the CRA-2004, Chapter 6.0, Section 6.4.3.1).

The ability of salt to creep, thereby healing fractures and filling porosity, is one of its fundamental advantages as a medium for geologic disposal of radioactive waste, and one reason it was recommended by the National Academy of Sciences (see the CCA, Chapter 1.0, Section 1.3).  Salt creep provides the mechanism for crushed salt compaction in the shaft seal system, yielding properties approaching those of intact salt within 200 years (see the CCA, Appendix SEAL, Appendix D, Section D5.2).  Salt creep will also cause the DRZ surrounding the shaft to heal rapidly around the concrete components of the seal system.  In the absence of elevated gas pressure in the repository, salt creep would also substantially compact the waste and heal the DRZ around the disposal region.  Fluid pressures can become large enough through the combined effect of salt creep reducing pore volumes, and gas generation from waste degradation processes, to maintain significant porosity (greater than 20%) within the disposal room throughout the performance period (see also the CRA-2004, Chapter 6.0, Section 6.4.3).

Characterization of the Salado indicates that fluid flow from the far field does not occur on time scales of interest in the absence of an artificially imposed hydraulic gradient (see the CRA-2004, Chapter 2.0, Section 2.1.3.4 for a description of Salado investigations).  This lack of fluid flow is the second fundamental reason for choosing salt as a medium for geologic disposal of radioactive waste.  Lack of fluid flow is a result of the extremely low permeability of evaporite rocks that make up the Salado.  Excavating the repository has disturbed the natural hydraulic gradient and rock properties, resulting in fluid flow.  Small quantities of interstitial brine present in the Salado move toward regions of low hydraulic potential, and brine seeps are observed in the underground repository.  The slow flow of brine from halite into more permeable anhydrite MBs, and then through the DRZ into the repository, is expected to continue as long as the hydraulic potential within the repository is below that of the far field.  The repository environment will also include gas, so the fluid flow must be modeled as a two-phase process. Initially, the gaseous phase will consist primarily of air trapped at the time of closure, although other gases may form from waste degradation.  In the PA, the gaseous phase pressure will rise due to creep closure, gas generation, and brine inflow, creating the potential for flow from the excavated region (see also the CRA-2004, Chapter 6.0, Section 6.4.3.2).

An understanding of waste degradation processes indicates that the gaseous phase in fluid flow and the repository’s pressure history will be far more important than if the initial air were the only gas present.  Waste degradation can generate significant additional gas by two processes (see also the CRA-2004, Chapter 6.0, Section 6.4.3.3 for historical perspective, and Leigh et al. 2005, Section 2.3 and Section 2.4 for changes):

1.     The generation of hydrogen (H2) gas by anoxic corrosion of steels, other iron (Fe)-based alloys, and aluminum (Al) and Al-based alloys

2.     The generation of carbon dioxide (CO2), hydrogen sulfide (H2S), and methane (CH4) by anaerobic microbial consumption of waste containing CPR materials

Coupling these gas-generation reactions to fluid-flow and salt-creep processes is complex.  Gas generation will increase fluid pressure in the repository, thereby decreasing the hydraulic gradient and deviatoric stress between the far field and the excavated region and inhibiting the processes of brine inflow and salt creep.  Anoxic corrosion will also consume brine as it breaks down water to oxidize steels and other Fe-based alloys and release H2.  Thus, corrosion has the potential to be a self-limiting process, in that as it consumes all water in contact with steels and other Fe-based alloys, it will cease.  Microbial reactions also require water, either in brine or the gaseous phase.  It is assumed that microbial reactions will result in neither the consumption nor production of water.

The total volume of gas generated by corrosion and microbial consumption may be sufficient to result in repository pressures that approach lithostatic.  Sustained pressures above lithostatic are not physically reasonable within the disposal system, because the more brittle anhydrite layers are expected to fracture if sufficient gas is present.  The conceptual model implemented in the PA causes anhydrite MB permeability and porosity to increase rapidly as pore pressure approaches and exceeds lithostatic.  This conceptual model for pressure-dependent fracturing approximates the hydraulic effect of pressure-induced fracturing and allows gas and brine to move more freely within the MBs at higher pressures (see the CRA-2004, Chapter 6.0, Section 6.4.5.2).

Overall, the behavior of the undisturbed disposal system will result in extremely effective isolation of the radioactive waste.  Concrete, clay, and asphalt components of the shaft seal system will provide an immediate and effective barrier to fluid flow through the shafts, isolating the repository until salt creep has consolidated the compacted crushed salt components and permanently sealed the shafts.  Around the shafts, the DRZ in halite layers will heal rapidly because the presence of the solid material within the shafts will provide rigid resistance to creep.  The DRZ around the shaft, therefore, will not provide a continuous pathway for fluid flow (see the CRA-2004, Chapter 6.0, Section 6.4.4).  Similarly, the panel closure will rigidly resist creep, leading to a build-up of compressive stress which in turn will cause a rapid elimination of the DRZ locally.  In PA, it is conservatively assumed that the DRZ does not heal around either the disposal region or the operations and experimental regions, and pathways for fluid flow may exist indefinitely to the overlying and underlying anhydrite layers (e.g., MB 139 and Anhydrites A and B).  Some quantity of brine will be present in the repository under most conditions and may contain actinides mobilized as both dissolved and colloidal species. Gas generation by corrosion and microbial degradation is expected to occur, and will result in elevated pressures within the repository.  These pressures will not significantly exceed lithostatic because the more brittle anhydrite layers will fracture and provide a pathway for gas to leave the repository.  Fracturing due to high gas pressures may enhance gas and brine migration from the repository, but gas transport will not contribute to the release of actinides from the disposal system.  Brine flowing out of the waste disposal region through anhydrite layers may transport actinides as dissolved and colloidal species.  However, the quantity of actinides that may reach the accessible environment boundary through the interbeds during undisturbed repository performance is insignificant and has no effect on the compliance determination.  No migration of radionuclides is expected to occur vertically through the Salado (see Section PA-7.0, and Ismail and Garner 2008).

The WIPP PA is required by the performance standards to consider scenarios that include intrusions into the repository by inadvertent and intermittent drilling for resources.  The probability of these intrusions is based on a future drilling rate. This rate was calculated using the method outlined in Section 33, which analyzes the past record of drilling events in the Delaware Basin.  Active institutional controls (AICs) are assumed to prevent intrusion during the first 100 years after closure (40 CFR § 194.41). Passive institutional controls (PICs) were assumed in the CCA to effectively reduce the drilling rate by two orders of magnitude for the 600-year period following 100 years of active control.  However, in certifying the WIPP, the EPA denied credit for the effectiveness of passive controls for 600 years (U.S. Environmental Protection Agency 1998a). Although the CRA-2009 PA does not include a reduced drilling intrusion rate to account for PICs, future PAs may do so.  Future drilling practices are assumed to be the same as current practice, also consistent with regulatory criteria.  These practices include the type and rate of drilling, emplacement of casing in boreholes, and the procedures implemented when boreholes are plugged and abandoned (see 40 CFR § 194.33).

Human intrusion by drilling may cause releases from the disposal system through five mechanisms:

1.     Cuttings, which include material intersected by the rotary drilling bit

2.     Cavings, which include material eroded from the borehole wall during drilling

3.     Spallings, which include solid material carried into the borehole during rapid depressurization of the waste disposal region

4.     DBRs, which include contaminated brine that may flow to the surface during drilling

5.     Long-term brine releases, which include the contaminated brine that may flow through a borehole after it is abandoned

The first four mechanisms immediately follow an intrusion event and are collectively referred to as direct releases.  The accessible environment boundary for these releases is the ground surface.  The fifth mechanism, An transport by long-term groundwater flow, begins when concrete plugs are assumed to degrade in an abandoned borehole and may continue throughout the regulatory period.  The accessible environment boundary for these releases is the lateral subsurface limit of the controlled area (CRA-2004, Chapter 6.0, Section 6.0.2.3).

Repository conditions prior to intrusion will be the same as those for undisturbed repository performance, and all processes active in undisturbed repository performance will continue to occur following intrusion.  Because an intrusion provides a pathway for radionuclides to reach the ground surface and enter the geological units above the Salado, additional processes will occur that don’t in the undisturbed repository performance.  These processes include the mobilization of radionuclides as dissolved and colloidal species in repository brine and groundwater flow, and subsequent An transport in the overlying units.  Flow and transport in the Culebra are of particular interest because it is the most transmissive unit above the repository.  Thus, the Culebra is a potential pathway for lateral migration of contaminated brine in the event of a drilling intrusion accompanied by significant flow up the intrusion borehole (see the CRA-2004, Chapter 6.0, Section 6.4.6.2).

In a rotary drilling operation, the volume of material brought to the surface as cuttings is calculated as the cylinder defined by the thickness of the unit and the diameter of the drill bit.  The quantity of radionuclides released as cuttings is therefore a function of the activity of the intersected waste and the diameter of the intruding drill bit.  The DOE uses a constant value of 0.31115 m (12.25 inches [in]), consistent with bits currently used at the WIPP depth in the Delaware Basin (see the CRA-2004, Chapter 6.0, Section 6.4.12.5).  The intersected waste activity may vary depending on the type of waste intersected.  The DOE considers random penetrations into remote-handled (RH) transuranic (TRU) (RH-TRU) waste and each of the 690 different waste streams (Leigh, Trone, and Fox 2005, Section 4.4) identified for contact-handled (CH) transuranic (TRU) (CH-TRU) waste (569 and 693 waste streams were used in the CCA and the CRA-2004, respectively; see the CRA-2004, Chapter 6.0, Section 6.0.2.3.1).

The volume of particulate material eroded from the borehole wall by the drilling fluids and brought to the surface as cavings may be affected by the drill bit diameter, effective shear resistance of the intruded material, speed of the drill bit, viscosity of the drilling fluid and rate at which it is circulated in the borehole, and other properties related to the drilling process.  The most important of these parameters, after drill bit diameter, is the effective shear resistance of the intruded material (Leigh et al. 2005, Section 7.2).  In the absence of data describing the reasonable and realistic future properties of degraded waste and magnesium oxide (MgO) backfill, the DOE used conservative parameter values based on the properties of fine-grained sediment.  Other properties are assigned fixed values consistent with current practice.  The quantity of radionuclides released as cavings depends on the volume of eroded material and its activity, which is treated in the same manner as the activity of cuttings (see also Section PA-4.5 and Section PA-6.8.2.1).

Unlike releases from cuttings and cavings, which occur with every modeled borehole intrusion, spalling releases will occur only if pressure in the waste-disposal region exceeds the hydrostatic pressure in the borehole.  At lower pressures, below about 8 megapascals (MPa), fluid in the waste-disposal region will not flow toward the borehole.  At higher pressures, gas flow toward the borehole may be sufficiently rapid to cause additional solid material to enter the borehole.  If spalling occurs, the volume of spalled material will be affected by the physical properties of the waste, such as its tensile strength and particle diameter. The DOE based the parameter values used in the PA on reasonable and conservative assumptions.  Since the CCA, a revised conceptual model for the spallings phenomena has been developed (see the CRA-2004, Appendix PA, Section PA-4.6 and Appendix PA, Attachment MASS, Section MASS-16.1.3). Model development, execution, and sensitivity studies necessitated implementing parameter values pertaining to waste characteristics, drilling practices, and physics of the process.  The parameter range for particle size was derived by expert elicitation (U.S. Department of Energy 1997).

The quantity of radionuclides released as spalled material depends on the volume of spalled waste and its activity.  Because spalling may occur at a greater distance from the borehole than cuttings and cavings, spalled waste is assumed to have the volume-averaged activity of CH-TRU waste, rather than the sampled activities of individual waste streams.  The low permeability of the region surrounding the RH-TRU waste means it is isolated from the spallings process and does not contribute to the volume or activity of spalled material (see also Section PA-4.6 and Section PA-6.8.2.2 for more description of the spallings model).

Radionuclides may be released to the accessible environment if repository brine enters the borehole during drilling and flows to the ground surface.  The quantity of radionuclides released by direct brine flow depends on the volume of brine reaching the ground surface and the concentration of radionuclides contained in the brine.  As with spallings, DBRs will not occur if repository pressure is below the hydrostatic pressure in the borehole.  At higher repository pressures, mobile brine present in the repository will flow toward the borehole.  If the volume of brine flowing from the repository into the borehole is small, it will not affect the drilling operation, and flow may continue until the driller reaches the base of the evaporite section and installs casing in the borehole (see also Section PA-4.7 and Section PA-6.8.2.3).

Actinides may be mobilized in repository brine in two principal ways:

1.     As dissolved species

2.     As colloidal species

The solubilities of actinides depend on their oxidation states, with the more reduced forms (for example, III and IV oxidation states) being less soluble than the oxidized forms (V and VI).  Conditions within the repository will be strongly reducing because of large quantities of metallic Fe in the steel containers and the waste, and—in the case of plutonium (Pu)—only the lower-solubility oxidation states (Pu(III) and Pu(IV)) will persist.  Microbial activity will also help create reducing conditions.  Solubilities also vary with pH.  The DOE is therefore emplacing MgO in the waste-disposal region to ensure conditions that reduce uncertainty and establish low An solubilities.  MgO consumes CO2 and buffers pH, lowering An solubilities in WIPP brines (see Appendix SOTERM-2009 and Appendix MgO-2009).  Solubilities in the PA are based on the chemistry of brines that might be present in the waste-disposal region, reactions of these brines with the MgO engineered barrier, and strongly reducing conditions produced by anoxic corrosion of steels and other Fe-based alloys.

The waste contains organic ligands that could increase An solubilities by forming complexes with dissolved An species.  However, these organic ligands also form complexes with other dissolved metals, such as magnesium (Mg), calcium (Ca), Fe, lead (Pb), vanadium (V), chromium (Cr), manganese (Mn), and nickel (Ni), that will be present in repository brines due to corrosion of steels and other Fe-based alloys.  The CRA-2009 PA speciation and solubility calculations include the effect of organic ligands but not the beneficial effect of competition with Fe, Pb, V, Cr, Mn, and Ni.  (Appendix SOTERM-2009, Section SOTERM-2.3.6 and Section SOTERM-4.6, and Brush and Xiong 2005).

Colloidal transport of actinides has been examined, and four types of colloids have been determined to represent the possible behavior at the WIPP.  These include microbial colloids, humic substances, An intrinsic colloids, and mineral fragments.  Concentrations of An mobilized as these colloidal forms are included in the estimates of total An concentrations used in PA (Appendix SOTERM-2009, Section SOTERM-3.8.1 and Section SOTERM-4.7, and Garner and Leigh 2005).

Long-term releases to the ground surface or groundwater in the Rustler Formation (hereafter referred to as the Rustler) or overlying units may occur after the borehole has been plugged and abandoned.  In keeping with regulatory criteria, borehole plugs are assumed to have properties consistent with current practice in the basin.  Thus, boreholes are assumed to have concrete plugs emplaced at various locations. Initially, concrete plugs effectively limit fluid flow in the borehole.  However, under most circumstances, these plugs cannot be expected to remain fully effective indefinitely.  For the purposes of PA, discontinuous borehole plugs above the repository are assumed to degrade 200 years after emplacement.  From then on, the borehole is assumed to fill with a silty-sand-like material containing degraded concrete, corrosion products from degraded casing, and material that sloughs into the hole from the walls.  Of six possible plugged borehole configurations in the Delaware Basin, three are considered either likely or adequately representative of other possible configurations; one configuration (a two-plug configuration) is explicitly modeled in the flow and transport model (see Section PA-3.7 and Appendix MASS-2009, Section MASS-16.3).

If sufficient brine is available in the repository, and if pressure in the repository is higher than in the overlying units, brine may flow up the borehole following plug degradation.  In principle, this brine could flow into any permeable unit or to the ground surface if repository pressure were high enough.  For modeling purposes, brine is allowed to flow only into the higher-permeability units and to the surface.  Lower-permeability anhydrite and mudstone layers in the Rustler are treated as if they were impermeable to simplify the analysis while maximizing the amount of flow into units where it could potentially contribute to disposal system releases.  Model results indicate that essentially all flow occurs into the Culebra, which has been recognized since the early stages of site characterization as the most transmissive unit above the repository and the most likely pathway for subsurface transport (see also the CRA-2004, Chapter 2.0, Section 2.2.1.4.1.2).

Site characterization activities in the units above the Salado have focused on the Culebra.  These activities have shown that the direction of groundwater flow in the Culebra varies somewhat regionally, but in the area that overlies the repository, flow is southward.  These characterization and modeling activities conducted in the units above the Salado confirm that the Culebra is the most transmissive unit above the Salado.  The Culebra is the unit into which actinides are likely to be introduced from long-term flow up an abandoned borehole.  Regional variation in the Culebra’s groundwater flow direction is influenced by the transmissivity observed, as well as the lateral (facies) changes in the lithology of the Culebra in the groundwater basin where the WIPP is located.  Site characterization activities have provided no evidence of karst groundwater systems in the controlled area, although groundwater flow in the Culebra is affected by the presence of fractures, fracture fillings, and vuggy pore features (see Appendix HYDRO-2009 and the CRA-2004, Chapter 2.0, Section 2.1.3.5).  Other laboratory and field activities have focused on the behavior of dissolved and colloidal actinides in the Culebra.

Basin-scale regional modeling of three-dimensional groundwater flow in the units above the Salado demonstrates that it is appropriate, for the purposes of estimating radionuclide transport, to conceptualize the Culebra as a two-dimensional confined aquifer (see the CRA-2004, Chapter 2.0, Section 2.2.1.1).  Uncertainty in the flow field is incorporated by using 100 different geostatistically based T fields, each of which is consistent with available head and transmissivity data (Appendix PA-2009, Appendix TFIELD-2009).

Groundwater flow in the Culebra is modeled as a steady-state process, but two mechanisms considered in the PA could affect flow in the future.  Potash mining in the McNutt Potash Zone (hereafter referred to as the McNutt) of the Salado, which occurs now in the Delaware Basin outside the controlled area and may continue in the future, could affect flow in the Culebra if subsidence over mined areas causes fracturing or other changes in rock properties (see the CRA-2004, Chapter 6.0, Section 6.3.2.3).  Climatic changes during the next 10,000 years may also affect groundwater flow by altering recharge to the Culebra (see the CRA-2004, Chapter 6.0, Section 6.4.9 and the CCA, Appendix CLI).

Consistent with regulatory criteria of 40 CFR § 194.32, mining outside the controlled area is assumed to occur in the near future, and mining within the controlled area is assumed to occur with a probability of 1 in 100 per century (adjusted for the effectiveness of AICs during the first 100 years after closure).  Consistent with regulatory guidance, the effects of mine subsidence are incorporated in PA by increasing the transmissivity of the Culebra over the areas identified as mineable by a factor sampled from a uniform distribution between 1 and 1000 (U.S. Environmental Protection Agency 1996a, p. 5229).  T fields used in PA are therefore adjusted and steady-state flow fields calculated accordingly; once for mining that occurs only outside the controlled area, and once for mining that occurs both inside and outside the controlled area (Appendix TFIELD-2009, Section 9.0).  Mining outside the controlled area is considered in both undisturbed and disturbed repository performance.

The extent to which the climate will change during the next 10,000 years and how such a change will affect groundwater flow in the Culebra are uncertain.  Regional three-dimensional modeling of groundwater flow in the units above the Salado indicates that flow velocities in the Culebra may increase by a factor of 1 to 2.25 for reasonably possible future climates (see the CCA, Appendix CLI).  This uncertainty is incorporated in PA by scaling the calculated steady-state specific discharge within the Culebra by a sampled parameter within this range.

Field tests have shown that the Culebra is best characterized as a double-porosity medium for estimating contaminant transport in groundwater (see the CRA-2004, Chapter 2.0, Section 2.2.1.4.1.2 and Appendix HYDRO-2009).  Groundwater flow and advective transport of dissolved or colloidal species and particles occurs primarily in a small fraction of the rock’s total porosity and corresponds to the porosity of open and interconnected fractures and vugs.  Diffusion and slower advective flow occur in the remainder of the porosity, which is associated with the low-permeability dolomite matrix.  Transported species, including actinides (if present), will diffuse into this porosity.

Diffusion from the advective porosity into the dolomite matrix will retard An transport through two mechanisms.  Physical retardation occurs simply because actinides that diffuse into the matrix are no longer transported with the flowing groundwater.  Transport is interrupted until they diffuse back into the advective porosity.  In situ tracer tests have demonstrated this phenomenon.  Chemical retardation also occurs within the matrix as actinides are sorbed onto dolomite grains.  The relationship between sorbed and liquid concentrations is assumed to be linear and reversible.  The distribution coefficients (Kds) that characterize the extent to which actinides will sorb on dolomite were based on experimental data (see the CRA-2004, Chapter 6.0, Section 6.4.6.2).

Human intrusion scenarios evaluated in the PA include both single intrusion events and combinations of multiple boreholes.  Two different types of boreholes are considered:  those that penetrate a pressurized brine reservoir in the underlying Castile Formation (hereafter referred to as the Castile), and those that do not.

The presence of a brine reservoir under the repository is speculative, but on the basis of current information cannot be ruled out.  A pressurized brine reservoir was encountered at the WIPP-12 borehole within the controlled area to the north of the disposal region, and other pressurized brine reservoirs associated with regions of deformation in the Castile have been encountered elsewhere in the Delaware Basin (see the CRA-2004, Chapter 2.0, Section 2.2.1.2.2).  Based on a geostatistical analysis of the geophysical data of brine encounters in the region, the DOE estimates that there is a 0.08 probability that a random borehole penetrating waste in the WIPP will also penetrate an underlying brine reservoir (see the CCA, Appendix MASS, Attachment 18-6).  Upon their review of the CCA, the EPA determined that the DOE should treat this probability as uncertain, ranging from 0.01 to 0.60 in the CCA PAVT.  The EPA also required the DOE to modify the assumptions concerning Castile properties to increase the brine reservoir volumes (U.S. Environmental Protection Agency 1998b; Technical Support Document for 194.23:  Parameter Justification Report, Section 5).  The EPA determined that changing the rock compressibility and porosity of the Castile effectively modified the sampled brine reservoir volume to include the possibility of larger brine reservoir volumes like those encountered by the WIPP-12 borehole.

The primary consequence of penetrating a pressurized reservoir is to provide an additional source of brine beyond that which might flow into the repository from the Salado.  Direct releases at the ground surface resulting from the first repository intrusion would be unaffected by additional Castile brine, even if it flowed to the surface, because brine moving straight up a borehole will not significantly mix with waste.  However, the presence of Castile brine could significantly increase radionuclide releases in two ways.  First, the volume of contaminated brine that could flow to the surface may be greater for a second or subsequent intrusion into a repository that has already been connected by a previous borehole to a Castile reservoir.  Second, the volume of contaminated brine that may flow up an abandoned borehole after plug degradation may be greater for combinations of two or more boreholes that intrude the same panel if one of the boreholes penetrates a pressurized reservoir.  Both processes are modeled in PA.

The DOE’s approach to demonstrating continued compliance is the PA, which is based on the criteria indicated in section 194.34.  The PA process comprehensively considers the FEPs relevant to disposal system performance (see Appendix SCR-2009).  Those FEPs shown by screening analyses to potentially affect performance are included in quantitative calculations using a system of loosely coupled computer models to describe the interaction of the repository with the natural system, both with and without human intrusion.  Uncertainty in parameter values is incorporated in the analysis by a Monte Carlo approach, in which multiple simulations (or realizations) are completed using sampled values for the imprecisely known input parameters (see the CRA-2004, Chapter 6.0, Section 6.1.5).  Distribution functions characterize the state of knowledge for these parameters, and each realization of the modeling system uses a different set of sampled input values.  A sample size of 300 results in 300 different values of each parameter.  Thus, there are 300 different sets (vectors) of input parameter values.  These 300 vectors were divided among 3 replicates.  Quality assurance activities demonstrate that the parameters, software, and analysis used in PA were the result of a rigorous process conducted under controlled conditions (40 CFR § 194.22).

Of the FEPs considered, exploratory drilling for natural resources was identified as the only disruption with sufficient likelihood and consequence of impacting releases from the repository. For each vector of parameters values, 10,000 possible futures (realizations) are simulated, where a single future is defined as a series of intrusion events that occur randomly in space and time (Section PA-2.2). Each of these futures is assumed to have an equal probability of occurring; hence a probability of 0.0001.  Cumulative radionuclide releases from the disposal system are calculated for each future, and CCDFs are constructed by sorting the releases from smallest to largest and then summing the probabilities across the future.  Mean CCDFs were then computed for the three replicates of sampled parameters (Section PA-2.2).

This section summarizes the results of the CRA-2009 PA and demonstrates that the WIPP continues to comply with the quantitative containment requirements in 40 CFR § 191.13(a).  The CRA-2009 PA is different than the original certification PA in the CCA and the PA in the CRA-2004 PABC because it includes additional information, changes, and new data required by 40 CFR § 194.15 recertification application requirements.  Table PA-1 details the changes and new information included in this PA.

The results of the CRA-2009 PA demonstrate that the repository continues to comply with the disposal standards.  The key metric for regulatory compliance is the mean CCDF.  Figure PA-1, which compares the overall mean CCDF for the CRA-2009 PA to the overall mean CCDF for the CRA-2004 PABC, demonstrates two key points.  First, the overall mean CCDF lies entirely below the limits specified in section 191.13(a).  Thus, the WIPP is in compliance with the containment requirements of Part 191.  Second, for any probability, the expected releases in the CRA-2009 PA are only slightly higher than those in the CRA-2004 PABC, primarily because of the increased drilling rate.

Figure PA-1.      Overall Mean CCDFs for Total Normalized Releases:  CRA-2009 PA and CRA-2004 PABC

Detailed results of the CRA-2009 PA are contained in Section PA-9.0, which describes sensitivity analyses conducted as the final step in the Monte Carlo analysis.  These sensitivity analyses indicate the relative importance of each sampled parameter in terms of its contribution to uncertainty in estimating disposal system performance.  Analyses also examine the sensitivity of intermediate performance measures to the sampled parameters.  Examples of such intermediate performance measures include the quantity of radionuclides released to the accessible environment by any one mechanism (for example, cuttings or DBRs), and other model results that describe conditions of interest, such as disposal region pressure.

Section PA-9.0 presents CCDF distributions for each replication of the analysis, mean CCDFs, and an overall mean CCDF with the 95% confidence interval estimated from the 3 independent mean distributions. All 300 individual CCDFs, as well as the overall mean CCDF determined from the 3 replicates, lie entirely below and to the left of the limits specified in section 191.13(a) (see Figure PA-79).  Thus, the WIPP continues to comply with the containment requirements of Part 191.  Comparing the results of the 3 replicates indicates that the sample size of 100 in each replicate is sufficient to generate a stable distribution of outcomes (see Figure PA-80).  Within the region of regulatory interest (that is, at probabilities greater than 10 -3/104 year [yr]), the mean CCDFs from each replicate are essentially indistinguishable from the overall mean.

As discussed in Section PA-9.1, Section PA-9.2, and Section PA-9.3, examining the normalized releases from cuttings and cavings, spallings, and DBRs provides insight into the relative importance of each release mode’s contribution to the mean CCDF’s location and the compliance determination.  Releases from cuttings and cavings dominate the mean CCDF at high probabilities, while DBRs dominate the mean CCDF at low probabilities.  Spallings are less important and have very little effect on the location of the mean.  Subsurface releases from groundwater transport are less than 10 -6 EPA units and make no contribution to the mean CCDF’s location.

Uncertainties characterized in the natural system and the interaction of waste with the disposal system environment show little variation between the location of the mean CCDFs of the three replicates, providing additional confidence in the compliance determination.  The natural and engineered barrier systems of the WIPP provide robust and effective containment of TRU waste even if the repository is penetrated by multiple borehole intrusions.

This section outlines the conceptual structure of the WIPP PA.  First, a discussion of the regulatory requirements is given.  The requirements of section 191.13 and section 194.34 (summarized in Section PA-2.2.1) lead to the identification of three main PA components:

1.      A probabilistic characterization of the likelihood for different futures to occur at the WIPP site over the next 10,000 years

2.      A procedure for estimating the radionuclide releases to the accessible environment associated with each possible future that could occur at the WIPP site over the next 10,000 years

3.      A probabilistic characterization of the uncertainty in the parameters used to estimate potential releases

The probabilistic methods employed in WIPP PA give rise to the CCDF specified in section 191.13(a) and the distributions specified by 40 CFR § 194.34(b).

The methodology employed in PA derives from the EPA’s standard for the geologic disposal of radioactive waste, Environmental Radiation Protection Standards for the Management and Disposal of Spent Nuclear Fuel, High-Level and Transuranic Radioactive Wastes (Part 191) (U.S. Environmental Protection Agency 1993), which is divided into three subparts.  40 CFR Part 191 Subpart A applies to a disposal facility prior to decommissioning and establishes standards for the annual radiation doses to members of the public from waste management and storage operations.  Part 191 Subpart B applies after decommissioning and sets probabilistic limits on cumulative releases of radionuclides to the accessible environment for 10,000 years (section 191.13) and assurance requirements to provide confidence that section 191.13 will be met (40 CFR § 191.14).  Part 191 Subpart B also sets limits on radiation doses to members of the public in the accessible environment for 10,000 years of undisturbed repository performance (section 191.15).  40 CFR Part 191 Subpart C limits radioactive contamination of groundwater for 10,000 years after disposal (section 191.24).  In this recertification application, the DOE must demonstrate a reasonable expectation that the WIPP will continue to comply with the requirements of Part 191 Subparts B and C.

The following is the central requirement in Part 191 Subpart B, and the primary determinant of the PA methodology (U.S. Environmental Protection Agency 1985, p. 38086).

§ 191.13 Containment Requirements:

(a) Disposal systems for spent nuclear fuel or high-level or transuranic radioactive wastes shall be designed to provide a reasonable expectation, based upon performance assessments, that cumulative releases of radionuclides to the accessible environment for 10,000 years after disposal from all significant processes and events that may affect the disposal system shall:

(1) Have a likelihood of less than one chance in 10 of exceeding the quantities calculated according to Table 1 (Appendix A); and

(2) Have a likelihood of less than one chance in 1,000 of exceeding ten times the quantities calculated according to Table 1 (Appendix A).

(b) Performance assessments need not provide complete assurance that the requirements of 191.13(a) will be met.  Because of the long time period involved and the nature of the events and processes of interest, there will inevitably be substantial uncertainties in projecting disposal system performance.  Proof of the future performance of a disposal system is not to be had in the ordinary sense of the word in situations that deal with much shorter time frames.  Instead, what is required is a reasonable expectation, on the basis of the record before the implementing agency, that compliance with 191.13(a) will be achieved.

Section 191.13(a) refers to “quantities calculated according to Table 1 (Appendix A),” which means a normalized radionuclide release to the accessible environment based on the type of waste being disposed, the initial waste inventory, and the size of release that may occur (U.S. Environmental Protection Agency 1985, Appendix A).  Table 1 of Appendix A specifies allowable releases (i.e., release limits) for individual radionuclides and is reproduced as Table PA-2.  The WIPP is a repository for TRU waste, which is defined as “waste containing more than 100 nanocuries of alpha-emitting TRU isotopes, with half-lives greater than twenty years, per gram of waste” (U.S. Environmental Protection Agency 1985, p. 38084).  The normalized release R for TRU waste is defined by

                                                                                                  (PA.1)

where Qi is the cumulative release of radionuclide i to the accessible environment during the 10,000-year period following closure of the repository (curies [Ci]), Li is the release limit for radionuclide i given in Table PA-2 (Ci), and C is the amount of TRU waste emplaced in the repository (Ci).  In the CRA-2009 PA, C = 2.32 ´ 106 Ci (Leigh and Trone 2005, Section 3).  Further, “accessible environment” means (1) the atmosphere, (2) land surfaces, (3) surface waters, (4) oceans, and (5) all of the lithosphere beyond the controlled area.  “Controlled area” means (1) a surface location, to be identified by PICs, that encompasses no more than 100 square kilometers (km2) and extends horizontally no more than 5 kilometers (km) in any direction from the outer boundary of the original radioactive waste’s location in a disposal system, and (2) the subsurface underlying such a location (section 191.12).

PAs are the basis for addressing the containment requirements.  To help clarify the intent of Part 191, the EPA promulgated 40 CFR Part 194 (2004), Criteria for the Certification and

Table PA-2.       Release Limits for the Containment Requirements (U.S. Environmental Protection Agency 1985, Appendix A, Table 1)

Radionuclide

Release Limit Li per 1000 MTHMa or Other Unit of Wasteb

Americium-241 or -243

100

Carbon-14

100

Cesium-135 or -137

1,000

Iodine-129

100

Neptunium-237

100

Pu-238, -239, -240, or -242

100

Radium-226

100

Strontium-90

1,000

Technetium-99

10,000

Thorium (Th) -230 or -232

10

Tin-126

1,000

Uranium (U) -233, -234, -235, -236, or -238

100

Any other alpha-emitting radionuclide with a half-life greater than 20 years

100

Any other radionuclide with a half-life greater than 20 years that does not emit alpha particles

1,000

a   Metric tons of heavy metal (MTHM) exposed to a burnup between 25,000 megawatt-days (MWd) per metric ton of heavy metal (MWd/MTHM) and 40,000 MWd/MTHM.

b   An amount of TRU wastes containing one million Ci of alpha-emitting TRU radionuclides with half-lives greater than 20 years.

Recertification of the Waste Isolation Pilot Plant’s Compliance with the Part 191 Disposal Regulations.  There, an elaboration on the intent of section 191.13 is prescribed.

§ 194.34 Results of performance assessments.

(a) The results of performance assessments shall be assembled into “complementary, cumulative distributions functions” (CCDFs) that represent the probability of exceeding various levels of cumulative release caused by all significant processes and events.

(b) Probability distributions for uncertain disposal system parameter values used in performance assessments shall be developed and documented in any compliance application.

(c) Computational techniques, which draw random samples from across the entire range of the probability distributions developed pursuant to paragraph (b) of this section, shall be used in generating CCDFs and shall be documented in any compliance application.

(d) The number of CCDFs generated shall be large enough such that, at cumulative releases of 1 and 10, the maximum CCDF generated exceeds the 99th percentile of the population of CCDFs with at least a 0.95 probability.

(e) Any compliance application shall display the full range of CCDFs generated.

(f) Any compliance application shall provide information which demonstrates that there is at least a 95% level of statistical confidence that the mean of the population of CCDFs meets the containment requirements of § 191.13 of this chapter.

The DOE’s methodology for PA uses information about the disposal system and waste to evaluate performance over the 10,000-year regulatory time period.  To accomplish this task, the FEPs with potential to affect the future of the WIPP are first defined (Section PA-2.3.1).  Next, scenarios that describe potential future conditions in the WIPP are formed from logical groupings of retained FEPs (Section PA-2.3.2).  The scenario development process results in a probabilistic characterization for the likelihood of different futures that could occur at the WIPP (Section PA-2.2.2).  Using the retained FEPs, models are developed to estimate the radionuclide releases from the repository (Section PA-2.2.3).  Finally, uncertainty in model parameters is characterized probabilistically (Section PA-2.2.4).

As discussed in Section PA-2.3.1, the CCA PA scenario development process for the WIPP identified exploratory drilling for natural resources as the only disruption with sufficient likelihood and consequence of impacting releases from the repository (see the CCA, Appendix SCR).  In addition, Part 194 specifies that the occurrence of mining within the LWB must be included in the PA.  This has not changed for the CRA-2009 PA.  As a result, the projection of releases over the 10,000 years following closure of the WIPP is driven by the nature and timing of intrusion events.

The collection of all possible futures xst forms the basis for the probability space (Sst, Sst, pst) characterizing aleatory uncertainty, where Sst = {xst : xst is a possible future of the WIPP}, Sst is a suitably restricted collection of sets of futures, called “scenarios” (Section PA-3.10), and pst is a probability measure for the elements of Sst .  A possible future, xst,i, is thus characterized by the collection of intrusion events that occur in that future:

                                                                                         (PA.2)

where

n      is the number of drilling intrusions

tj      is the time (year) of the jth intrusion

lj      designates the location of the jth intrusion

ej     designates the penetration of an excavated or nonexcavated area by the jth intrusion

bj     designates whether or not the jth  intrusion penetrates pressurized brine in the Castile Formation

pj     designates the plugging procedure used with the jth intrusion (i.e., continuous plug, two discrete plugs, three discrete plugs)

aj     designates the type of waste penetrated by the jth intrusion (i.e., no waste, CH-TRU waste, RH-TRU waste, and, for CH-TRU waste, the waste streams encountered)

tmin   is the time at which potash mining occurs within the LWB

The subscript st indicates that aleatory (i.e., stochastic) uncertainty is being considered; the letters st are retained for historical context. The subscript i indicates that the future xst is one of many sample elements from Sst.

The probabilistic characterization of n, tj, lj,and ej is based on the assumption that drilling intrusions will occur randomly in time and space at a constant average rate (i.e., follow a Poisson process); the probabilistic characterization of bj derives from assessed properties of brine pockets; the probabilistic characterization of aj derives from the volumes of waste emplaced in the WIPP in relation to the volume of the repository; and the probabilistic characterization of pj derives from current drilling practices in the sedimentary basin (i.e., the Delaware Basin) in which the WIPP is located.  A vector notation is used for aj because it is possible for a given drilling intrusion to miss the waste or to penetrate different waste types (CH-TRU and RH-TRU) as well as to encounter different waste streams in the CH-TRU waste. Further, the probabilistic characterization for tmin follows from the criteria in Part 194 that the occurrence of potash mining within the LWB should be assumed to occur randomly in time (i.e., follow a Poisson process with a rate constant of lm = 10 -4 yr -1), with all commercially viable potash reserves within the LWB extracted at time tmin . In practice, the probability measure pst is defined by specifying probability distributions for each component of xst, as discussed further in Section PA-3.0.

Based on the retained FEPs (Section PA-2.3.1), release mechanisms include direct transport of material to the surface at the time of a drilling intrusion (i.e., cuttings, spallings, and brine flow) and release subsequent to a drilling intrusion due to brine flow up a borehole with a degraded plug (i.e., groundwater transport). The quantities of releases are determined by the state of the repository through time, which is determined by the type, timing, and sequence of prior intrusion events.  For example, pressure in the repository is an important determinant of spallings, and the amount of pressure depends on whether the drilling events that have occurred penetrated brine pockets and how long prior to the current drilling event the repository was inundated.

Computational models for estimating releases were developed using the retained FEPs; these models are summarized in Figure PA-2.  These computational models implement the conceptual models representing the repository system as described in section 194.23 and the mathematical models for physical processes presented in Section PA-4.0.  Most of the computational models involve the numerical solution of partial differential equations (PDEs) used to represent processes such as material deformation, fluid flow, and radionuclide transport.

The collection of computation models can be represented abstractly as a function f (xst|vsu), which quantifies the release that could result from the occurrence of a specific future xst and a specific set of values for model parameters vsu.  Because the future of the WIPP is unknown, the values of f (xst|vsu) are uncertain.  Thus, the probability space (Sst, Sst, pst), together with the function f (xst|vsu), give rise to the CCDF specified in section 191.13(a), as illustrated in Figure PA-3.  The CCDF represents the probability that a release from the repository greater than R will be observed, where R is a point on the abscissa (x-axis) of the graph (Figure PA-3).

Formally, the CCDF depicted in Figure PA-3 results from an integration over the probability space (Sst, Sst, pst):

                                                 (PA.3)

where dR( f (xst|vsu)) = 1 if f (xst|vsu) > R, dR( f (xst|vsu)) = 0 if f (xst|vsu) £ R, and dst(xst|vsu) is the probability density function associated with the probability space (Sst, Sst, pst).  In practice, the integral in Equation (PA.3) is evaluated by a Monte Carlo technique, where a random sample xst,i, i = 1, nR, is generated from Sst consistent with the probability distribution pst.  Using this random sample, Equation (PA.3) is numerically evaluated as

                                                  (PA.4)

The models in Figure PA-2 are too complex to permit a closed-form evaluation of the integral in Equation (PA.4) that defines the CCDF specified in Part 191.  In WIPP PA, these probability distribution functions (PDFs) are constructed using Monte Carlo simulation to sample the entire possible set of release outcomes. As long as the sampling is conducted properly and a sufficient number of samples is collected, the PDF of the sample should successfully approximate the PDF of the sample “universe” of all possible releases.

Figure PA-2.  Computational Models Used in PA

Figure PA-3.  Construction of the CCDF Specified in 40 CFR Part 191 Subpart B

In PA, the number of samples nR used to construct a CCDF is 10,000. However, the models in Figure PA-2 are also too computationally intensive to permit their evaluation for each of these 10,000 futures.  Due to this constraint, the models in Figure PA-2 are evaluated for a relatively small number of specific scenarios, and the results of these evaluations are used to construct CCDFs.  The representative scenarios are labeled E0, E1, E2, and E1E2, and are defined in Section PA-3.10; the procedure for constructing a CCDF from these scenarios is described in Section PA-6.0.

If the parameters used in the process-level models of Figure PA-2 were precisely known and if the models could accurately predict the future behavior of the repository, the evaluation of repository performance alone would be sufficient to answer the first three questions related to repository performance. However, the models do not perfectly represent the dynamics of the system and their parameters are not precisely known. Therefore, it is necessary to estimate the confidence one has in the CCDFs being constructed. The confidence in the CCDFs is established using Monte Carlo methods to evaluate how the uncertainty in the model parameters impacts the CCDFs or releases. The probabilistic characterization of the uncertainty in the model parameters is the outcome of the data development effort for the WIPP (summarized in Section 8.0 in Fox 2008).

Formally, uncertainty in the parameters that underlie the WIPP PA can be characterized by a second probability space (Ssu, Ssu, psu), where the sample space Ssu is defined by

                                  Ssu = {vsuvsu is a sampled vector of parameter values}                (PA.5)

The subscript su indicates that epistemic (i.e., subjective) uncertainty is being considered; the letters su are retained for historical context.  An element vsu Î Ssu is a vector vsu = vsu,1, vsu,2, …, vsu,N) of length N, where each element vsu,k is an uncertain parameter used in the models to estimate releases.  In practice, the probability measure psu is defined by specifying probability distributions for each element of vsu, discussed further in Section PA-5.0.

If the actual value for vsu were known, the CCDF resulting from evaluation of Equation (PA.4) could be determined with certainty and compared with the criteria specified in Part 191.  However, given the complexity of the WIPP site, the 10,000-year period under consideration, and the state of knowledge about the natural and engineered system, values for vsu are not known with certainty.  Rather, the uncertainty in vsu is characterized probabilistically, as described above, leading to a distribution of CCDFs (Figure PA-4) with each CCDF resulting from one of many vectors of values of vsu.  The uncertainty associated with the parameters is termed epistemic uncertainty, and has been referred to in WIPP PA documentation as subjective uncertainty.

WIPP PA uses a Monte Carlo procedure for evaluating the effects of epistemic uncertainty on releases. The procedure involves sampling the distributions assigned to the uncertain parameters and generating a CCDF of releases based on the results of the process-level models generated using those parameters values.  By repeating this process many times, a distribution of the

Figure PA-4.      Distribution of CCDFs Resulting from Possible Values for the Sampled Parameters

CCDFs can be constructed. The requirements of section 191.13 are evaluated, in part, using the mean probability of release.  The overall mean probability curve is created by averaging across the CCDFs for releases, i.e., averaging the CCDFs across vertical slices (Figure PA-4) (a formal definition is provided in Helton et al. 1998).  In addition, confidence limits on the mean are computed using standard t-statistics (Figure PA-5).  The proximity of these curves to the boundary line in Figure PA-3 indicates the confidence with which Part 191 will be met. Confidence is also established by examining the distribution of the CCDFs in relation to the release limits (Figure PA-6).

WIPP PA uses a stratified sampling design called LHS (McKay, Beckman, and Conover 1979) to generate a sample vsu, i = 1, …, nLHS, from Ssu consistent with the probability distribution psu.  LHS is an efficient scheme for sampling the range of a distribution using a relatively small sample. Based on order statistics, the sample size of nLHS = 300 replicates would provide coverage of 99% of the CCDF distribution with a confidence of 95%.

In Part 194, the EPA decided that the statistical portion of the determination of compliance with Part 191 will be based on the sample mean.  The LHS sample sizes should be demonstrated operationally to improve (reduce the size of) the confidence interval for the estimated mean.  The underlying principle is to show convergence of the mean (U.S. Environmental Protection Agency 1996b, p. 8-41).

Figure PA-5.      Example Mean Probability of Release and the Confidence Limits on the Mean from CRA-2009 PA

Figure PA-6.  Example CCDF Distribution From CRA-2009 PA (Replicate 1)


The DOE has chosen to demonstrate repeatability of the mean and to address the associated criteria of Part 194 using an operational approach of multiple replication, as proposed by Iman (1982). The complete set of PA calculations was repeated three times with all aspects of the analysis identical except for the random seed used to initiate the LHS procedure.  Thus, PA results are available for 3 replicates, each based on an independent set of 100 LHS vectors drawn from identical distributions for imprecisely known parameters and propagated through an identical modeling system.  This technique of multiple replication allows the adequacy of the sample size chosen in the Monte Carlo analysis to be evaluated and provides a suitable measure of confidence in the mean CCDF estimation used to demonstrate compliance with section 191.13(a).

This section addresses scenarios formed from FEPs that were retained for PA calculations, and introduces the specification of scenarios for consequence analysis.

The EPA has provided criteria concerning the scope of PAs in 40 CFR § 194.32.  In particular, criteria relating to the identification of potential processes and events that may affect disposal system performance are provided in 40 CFR § 194.32(e), which states

Any compliance application(s) shall include information which:

(1)  Identifies all potential processes, events or sequences and combinations of processes and events that may occur during the regulatory time frame and may affect the disposal system;

(2)  Identifies the processes, events or sequences and combinations of processes and events included in performance assessments; and

(3)  Documents why any processes, events or sequences and combinations of processes and events identified pursuant to paragraph (e)(1) of this section were not included in performance assessment results provided in any compliance application.

Section 32 of this application fulfills these criteria by documenting the DOE’s identification, screening, and screening results of all potential processes and events consistent with the criteria specified in section 194.32(e).

The first two steps in scenario development involve identifying and screening FEPs that are potentially relevant to the performance of the disposal system.  The CRA-2004, Chapter 6.0, Section 6.2 discusses the development of a comprehensive initial FEPs set used in the CCA, the methodology and criteria used for screening, the method used to reassess the CCA FEPs for the CRA-2004, and a summary of the FEPs retained for scenario development.  Changes to FEPs since the CRA-2004 are outlined in Section 32 and Appendix SCR-2009 of this application.

The original FEPs generation and screening were documented in the CCA, and the resulting FEPs list became the FEPs compliance baseline.  The baseline contained 237 FEPs and was documented the CCA, Appendix SCR.  The EPA compliance review of FEPs was documented in EPA’s Technical Support Document 194.32:  Scope of PA (U.S. Environmental Protection Agency 1998c).  The EPA numbered each FEP with a different scheme than the DOE used for the CCA.  The DOE has since adopted the EPA’s numbering scheme.

Logic diagrams illustrate the formation of scenarios for consequence analysis from combinations of events that remain after FEP screening (Cranwell et al. 1990) (Figure PA-7).  Each scenario shown in Figure PA-7 is defined by a combination of occurrence and nonoccurrence for all potentially disruptive events.  Disruptive events are defined as those that create new pathways or significantly alter existing pathways for fluid flow and, potentially, radionuclide transport within the disposal system.  Each of these scenarios also contains a set of features and nondisruptive events and processes that remain after FEP screening.  As shown in Figure PA-7, undisturbed repository performance (UP) and disturbed repository performance (DP) scenarios are considered in consequence modeling for the WIPP PA.  The UP scenario is used for compliance assessments (40 CFR § 194.54 and 40 CFR § 194.55).  Important aspects of UP and DP scenarios are summarized in this section.

The UP scenario is defined in 40 CFR § 191.12 to mean “the predicted behavior of a disposal system, including consideration of the uncertainties in predicted behavior, if the disposal system is not disrupted by human intrusion or the occurrence of unlikely natural events.”  For compliance assessments with respect to the Individual and Groundwater Protection Requirements (section 191.15; Appendix IGP-2009), it is only necessary to consider the UP scenario.  The UP scenario is also considered with DP scenario for PA with respect to the containment requirements (section 191.13).

No potentially disruptive natural events and processes are likely to occur during the regulatory time frame.  Therefore, all naturally occurring events and processes retained for scenario construction are nondisruptive and are considered part of the UP scenario.  Mining outside the LWB is assumed at the end of AIC for all scenarios.  The M scenario involves future mining within the controlled area. The disturbed repository E scenario involves at least one deep drilling event that intersects the waste disposal region. The M scenario and the E scenario may both occur in the future. The DOE calls a future in which both of these events occur the ME scenario.  More detailed descriptions are found in Section PA-2.3.2.2

The only natural features and waste- and repository-induced FEPs retained after screening that are excluded in the UP scenario, but included in the DP scenario, are those directly associated with the potential effects of future deep drilling within the controlled area.  Among the most significant FEPs that will affect the UP scenario within the disposal system are excavation-induced fracturing, gas generation, salt creep, and MgO in the disposal rooms.

·       The repository excavation and consequent changes in the rock stress field surrounding the excavated opening will create a DRZ immediately adjacent to excavated openings.  The DRZ will exhibit mechanical and hydrological properties different than those of the intact rock.

Figure PA-7.  Logic Diagram for Scenario Analysis

·       Organic material in the waste may degrade because of microbial activity, and brine will corrode metals in the waste and waste containers, with concomitant generation of gases.  Gas generation may result in pressures sufficient to both maintain or develop fractures and change the fluid flow pattern around the waste disposal region.

·       At the repository depth, salt creep will tend to heal fractures and reduce the permeability of the DRZ and the crushed salt component of the long-term shaft seals to near that of the host rock salt.

·       The MgO engineered barrier emplaced in the disposal rooms will react with CO2 and maintain mildly alkaline conditions.  Metal corrosion in the waste and waste containers will maintain reducing conditions.  These effects will maintain low radionuclide solubility.

Radionuclides can become mobile as a result of waste dissolution and colloid generation following brine flow into the disposal rooms.  Colloids may be generated from the waste (humics, mineral fragments, and An intrinsic colloids) or from other sources (humics, mineral fragments, and microbes).

Conceptually, there are several pathways for radionuclide transport within the undisturbed disposal system that may result in releases to the accessible environment (Figure PA-8).  Contaminated brine may migrate away from the waste-disposal panels if pressure within the panels is elevated by gas generated from corrosion or microbial consumption.  Radionuclide transport may occur laterally, through the anhydrite interbeds toward the subsurface boundary of the accessible environment in the Salado, or through access drifts or anhydrite interbeds to the base of the shafts.  In the latter case, if the pressure gradient between the panels and overlying strata is sufficient, contaminated brine may migrate up the shafts.  As a result, radionuclides may be transported directly to the ground surface, or laterally away from the shafts through permeable strata such as the Culebra, toward the subsurface boundary of the accessible environment.  These conceptual pathways are shown in Figure PA-8.

Figure PA-8.   Conceptual Release Pathways for the UP Scenario

The modeling system described in Section PA-4.0 includes potential radionuclide transport along other pathways, such as migration through Salado halite.  However, the natural properties of the undisturbed system make radionuclide transport to the accessible environment via these other pathways unlikely.

Assessments for compliance with section 191.13 need to consider the potential effects of future disruptive natural and human-initiated events and processes on the performance of the disposal system.  No potentially disruptive natural events and processes are considered sufficiently likely to require inclusion in analyses of either the UP or DP scenario.  The only future human-initiated events and processes retained after FEP screening are those associated with mining and deep drilling (but not the subsequent use of a borehole) within the controlled area or LWB when institutional controls cannot be assumed to eliminate the possibility of such activities (Section PA-3.2 and the CRA-2004, Chapter 6.0, Section 6.4.12.1).  In total, 21 disturbed repository FEPs associated with future mining and deep drilling have been identified.  These FEPs were assigned a screening designator of the DP scenario.

For evaluating the consequences of disturbed repository performance, the DOE has defined the mining scenario (M), the deep drilling scenario (E), and a mining and drilling scenario (ME).  These scenarios are described in the following sections.

The M scenario involves future mining within the controlled area.  Consistent with the criteria stated by the EPA in 40 CFR § 194.32(b) for PA calculations, the effects of potential future mining within the controlled area are limited to changes in hydraulic conductivity of the Culebra that result from subsidence (as described in Section PA-3.9).

Radionuclide transport may be affected in the M scenario if a head gradient between the waste-disposal panels and the Culebra causes brine contaminated with radionuclides to move from the waste-disposal panels to the base of the shafts and up to the Culebra.  The changes in the Culebra T field may affect the rate and direction of radionuclide transport within the Culebra.  Features of the M scenario are illustrated in Figure PA-9.

The three disturbed repository FEPs labeled M in the CRA-2004, Chapter 6.0, Table 6-9 are related to the occurrence and effects of future mining.  The modeling system used for the M scenario is similar to that developed for the UP scenario, but with a modified Culebra T field in the controlled area to account for the mining effects.

The disturbed repository E scenario involves at least one deep drilling event that intersects the waste disposal region.  The EPA provides criteria for analyzing the consequences of future drilling events in PA in 40 CFR § 194.33(c).

Performance assessments shall document that in analyzing the consequences of drilling events, the Department assumed that:

(1) Future drilling practices and technology will remain consistent with practices in the Delaware Basin at the time a compliance application is prepared.  Such future drilling practices shall include, but shall not be limited to: the types and amounts of drilling fluids; borehole depths, diameters, and seals; and the fraction of such boreholes that are sealed by humans; and

Figure PA-9.  Conceptual Release Pathways for the Disturbed Repository M Scenario

(2) Natural processes will degrade or otherwise affect the capability of boreholes to transmit fluids over the regulatory time frame.

Consistent with these criteria, there are several pathways for radionuclides to reach the accessible environment in the E scenario.  Before any deep drilling intersects the waste, potential release pathways are identical to those in the undisturbed repository scenario.

If a borehole intersects the waste in the disposal rooms, releases to the accessible environment may occur as material entrained in the circulating drilling fluid is brought to the surface. Particulate waste brought to the surface may include cuttings, cavings, and spallings.  Cuttings are the materials cut by the drill bit as it passes through waste.  Cavings are the materials eroded by the drilling fluid in the annulus around the drill bit.  Spallings are the materials forced into the circulating drilling fluid if there is sufficient pressure in the waste disposal panels.  During drilling, contaminated brine may flow up the borehole and reach the surface, depending on fluid pressure within the waste disposal panels.

When abandoned, the borehole is assumed to be plugged in a manner consistent with current practices in the Delaware Basin as prescribed in 40 CFR § 194.33(c)(1).  An abandoned intrusion borehole with degraded casing and/or plugs may provide a pathway for fluid flow and contaminant transport from the intersected waste panel to the ground surface if the fluid pressure within the panel is sufficiently greater than hydrostatic.  Additionally, if brine flows through the borehole to overlying units, such as the Culebra, it may carry dissolved and colloidal actinides that can be transported laterally to the accessible environment by natural groundwater flow in the overlying units.

Alternatively, the units intersected by an intrusion borehole may provide sources for brine flow to a waste panel during or after drilling. For example, in the northern Delaware Basin, the Castile, which underlies the Salado, contains isolated volumes of brine at fluid pressures greater than hydrostatic (as discussed in the CRA-2004, Chapter 2.0, Section 2.2.1.2.2).  The WIPP-12 penetration of one of these reservoirs provided data on one brine reservoir within the controlled area.  The location and properties of brine reservoirs cannot be reliably predicted; thus, the possibility of a deep borehole penetrating both a waste panel and a brine reservoir is accounted for in consequence analysis of the WIPP, as discussed in the CRA-2004, Chapter 6.0, Section 6.4.8.  Such a borehole could provide a connection for brine flow from the Castile to the waste panel, thus increasing fluid pressure and brine volume in the waste panel.

A borehole that is drilled through a disposal room pillar, but does not intersect waste, could also penetrate the brine reservoir underlying the waste disposal region.  Such an event would, to some extent, depressurize the brine reservoir, and thus would affect the consequences of any subsequent reservoir intersections.  The PA does not take credit for possible brine reservoir depressurization.

The DOE has distinguished two types of deep drilling events by whether or not the borehole intersects a Castile brine reservoir.  A borehole that intersects a waste disposal panel and penetrates a Castile brine reservoir is designated an E1 event.  A borehole that intersects a waste panel but does not penetrate a Castile brine reservoir is designated an E2 event. The consequences of deep drilling intrusions depend not only on the type of a drilling event, but on whether the repository was penetrated by an earlier E2 event or flooded due to an earlier E1 event.  The PA also does not take credit for depressurization of brine reservoirs from multiple drilling intrusions.  These scenarios are described in order of increasing complexity in the following sections.

The E2 scenario is the simplest scenario for inadvertent human intrusion into a waste disposal panel.  In this scenario, a panel is penetrated by a drill bit; cuttings, cavings, spallings, and brine flow releases may occur; and brine flow may occur in the borehole after it is plugged and abandoned.  Sources for brine that may contribute to long-term flow up the abandoned borehole are the Salado or, under certain conditions, the units above the Salado.  An E2 scenario may involve more than one E2 drilling event, although the flow and transport model configuration developed for the E2 scenario evaluates the consequences of futures that have only one E2 event.  Features of the E2 scenario are illustrated in Figure PA-10.

Figure PA-10.    Conceptual Release Pathways for the Disturbed Repository Deep Drilling E2 Scenario

Any scenario with exactly one inadvertent penetration of a waste panel that also penetrates a Castile brine reservoir is called E1.  Features of this scenario are illustrated in Figure PA-11.

Sources of brine in the E1 scenario are the brine reservoir, the Salado, and, under certain conditions, the units above the Salado.  However, the brine reservoir is conceptually the dominant source of brine in this scenario.  The flow and transport model configuration developed for the E1 scenario evaluates the consequences of futures that have only one E1 event.

Figure PA-11.    Conceptual Release Pathways for the Disturbed Repository Deep Drilling E1 Scenario

The E1E2 scenario is defined as all futures with multiple penetrations of a waste panel of which at least one intrusion is an E1.  One example of this scenario, with a single E1 event and a single E2 event penetrating the same panel, is illustrated in Figure PA-12.  However, the E1E2 scenario can include many possible combinations of intrusion times, locations, and types of event (E1 or E2).  The sources of brine in this scenario are those listed for the E1 scenario, and multiple E1 sources may be present.  The E1E2 scenario has a potential flow path not present in the E1 or E2 scenarios: flow from an E1 borehole through the waste to another borehole.  This flow path has the potential to (1) bring large quantities of brine in direct contact with waste and (2) provide a

Figure PA-12.    Conceptual Release Pathways for the Disturbed Repository Deep Drilling E1E2 Scenario

less restrictive path for this brine to flow to the units above the Salado (via multiple boreholes) compared to either the individual E1 or E2 scenarios.  It is both the presence of brine reservoirs and the potential for flow through the waste to other boreholes that make this scenario different from combinations of E2 boreholes in terms of potential consequences.  Estimates from flow and transport modeling are used to determine the extent of flow between boreholes and whether modeled combinations of E1 and E2 boreholes at specific locations in the repository should be treated as E1E2 scenarios or as independent E1 and E2 scenarios in the consequence analysis.

The M scenario and the E scenario may both occur in the future.  The DOE calls a future in which both of these events occur the ME scenario.  The occurrence of both mining and deep drilling do not create processes beyond those already described separately for the M and E scenarios.  For example, the occurrence of mining does not influence any of the interactions between deep boreholes and the repository or brine reservoirs, nor does the occurrence of drilling impact the effects of mining on Culebra hydrogeology.

The scenarios described in Section PA-2.3.2.1, Section PA-2.3.2.2, and Section PA-2.3.2.3 have been retained for consequence analysis to determine compliance with the containment requirements in section 191.13.  The modeling systems used to evaluate the consequences of these undisturbed and disturbed scenarios are discussed in Section PA-2.3.3.

Calculating scenario consequences requires quantitative modeling.  This section discusses the conceptual and computational models and some parameter values used to estimate the consequence of the scenarios described in Section PA-2.3.2.  Additional discussion of conceptual models and modeling assumptions is provided in Section PA-4.0.  Additional descriptions of sampled parameter values are included in Fox (2008).

A single modeling system was used to represent the disposal system and calculate the CCDFs.  The modeling system, however, can be conveniently described in terms of various submodels, with each describing a part of the overall system.  The models used in the WIPP PA, as in other complex analyses, exist at four different levels.

1.     Conceptual models are a set of qualitative assumptions that describe a system or subsystem for a given purpose.  At a minimum, these assumptions concern the geometry and dimensionality of the system, initial and boundary conditions, time dependence, and the nature of the relevant physical and chemical processes.  The assumptions should be consistent with one another and with existing information within the context of the given purpose.

2.     Mathematical models represent the processes at the site.  The conceptual models provide the context within which these mathematical models must operate, and define the processes they must characterize.  The mathematical models are predictive in the sense that, once provided with the known or assumed properties of the system and possible perturbations to the system, they predict the response of the system.  The processes represented by these mathematical models include fluid flow, mechanical deformation, radionuclide transport in groundwater, and removal of waste through intruding boreholes.

3.     Numerical models are developed to approximate mathematical model solutions because most mathematical models do not have closed-form solutions.

4.     The complexity of the system requires computer codes to solve the numerical models.  The implementation of the numerical model in the computer code with specific initial and boundary conditions and parameter values is generally referred to as the computational model.

Data are descriptors of the physical system being considered, normally obtained by experiment or observation. Parameters are values necessary in mathematical, numerical, or computational models.  The distinction between data and parameters can be subtle. Parameters are distinct from data, however, for three reasons:  (1) Data may be evaluated, statistically or otherwise, to generate model parameters to account for uncertainty in data.  (2) Some parameters have no relation to the physical system, such as the parameters in a numerical model to determine when an iterative solution scheme has converged.  (3) Many model parameters are applied at a different scale than one directly observed or measured in the physical system.  The distinction between data and parameter values is described further in Fox (2008) and Tierney (1990), where distribution derivations for specific parameters are given.  The interpretation and scaling of experimental and field data are discussed in Fox (2008) for individual and sampled parameters, as appropriate.


The PA for the WIPP identifies uncertainty in parameters and uncertainty in future events as distinctly different entities and requires sampling to be conducted in two dimensions. One dimension focuses on characterizing the uncertainty in terms of the probability that various possible futures will occur at the WIPP site over the next 10,000 years.  The other dimension characterizes the uncertainty due to lack of knowledge about the precise values of model parameters appropriate for the WIPP repository.  Each dimension of the analysis is characterized by a probability space.  Monte Carlo methods are used with the WIPP PA modeling system to sample each of the two probability spaces. 

Characterizing the probability distribution for the first dimension of the PA depends on identifying the kinds of events that could impact releases from the repository over the next 10,000 years. Screening analyses of possible future events concluded that the only significant events with the potential to affect radionuclide releases to the accessible environment are drilling and mining within the LWB (Appendix SCR-2009, Section SCR-5.0).  Consequently, modeling the future states of the repository focuses on representing the occurrences and effects of these two events. CCDFGF uses stochastic processes to simulate intrusion events by drilling and the occurrence of mining for natural resources. CCDFGF assembles the results from the deterministic models and selects the most appropriate scenario data provided by these models to use as the simulation of a 10,000-year future progresses.  Ten thousand potential futures are simulated and used to create distributions of potential releases, and then compiled into a single CCDF of potential releases.

WIPP PA is required not only to estimate the likelihood of future releases, but to establish confidence in those estimates.  Confidence is established using the second dimension of the analysis, which is based on the evaluation of uncertainty in the values of some of the parameters of the deterministic models.  This uncertainty is assumed to represent a lack of knowledge about the true values of the parameters, and is labeled epistemic uncertainty.  Epistemic uncertainty can be viewed as the representation of potential systematic errors in the results.  The impact of epistemic uncertainty on the results is determined by generating 300 sets of parameter values using a stratified random sampling design, LHS, and then running the deterministic models and CCDFGF with each set of sampled parameters.  Thus, 300 CCDFs are generated by CCDFGF. One set of parameters is often referred to as a vector.  The 300 simulations are organized as 3 replicates of 100 vectors each.  Because the uncertainty assigned to the parameters represents a lack of knowledge, this epistemic uncertainty could theoretically be reduced by collecting data to improve knowledge about the parameters.  Epistemic uncertainty is represented in the projections of potential releases from the repository by the variability among the 300 CCDFs.

The WIPP PA modeling system consists of a set of loosely coupled deterministic models (BRAGFLO, PANEL, NUTS, SECOTP2D, and CUTTINGS_S) that provide scenario-specific results to the code CCDFGF (Figure PA-2).  CCDFGF is, in contrast, a stochastic simulation model used to simulate potential futures of repository performance where drilling and mining intrusions can impact the state of the repository and produce release events. CCDFGF implements intrusions as stochastic events, thus incorporating the aleatory uncertainty associated with projections of future events.  This section describes how aleatory uncertainty is implemented in PA.  Epistemic uncertainty is discussed in Section PA-6.0.

As discussed in Section PA-2.2.2, aleatory uncertainty is defined by the possible futures xst,i  conditional on the set i of parameters used in Equation (PA.2).  Section PA-3.2, Section PA-3.3, Section PA-3.4, Section PA-3.5, Section PA-3.6, Section PA-3.7, Section PA-3.8, and Section PA-3.9 describe the individual components tj, ej, lj, bj, pj, aj, and tmin of xst,i and their associated probability distributions.  The concept of a scenario as a subset of the sample space of xst,i is discussed in Section PA-3.10.  The procedure used to sample the individual elements xst,i is described in Section PA-6.5.

AICs and PICs will be implemented at the WIPP site to deter human activity detrimental to repository performance.  AICs and PICs are described in detail in the CRA-2004, Chapter 7.0 and in appendices referenced in Chapter 7.0.  In this section, the impact of AICs and PICs on PA is described.

AICs will be implemented at the WIPP after final facility closure to control site access and ensure that activities detrimental to disposal system performance do not occur within the controlled area.  The AICs will preclude human intrusion in the disposal system.  A 100-year limit on the effectiveness of AICs in PA is established in 40 CFR § 191.14(a).  Because of the regulatory restrictions and the nature of the AICs that will be implemented, PA assumes there are no inadvertent human intrusions or mining in the controlled area for 100 years following repository closure.

PICs are designed to deter inadvertent human intrusion into the disposal system.  Only minimal assumptions were made about the nature of future society when designing the PICs to comply with the assurance requirements.  The preamble to Part 194 limits any credit for PICs in deterring human intrusion to 700 years after disposal (U.S. Environmental Protection Agency 1996a, p. 5231).  Although the DOE originally took credit for PICs in the CCA PA, but has not taken credit since.  Not including PICs is a conservative implementation, as no credit is taken for a beneficial process.

As described in Section PA-2.3.2.2, drilling intrusions in PA are assumed to occur randomly in time and space following a Poisson process.  Specifically, the drilling rate considered within the area marked by a berm as part of the system for PICs (Fox 2008, Table 46) is 5.85 ´ 10-3 intrusions per square kilometer per year (km2/yr).  AICs are assumed to prevent any drilling intrusions for the first 100 years after the decommissioning of the WIPP (Section PA-3.2).  In the computational implementation of PA, it is convenient to represent the Poisson process for drilling intrusions by its corresponding rate term ld( t) for intrusions into the area marked by the berm.  Specifically,

                                                                                      (PA.6)

where 0.6285 km2 is the area enclosed by the berm (Fox 2008, Table 45) and t is elapsed time since decommissioning the WIPP.

The function ld( t) defines the parameter of the exponential distribution that gives rise to the times of intrusions, tj of Equation (PA.2).  In the computational implementation of the analysis, the exponential distribution is sampled to define the times between successive drilling intrusions (Figure PA-13, Section PA-6.5).  A key assumption of the exponential distribution is that events are independent of each other, so the occurrence of one has no effect on the occurrence of the next event.  The process giving rise to such events is sometimes called a Poisson process because the distribution of such events over a fixed interval of time is a Poisson distribution. Due to the 10,000-year regulatory period specified in section 191.13, tj is assumed to be bounded above by 10,000 years in the definition of xst,i.  Further, tj is bounded below by 100 years as defined in Equation (PA.6).

Figure PA-13.  CDF for Time Between Drilling Intrusions

The variable ej is a designator for whether or not the jth drilling intrusion penetrates an excavated, waste-filled area of the repository:  ej = 0 or 1 implies penetration of nonexcavated or excavated area, respectively.  The corresponding probabilities P[ej = 0] and P[ej = 1] for ej = 0 and ej = 1 are

                                                     (PA.7)

                                                                                (PA.8)

where 0.1273 km2 and 0.6285 km2 are the excavated area of the repository and the area of the berm, respectively (Fox 2008, Table 45).

Locations of drilling intrusions through the excavated, waste-filled area of the repository are discretized to the 144 locations in Figure PA-14.  Assuming that a drilling intrusion occurs within the excavated area, it is assumed to be equally likely to occur at each of these 144 locations.  Thus, the probability pLk that drilling intrusion j will occur at location lk, k = 1, 2, ¼, 144 in Figure PA-14 is

                           (PA.9)

Figure PA-14.  Discretized Locations for Drilling Intrusions

The conceptual models for the Castile include the possibility that pressurized brine reservoirs underlie the repository (Section PA-4.2.10).  The variable bj is a designator for whether or not the jth drilling intrusion penetrates pressurized brine, where bj = 0 signifies nonpenetration and bj = 1 signifies penetration of pressurized brine.  In PA, the probability pB1 = P[b= 1] is sampled from a uniform distribution ranging from 0.01 to 0.60 (see GLOBAL:PBRINE in Table PA-19).

Three borehole plugging patterns, pk, are considered in PA:  (1) p1, a full concrete plug through the Salado to the Bell Canyon Formation (hereafter referred to as Bell Canyon), (2) p2, a two-plug configuration with concrete plugs at the Rustler/Salado interface and the Castile/Bell Canyon interface, and (3) p3, a three-plug configuration with concrete plugs at the Rustler/ Salado, Salado/Castile, and Castile/Bell Canyon interfaces.  The probability that a given drilling intrusion will be sealed with plugging pattern pk, k= 1, 2, 3, is given by pPLk, where pPL1 = P[k = 1] = 0.015, pPL2 = P[k = 2] = 0.696, pPL3 = P[k = 3] = 0.289 (Fox 2008, Table 46).

The waste intended for disposal at the WIPP is represented by 767 distinct waste streams, with 690 of these waste streams designated as CH-TRU waste and 77 designated as RH-TRU waste (Leigh, Trone, and Fox 2005, Section 4.4).  For the CRA-2009 PA, the 77 separate RH-TRU waste streams are represented by a single, combined RH-TRU waste stream.  The activity levels for the waste streams are given in Fox (2008, Table 48 and Table 49).  Each waste container emplaced in the repository contains waste from a single CH-TRU waste stream.  Waste packaged in 55-gallon (gal) drums is stacked 3 drums high within the repository.  Although waste in other packages (e.g., standard waste boxes, 10-drum overpacks, etc.) may not be stacked 3 high, PA assumes that each drilling intrusion into CH-TRU waste intersects 3 different waste streams.  In contrast, all RH-TRU waste is represented by a single waste stream, and so each drilling intrusion through RH-TRU waste is assumed to intersect this single waste stream.  Appendix MASS-2009, Section MASS-21.0 examines the sensitivity of PA results to the assumption that three waste streams are intersected by each drilling intrusion into CH-TRU waste.

The vector aj characterizes the type of waste penetrated by the jth drilling intrusion.  Specifically,

                                                        aj = 0 if ej = 0                                                (PA.10)

(i.e., if the ith drilling intrusion does not penetrate an excavated area of the repository);

                                      aj = 1 if ej = 1 and RH-TRU is penetrated                          (PA.11)

                         aj = [iCHj1, iCHj2, iCHj3] if ej = 1 and CH-TRU is penetrated           (PA.12)

where iCHj1, iCHj2, and iCHj3 are integer designators for the CH-TRU waste streams intersected by the jth drilling intrusion (i.e., each of iCHj1, iCHj2, and iCHj3 is an integer between 1 and 690).

Whether the jth intrusion penetrates a nonexcavated or excavated area is determined by the probabilities pE0 and pE1 discussed in Section PA-3.4.  The type of waste penetrated is determined by the probabilities pCH and pRH.  The excavated area used for disposal of CH-TRU waste (aCH) is 1.115 ´ 105 square meters (m2) and the area used for disposal of RH-TRU waste (aRH) is 1.576 ´ 104 m2 (Fox 2008, Table 45), for a total disposal area of aEX = aCH + aRH = 1.273 ´ 105 m2.  Given that the jth intrusion penetrates an excavated area, the probabilities pCH and pRH of penetrating CH-TRU and RH-TRU waste are given by


                                                                                                                                         (
PA.13)


                                                                                                                                         (
PA.14)

As indicated in this section, the probabilistic characterization of aj depends on a number of individual probabilities.  Specifically, pEx0 and pEx1 determine whether a nonexcavated or excavated area is penetrated (Section PA-3.5); pCH and pRH determine whether CH-TRU or RH-TRU waste is encountered, given penetration of an excavated area; and the individual waste stream probabilities in Fox (2008, Table 48 and Table 49), determine the specific waste streams iCHj1, iCHj2, and iCHj3 encountered given a penetration of CH-TRU waste. The probability of encountering a particular CH-TRU waste stream is computed as the ratio of the volume of that waste stream to the volume of CH-TRU waste.

Full mining of known potash reserves within the LWB is assumed to occur at time tmin.  The occurrence of mining within the LWB in the absence of institutional controls is specified as following a Poisson process with a rate of lm = 1 ´ 10 -4 yr -1 (Fox 2008, Table 46).  However, this rate can be reduced by AICs and PICs.  Specifically, AICs are assumed to result in no possibility of mining for the first 100 years after decommissioning of the WIPP.  In PA, PICs do not affect the mining rate. Thus, the mining rate lm( t) is

                                                                                   (PA.15)

                                                                (PA.16)

where t is the elapsed time since decommissioning of the WIPP.

In the computational implementation of the analysis, lm( t) is used to define the distribution of time to mining.  The use of lm( t) to characterize tmin is analogous to the use of ld to characterize the tj, except that only one mining event is assumed to occur (i.e., xst,i contains only one value for tmin) in order to be consistent with guidance given in Part 194 that mining within the LWB should be assumed to remove all economically viable potash reserves.  Due to the 10,000-year regulatory period specified in section 191.13, tmin is assumed to be bounded above by 10,000 years in the definition of xst,i.

A scenario is a subset of the sample space for aleatory uncertainty.  The underlying goal of scenario definition is to define the state of repository conditions prior to and following intrusion events.  Scenarios are specific cases of inputs or system states that are selected to cover the range of possible cases.  Given the complexity of the futures xst,i (see Equation (PA.2)), many different scenarios can be defined.  The computational complexity of the function f (xst|vsu) in Section PA-2.2.3 limits evaluation to only a few intrusion scenarios.  As presented in Section PA-2.3.2, PA considers four fundamental intrusion scenarios:

           E0  =  no drilling intrusion through an excavated area of the repository

           E1  =  a drilling intrusion through an excavated area of the repository that penetrates pressurized brine in the Castile

           E2  =  a drilling intrusion through an excavated area of the repository that does not penetrate pressurized brine in the Castile

      E1E2  =  two or more previous intrusions, at least one of which is an E1 intrusion

These definitions of intrusion scenarios capture the most important events impacting the state of the repository:  whether or not the repository is inundated by the penetration of a brine pocket, and whether or not there exists a possible route of release upward via a borehole. The state of the repository is also designated as E0, E1, E2, or E1E2. Scenarios for some of the process-level models consist of a single intrusion scenario occurring at specific times.  CCDFGF is used to simulate multiple intrusions over 10,000 years.

If only the intrusion scenarios controlled the state of the repository, then the state would be defined by the sequence of drilling events alone.  However, CCDFGF also considers the impact of plugging pattern on boreholes.  A borehole with a full plugging pattern that penetrates the waste area is also assumed to have no impact, and leaves the repository in its previous state, including the undisturbed state (see Section PA-6.8.4.1 and Figure PA-41 for more details).  Thus, an E2 intrusion event into an E0 repository will result in an E0 state if a full plugging pattern is used, or an E2 state otherwise.  An E1 intrusion subsequent to an E2 intrusion will leave the repository in an E1E2 state, where it will remain, regardless of subsequent intrusions.  It is therefore important to distinguish between the type of intrusion, listed above, and the state of the repository.

The probability that no excavated area will be penetrated during the 10,000-year interval can be computed using a distribution of the number of penetration events and the probability that a drilling event will penetrate the excavated area.  For the Poisson distribution of drilling events, the probability of there being n events in the 10,000-year history is

                                                                         (PA.17)

where ld is the mean drilling rate per year in the period following the period of AICs, 9,900 is the number of years in which drilling can occur after the institutional control period of 100 years, and n is the number of drilling events.  The probability of having n events all within the nonexcavated area is pEx0n, or specifically 0.797n.  Thus, the probability of having only events in the nonexcavated area over 10,000 years, i.e., having no drilling intrusions into the excavated area, is just the sum across all n of the products of the probability of having exactly n drilling events and the probability that all n events penetrate the unexcavated area:

                                                                                                (PA.18)

The calculated probability becomes

                                                                             (PA.19)

This probability is the lower bound on the probability of the repository being in an E0 state, given that it does not include the consideration of the plugging pattern.

The probability of a single E1, E2, or E1E2 intrusion over 10,000 years is relatively small. Assuming that pB1 takes on its mean value of 0.305 (see Section PA-3.6), and ignoring the impact of the plugging pattern, for a constant rate of drilling, ld, these equations are

                                                                    (PA.20)

and

                                                                   (PA.21)

respectively, where (pEx1 × ld)represents the annual rate of drilling into the excavated region of the repository which is multiplied by 9900 to give the rate per 9,900 years.  The probability of an intrusion into the excavated area is subsequently multiplied by the probability of hitting or missing a brine pocket.  In this form, it can be seen that the term for the probability for intrusion is equivalent to the PDF of the Poisson distribution for n = 1:

                                                                                                                  (PA.22)

The expressions defining the probability of being in the E0 state after 10,000 years and of having a single E1 or E2 intrusion event after 10,000 years are relatively simple because the scenarios E0, E1, and E2 are relatively simple.  The scenario E1E2 is more complex and, as a result, computing its probability is also more complex.  Closed-form formulas for the probabilities of quite complex scenarios can be derived, but they are very complicated and involve large numbers of iterated integrals (Helton 1993).  These probabilities of single E1 and E2 intrusions are relevant to the scenarios used by the process-level models.

CCDFGF simulates histories that can have many intrusion events (WIPP Performance Assessment 2003a).  The process-level models evaluate the releases at a small number of specific times for each of the four intrusion scenarios.  Releases from the repository are calculated using results from these fundamental scenarios (Section PA-6.7 and Section PA-6.8).  Releases for an arbitrary future are estimated from the results of these fundamental scenarios (Section PA-6.8); these releases are used to construct CCDFs by Equation (PA.4).

Previous WIPP PAs have used the Monte Carlo approach to construct the CCDF indicated in Equation (PA.4). The Monte Carlo approach generates releases for 10,000 possible futures.  CCDFs are constructed by treating the 10,000 releases values as order statistics; each release is assigned a probability of 1 ´ 10-4, and the CCDF can be constructed by plotting the complement of the sum of the probabilities ordered by the release value.  The CRA-2009 PA uses the same approach as the CRA-2004 PA.


This section describes how releases to the accessible environment are estimated for a particular future in PA.

The function f(xst,i) estimates the radionuclide releases to the accessible environment associated with each of the possible futures (xst,i) that could occur at the WIPP site over the next 10,000 years.  In practice, f(xst,i) is quite complex and is constructed by the models implemented in computer programs used to simulate important processes and releases at the WIPP.  In the context of these models, f(xst,i) has the form

                    (PA.23)

where

                                                    xst,i   ~  particular future under consideration

                                                  ~  future involving no drilling intrusions but a mining event at the same time tmin as in xst

                                           fC(xst,i) ~  cuttings and cavings release to accessible environment for xst,i calculated with CUTTINGS_S

                                           fB(xst,i) ~ two-phase flow in and around the repository calculated for xst,i with BRAGFLO; in practice, fB(xst,i) is a vector containing a large amount of information, including pressure and brine saturation in various geologic members

      ~ spallings release to accessible environment for xst,i calculated with the spallings model contained in DRSPALL and CUTTINGS_S; this calculation requires repository conditions calculated by fB(xst,i) as input

     ~ DBR to accessible environment for xst,i also calculated with BRAGFLO; this calculation requires repository conditions calculated by fB(xst,i) as input

        ~ release through anhydrite MBs to accessible environment for xst,i calculated with NUTS; this calculation requires flows in and around the repository calculated by fB(xst,i) as input

        ~ release through Dewey Lake to accessible environment for xst,i calculated with NUTS; this calculation requires flows in and around the repository calculated by fB(xst,i) as input

           ~ release to land surface due to brine flow up a plugged borehole for xst,i calculated with NUTS; this calculation requires flows in and around the repository calculated by fB(xst,i) as input

                                 ~ flow field in the Culebra calculated for xst,0 with MODFLOW; xst,0 is used as an argument to fMF because drilling intrusions are assumed to cause no perturbations to the flow field in the Culebra

         ~ release to Culebra for xst,i calculated with NUTS or PANEL as appropriate; this calculation requires flows in and around the repository calculated by fB(xst,i) as input

 ~ groundwater transport release through Culebra to accessible environment calculated with SECOTP2D.  This calculation requires MODFLOW results (i.e., fMF(xst,i)) and NUTS or PANEL results (i.e., ) as input

The remainder of this section describes the mathematical structure of the mechanistic models that underlie the component functions of f(xst,i) in Equation (PA.23).

The Monte Carlo CCDF construction procedure, implemented in the code CCDFGF (WIPP Performance Assessment 2003a), uses a sample of size nS = 10,000 in PA.  The individual programs that estimate releases do not run fast enough to allow this many evaluations of f.  As a result, a two-step procedure is being used to evaluate f in calculating the summation in Equation (PA.23).  First, f and its component functions are evaluated with the procedures (i.e., models) described in this section for a group of preselected futures.  Second, values of f(xst) for the randomly selected futures xst,i used in the numerical evaluation of the summation in Equation (PA.23) are then constructed from results obtained in the first step.  These constructions are described in Section PA-6.7 and Section PA-6.8, and produce the evaluations of f(xst) that are actually used in Equation (PA.23).

For notational simplicity, the functions on the right-hand side of Equation (PA.23) will typically be written with only xst as an argument (e.g., fSP(xst) will be used instead of fSP[xst, fB(xst)]).  However, the underlying dependency on the other arguments will still be present.

The major topics considered in this chapter are two-phase flow in the vicinity of the repository as modeled by BRAGFLO (i.e., fB) (Section PA-4.2), radionuclide transport in the vicinity of the repository as modeled by NUTS (i.e., fMB, fDL, fS, fNP) (Section PA-4.3), radionuclide transport in the vicinity of the repository as modeled by PANEL (i.e., fNP) (Section PA-4.4), cuttings and cavings releases to the surface as modeled by CUTTINGS_S (i.e., fC) (Section PA-4.5), spallings releases to the surface as modeled by DRSPALL and CUTTINGS_S (i.e., fSP) (Section PA-4.6), DBRs to the surface as modeled by BRAGFLO (i.e., fDBR) (Section PA-4.7), brine flow in the Culebra as modeled by MODFLOW (i.e., fMF) (Section PA-4.8), and radionuclide transport in the Culebra as modeled by SECOTP2D (i.e., fST) (Section PA-4.9).

Quantifying the effects of gas and brine flow on radionuclide transport from the repository requires a two-phase (brine and gas) flow code.  The two-phase flow code BRAGFLO is used to simulate gas and brine flow in and around the repository (Nemer 2007b and 2007c).  Additionally, the BRAGFLO code incorporates the effects of disposal room consolidation and closure, gas generation, and rock fracturing in response to gas pressure.  This section describes the mathematical models on which BRAGFLO is based, the representation of the repository in the model, and the numerical techniques employed in the solution.

Two-phase flow in the vicinity of the repository is represented by the following system of two conservation equations, two constraint equations, and three equations of state:

Gas Conservation

                               Ñ×           (PA.24)

Brine Conservation

                                 Ñ×             (PA.25)

Saturation Constraint

                                                                                                                      (PA.26)

Capillary Pressure Constraint

                                                                                                        (PA.27)

Gas Density

       rg (determined by Redlich-Kwong-Soave (RKS) equation of state; see Equation (PA.51))
                                                                                                                                         (PA.28)

Brine Density

                                                                                              (PA.29)

Formation Porosity

                                                                                                  (PA.30)

where

          g  =  acceleration due to gravity (meters per second squared [m])

          h  =  vertical distance from a reference location (m)

         krl  =  relative permeability (dimensionless) to fluid l, l = b (brine), g (gas)

         Pc  =  capillary pressure in Pascals (Pa)

         Pl  = pressure of fluid l (Pa)

         qrl  =  rate of production (or consumption, if negative) of fluid l due to chemical reaction (kilograms per cubic meter per seconds [kg/m3/s])

         ql  =  rate of injection (or removal, if negative) of fluid l (kg/m3/s)

         Sl =  saturation of fluid l (dimensionless)

           t  =  time (s)

         a  =  geometry factor (m)

         rl =  density of fluid l (kg/m3)

         ml =  viscosity of fluid l (Pa s)

          f  =  porosity (dimensionless)

         f0 =  reference (i.e., initial) porosity (dimensionless)

       Pb0   =  reference (i.e., initial) brine pressure (Pa), constant in Equation (PA.29) and spatially variable in Equation (PA.30)

         r0  =  reference (i.e., initial) brine density (kg/m3)

         c f  =  pore compressibility (Pa-1)

         cb  =  brine compressibility (Pa-1)

         K  =  permeability of the material (m2), isotropic for PA (Howarth and Christian-Frear 1997)

For the brine transport Equation (PA.25), the intrinsic permeability of the material is used.  For the gas transport Equation (PA.24), the permeability K is modified to account for the Klinkenberg effect (Klinkenberg 1941).  Specifically,

                                                                                                          (PA.31)

where a and b are gas and formation-dependent constants.  Values of a = -0.3410 and = 0.2710 were determined from data obtained for MB 139 (Christian-Frear 1996), with these values used for all regions in Figure PA-15.

The conservation equations are valid in one (i.e., Ñ = [/x]), two (i.e., Ñ = [ / x,  / y]), and three (i.e., Ñ = [ / x,  / y,  / z]) dimensions.  In PA, the preceding system of equations is used to model two-phase fluid flow within the two-dimensional region shown in Figure PA-15.  The details of this system are now discussed.

The a term in Equation (PA.24) and Equation (PA.25) is a dimension-dependent geometry factor and is specified by

a  =  area normal to flow direction in one-dimensional flow (i.e., DyDz; units = m2)
    =  thickness normal to flow plane in two-diemensional flow (i.e., Dz; units = m)
    =  1 in three-dimensional flow (dimensionless)                                              (
PA.32)

PA uses a two-dimensional geometry to compute two-phase flow in the vicinity of the repository, and as a result, a is the thickness of the modeled region (i.e., Dz) normal to the flow plane (Figure PA-15).  Due to the use of the two-dimensional grid in Figure PA-15, a is spatially dependent, with the values used for a defined in the column labeled “Dz.”  Specifically, a increases with distance away from the repository edge in both directions to incorporate the increasing pore volume through which fluid flow occurs.  The method used in PA, called rectangular flaring, is illustrated in Figure PA-16 and ensures that the total volume surrounding the repository is conserved in the numerical grid.  The equations and method used to determine a  for the grid shown in Figure PA-15 are described in detail by Stein (2002).

The h term in Equation (PA.24) and Equation (PA.25) defines vertical distance from a reference point.  In PA, this reference point is taken to be the center of MB 139 at the location of the shaft (i.e., (xref, yref) = (23664.9 m, 378.685 m), which is the center of cell 1266 in Figure PA-17).  Specifically, h is defined by

                                                                     (PA.33)

where q is the inclination of the formation in which the point (x, y) is located.  In PA, the Salado is modeled as having an inclination of 1 degree from north to south, and all other formations are modeled as being horizontal.  Thus, q = 1 degree for points within the Salado, and q = 0 degrees otherwise.  Treating the Salado as an inclined formation and treating the Castile, Castile brine reservoir, Rustler, and overlying units as horizontal creates discontinuities in the grid at the lower and upper boundaries of the Salado.  However, this treatment does not create a computational problem, since the Salado is isolated from vertical flow; its upper boundary adjoins the impermeable Los Medaños Member (formerly referred to as the Unnamed Member) at the base of the Rustler, and its lower boundary adjoins the impermeable Castile.

In the solution of Equations (PA.24) through (PA.30), Sb and Sg are functions of location and time.  Thus, Pc, krb, and krg are functions of the form Pc( x, y, t), krb( x, y, t), and krg( x, y, t).  In the computational implementation of the solution of the preceding equations, flow of phase l out of a computational cell (Figure PA-17) cannot occur when Sl( x, y, t) £ Slr( x, y, t), where Slr denotes the residual saturation for phase l.  The values used for Slr, l = b, g are summarized in Table PA-3.

Values for f0 and c f (Equation (PA.30)) are also given in Table PA-3.  Initial porosity f0 for the DRZ is a function of the uncertain parameter for initial halite porosity f0H (HALPOR; see Table PA-19) and is given by Martell (1996a) and Bean et al. (1996), Section 4):

                                                                f0 = f0H + 0.0029                                            (PA.34)

Initial porosity f0 of the Castile brine reservoir is calculated from the uncertain sampled parameter for the bulk Castile rock compressibility (BPCOMP; see Table PA-19), according to the following relationship:

                                                                                                          (PA.35)

where 1.0860 ´ 10-10 is a scaling constant that ensures that the productivity ratio, PR,  remains constant at  2.0 ´ 10-3 m3/Pa.  The productivity ratio PR is computed by

                                                                                                          (PA.36)

where V is the volume of the grid block representing the Castile brine reservoir in Figure PA-15. Because of this relationship, the initial porosity of the brine reservoir ranges from 0.1842 to 0.9208.  This range of porosity is not meant to represent an actual reservoir, but rather allows a reservoir to supply a volume of brine to the repository in the event of an E1 intrusion consistent with observed brine flows in the Delaware Basin.

The compressibility c f in Equation (PA.30) and Table PA-3 is pore compressibility.  Compressibility is treated as uncertain for Salado anhydrite, Salado halite, and regions of pressurized brine in the Castile.  However, the sampled value for each of these variables corresponds to bulk compressibility rather than to the pore compressibility actually used in the calculation.  Assuming all of the change in volume during compression occurs in the pore volume, the conversion from bulk compressibility cr to pore compressibility c f is approximated by

                                                                                                                          (PA.37)

where f0 is the initial porosity in the region under consideration.

The primary model used in PA for capillary pressure Pc and relative permeability krl is a modification of the Brooks-Corey model (Brooks and Corey 1964).  In this model, Pc, krb, and krg are defined by

                                                                                                          (PA.38)


Figure PA-15.  Computational Grid Used in BRAGFLO for PA


Figure PA-16.  Definition of Element Depth in BRAGFLO Grid


 

Figure PA-17.  Identification of Individual Cells in BRAGFLO Grid

Table PA-3.  Parameter Values Used in Representation of Two-Phase Flow

Region

Material

Material Description

Brooks-Corey Pore Distribution (PORE_DIS)a
l

Threshold Pressure Linear Parameter (PCT_A)a
a

Threshold Pressure Exponential Parameter
(PCT_EXP)a
h

Residual Brine Saturation
(SAT_RBRN)a
Sbr

Residual Gas Saturation
(SAT_RGAS)a
Sgr

Porosity
(POROSITY)a
f0

Pore Compressibilitya
c
f  , Pa-1

Intrinsic Permeability
(PRMX_LOG)a
k, m2

Salado

S_HALITE

Undisturbed halite

0.7

0.56

-0.346

0.3

0.2

HALPORb

f(HALCOMP)b,d

10x, x = HALPRMb

DRZ

DRZ_0

DRZ, -5 to 0 years

0.7

0.0

0.0

0.0

0.0

f(HALPOR)b,c

f(HALCOMP)b,d

1.0 ´ 10-17

 

DRZ_1

DRZ, 0 to 10,000 years

0.7

0.0

0.0

0.0

0.0

f(HALPOR)b,c

f(HALCOMP)b,d

10x, x = DRZPRMb

MB 138

S_MB138

Anhydrite MB in Salado

ANHBCEXPb

0.26

-0.348

ANRBSATb

ANRGSSATb

0.011

f(ANHCOMP)b,d

10x, x = ANHPRMb

Anhydrite AB

S_ANH_AB

Anhydrite layers A and B in Salado

ANHBCEXPb

0.26

-0.348

ANRBSATb

ANRGSSATb

0.011

f(ANHCOMP)b,d

10x, x = ANHPRMb

MB 139

S_MB139

Anhydrite MB in Salado

ANHBCEXPb

0.26

-0.348

ANRBSATb

ANRGSSATb

0.011

f(ANHCOMP)b,d

10x, x = ANHPRMb

Waste Panel

CAVITY_1

Single waste panel, -5 to 0 years

NAe

NAe

NAe

0.0

0.0

1.0

0.0

1.0 ´ 10-10

WAS_AREA

Single waste panel, 0 to 10,000 years

2.89

0.0

0.0

WRBRNSATb

WRGSSATb

0.848f

0.0

2.4 ´ 10-13

Rest of Repository (RoR)

CAVITY_2

RoR, -5 to 0 years

NAe

NAe

NAe

0.0

0.0

1.0

0.0

1.0 ´ 10-10

REPOSIT

RoR, 0 to 10,000 years

2.89

0.0

0.0

WRBRNSATb

WRGSSATb

0.848f

0.0

2.4 ´ 10-13

Ops and Exp

CAVITY_3

Operations area, -5 to 0 years

NAe

NAe

NAe

0.0

0.0

1.0

0.0

1.0 ´ 10-10

OPS_AREA

Operations area, 0 to 10,000 years

NAe

NAe

NAe

0.0

0.0

0.18

0.0

1.0 ´ 10-11

Exp

CAVITY_3

Experimental area, -5 to 0 years

NAe

NAe

NAe

0.0

0.0

1.0

0.0

1.0 ´ 10-10

a   Parenthetical parameter names are property names for the corresponding material, as indicated in Table PA-19.

b   Uncertain variable; see Table PA-19.

c   See Equation (PA.34).

d   See Equation (PA.37); f0 can also be defined by an uncertain variable.

e   These materials are using relative permeability model = 11; see Table PA-4.

f    Initial value of porosity f0; porosity changes dynamically to account for creep closure (see Section PA-4.2.3).

g   See Equation (PA.35).

 

Table PA-3.  Parameter Values Used in Representation of Two-Phase Flow (Continued)

Region

Material

Material Description

Brooks-Corey Pore Distribution (PORE_DIS)a
l

Threshold Pressure Linear Parameter (PCT_A)a
a

Threshold Pressure Exponential Parameter
(PCT_EXP)a
h

Residual Brine Saturation
(SAT_RBRN)a
Sbr

Residual Gas Saturation
(SAT_RGAS)a
Sgr

Porosity
(POROSITY)a
f0

Pore Compressibilitya
c
f  , Pa-1

Intrinsic Permeability
(PRMX_LOG)a
k, m2

EXP_AREA

Experimental area, 0 to 10,000 years

NAe

NAe

NAe

0.0

0.0

0.18

0.0

1.0 ´ 10-11

Castile

IMPERM_Z

Castile

0.7

0.0

0.0

0.0

0.0

0.005

0.0

1.0 ´ 10-35

Castile Brine Reservoir

CASTILER

Brine Reservoir in Castile

0.7

0.56

-0.346

0.2

0.2

f(BPCOMP)b,g

f(BPCOMP)b,d

10x, x = BPPRMb

Culebra

CULEBRA

Culebra Member of Rustler

0.6436

0.26

-0.348

0.08363

0.07711

0.151

6.622517 ´ 10-10

7.72681 ´ 10-14

Magenta

MAGENTA

Magenta Member of Rustler

0.6436

0.26

-0.348

0.08363

0.07711

0.138

1.915942 ´ 10-9

6.309576 ´ 10-16

Dewey Lake

DEWYLAKE

Dewey Lake Redbeds

0.6436

0.0

0.0

0.08363

0.07711

0.143

6.993007 ´ 10-8

5.011881 ´ 10-17

Santa Rosa

SANTAROS

Santa Rosa Formation

0.6436

0.0

0.0

0.08363

0.07711

0.175

5.714286 ´ 10-8

1.0 ´ 10-10

Los Medaños

UNNAMED

Los Medaños Member of Rustler

0.7

0.0

0.0

0.2

0.2

0.181

0.0

1.0 ´ 10-35

Tamarisk

TAMARISK

Tamarisk Member of Rustler

0.7

0.0

0.0

0.2

0.2

0.064

0.0

1.0 ´ 10-35

Forty-niner

FORTYNIN

Forty-niner Member of Rustler

0.7

0.0

0.0

0.2

0.2

0.082

0.0

1.0 ´ 10-35

DRZ_PCS

DRZ_0

DRZ, -5 to 0 years

0.7

0.0

0.0

0.0

0.0

f(HALPOR)b,c

f(HALCOMP)b,d

1.0 ´ 10-17

DRZ_PCS

DRZ_PCS

DRZ above the panel closures, 0 to 10,000 years

0.7

0.0

0.0

0.0

0.0

f(HALPOR)b,c

f(HALCOMP)b,d

10x, x = DRZPCPRMb

a   Parenthetical parameter names are property names for the corresponding material, as indicated in Table PA-19.

b   Uncertain variable; see Table PA-19.

c   See Equation (PA.34).

d   See Equation (PA.37); f0 can also be defined by an uncertain variable.

e   These materials are using relative permeability model = 11; see Table PA-4.

f    Initial value of porosity f0; porosity changes dynamically to account for creep closure (see Section PA-4.2.3).

g   See Equation (PA.35).

 

Table PA-3.  Parameter Values Used in Representation of Two-Phase Flow (Continued)

Region

Material

Material Description

Brooks-Corey Pore Distribution (PORE_DIS)a
l

Threshold Pressure Linear Parameter (PCT_A)a
a

Threshold Pressure Exponential Parameter
(PCT_EXP)a
h

Residual Brine Saturation
(SAT_RBRN)a
Sbr

Residual Gas Saturation
(SAT_RGAS)a
Sgr

Porosity
(POROSITY)a
f0

Pore Compressibilitya
c
f  , Pa-1

Intrinsic Permeability
(PRMX_LOG)a
k, m2

CONC_PCS

CAVITY_4

Concrete portion of panel closures, -5 to 0 years

NAe

NAe

NAe

0.0

0.0

1.0

0.0

1.0 ´ 10-10

CONC_PCS

Concrete portion of panel closures, 0 to 10,000 years

CONBCEXPb

0.0

0.0

CONBRSATb

CONGSSATb

0.05

1.2 ´ 10-9

10x, x = CONPRMb

DRF_PCS

CAVITY_4

Drift adjacent to panel closures, -5 to 0 years

NAe

NAe

NAe

0.0

0.0

1.0

0.0

1.0 ´ 10-10

DRF_PCS

Drift adjacent to panel closures, 0 to 10,000 years

2.89

0.0

0.0

WRBRNSATb

WRGSSATb

0.848f

0.0

2.4 ´ 10-13

CONC_MON

CAVITY_4

Concrete monolith portion of shaft seals, -5 to 0 years

NAe

NAe

NAe

0.0

0.0

1.0

0.0

1.0 ´ 10-10

CONC_MON

Concrete monolith portion of shaft seals, 0 to 10,000 years

0.94

0.0

0.0

SHURBRNb

SHURGASb

0.05

1.2 ´ 10-9

1.0 ´ 10-14

Upper Shaft

CAVITY_4

Upper portion of shaft seals, -5 to 0 years

NAe

NAe

NAe

0.0

0.0

1.0

0.0

1.0 ´ 10-10

SHFTU

Upper portion of shaft seals, 0 to 10,000 years

CONBCEXPb

0.0

0.0

SHURBRNb

SHURGASb

0.005

2.05 ´ 10-8

10x, x = SHUPRMb

a   Parenthetical parameter names are property names for the corresponding material, as indicated in Table PA-19.

b   Uncertain variable; see Table PA-19.

c   See Equation (PA.34).

d   See Equation (PA.37); f0 can also be defined by an uncertain variable.

e   These materials are using relative permeability model = 11; see Table PA-4.

f    Initial value of porosity f0; porosity changes dynamically to account for creep closure (see Section PA-4.2.3).

g   See Equation (PA.35).

 

Table PA-3.  Parameter Values Used in Representation of Two-Phase Flow (Continued)

Region

Material

Material Description

Brooks-Corey Pore Distribution (PORE_DIS)a
l

Threshold Pressure Linear Parameter (PCT_A)a
a

Threshold Pressure Exponential Parameter
(PCT_EXP)a
h

Residual Brine Saturation
(SAT_RBRN)a
Sbr

Residual Gas Saturation
(SAT_RGAS)a
Sgr

Porosity
(POROSITY)a
f0

Pore Compressibilitya
c
f  , Pa-1

Intrinsic Permeability
(PRMX_LOG)a
k, m2

Lower Shaft

CAVITY_4

Lower portion of shaft seals, -5 to 0 years

NAe

NAe

NAe

0.0

0.0

1.0

0.0

1.0 ´ 10-10

SHFTL_T1

Lower portion of shaft seals, 0 to 200 years

CONBCEXPb

0.0

0.0

SHURBRNb

SHURGASb

0.005

4.28 ´ 10-9

10x, x = SHLPRM1b

SHFTL_T2

Lower portion of shaft seals, 200 to 10,000 years

CONBCEXPb

0.0

0.0

SHURBRNb

SHURGASb

0.005

4.28 ´ 10-9

10x, x = SHLPRM2b

Borehole plugs

CONC_PLG

Concrete borehole plug, before plug degradation

0.94

0.0

0.0

0.0

0.0

0.32

1.1875 ´ 10-9

10x, x = PLGPRMb

BH_SAND

Borehole after plug degradation, 200 years after intrusion

0.94

0.0

0.0

0.0

0.0

0.32

0.0

10x, x = BHPRMb

Upper Borehole

BH_OPEN

Borehole above repository before plug degradation

0.7

0.0

0.0

0.0

0.0

0.32

0.0

1.0 ´ 10-9

BH_SAND

Borehole after plug degradation, 200 years after intrusion

0.94

0.0

0.0

0.0

0.0

0.32

0.0

10x, x = BHPRMb

Lower Borehole

BH_OPEN

Borehole below repository before creep closure

0.7

0.0

0.0

0.0

0.0

0.32

0.0

1.0 ´ 10-9

BH_CREEP

Borehole below repository after creep closure, 1,000 years after intrusion

0.94

0.0

0.0

0.0

0.0

0.32

0.0

10x/10, x = BHPRMa

a   Parenthetical parameter names are property names for the corresponding material, as indicated in Table PA-19.

b   Uncertain variable; see Table PA-19.

c   See Equation (PA.34).

d   See Equation (PA.37); f0 can also be defined by an uncertain variable.

e   These materials are using relative permeability model = 11; see Table PA-4.

f    Initial value of porosity f0; porosity changes dynamically to account for creep closure (see Section PA-4.2.3).

g   See Equation (PA.35).


                                                                                                              (PA.39)

                                                                                     (PA.40)

where

l       =  pore distribution parameter (dimensionless)

Pt( k) =  capillary threshold pressure (Pa) as a function of intrinsic permeability k (Webb 1992)

                                                                                                              (PA.41)

   =  effective brine saturation (dimensionless) without correction for residual gas saturation

                                                                                     (PA.42)

  =  effective brine saturation (dimensionless) with correction for residual gas saturation

                                                                           (PA.43)

The values used for l, a, h, Sbr, Sgr, and k are summarized in Table PA-3.  The statement that the Brooks-Corey model is in use means that Pc, krb, and krg are defined by Equation (PA.38) through Equation (PA.40).

In the anhydrite MBs, either the Brooks-Corey model or the van Genuchten-Parker model is used as determined by the subjectively uncertain parameter ANHBCVGP (see Table PA-19).  A linear model is used to represent two-phase flow in an open borehole (i.e., for the first 200 years after a drilling intrusion for boreholes with two-plug or three-plug configurations, in the open cavities [CAVITY_1, . . , CAVITY_4], and for the experimental and operations areas).  This is discussed further below.

In the van Genuchten-Parker model, Pc, krb, and krg are defined by (van Genuchten 1978)

                                                                                               (PA.44)

                                                                                        (PA.45)

                                                                                     (PA.46)

where m = l/(1 + l) and the capillary pressure parameter PVGP is determined by requiring that the capillary pressures defined in Equation (PA.38) and Equation (PA.44) are equal at an effective brine saturation of Se2 = 0.5 (Webb 1992).  The van Genuchten-Parker model is only used for the anhydrite MBs in the Salado and uses the same values for l, Sbr, and Sgr as the Brooks-Corey model (Table PA-3).

In the linear model used for the open borehole (RELP_MOD = 5), Pc, krb, and krg are defined by

                                                        Pc = 0, krb = Se1, krg = 1 – Se1                                     (PA.47)

Another linear model (RELP_MOD = 11) is used for the open cavities (CAVITY_1, . . . , CAVITY_4) for the −5 to 0 year portion of the simulation (see Section PA-4.2.2) and the experimental and operations areas (t = 0 to 10,000 years) which, in PA, are modeled without a time-dependent creep closure:

                                                                                                     (PA.48)

                                                                         (PA.49)

                                                                                             (PA.50)

where l = gas or brine and tol is a tolerance (slope) over which the relative permeability changes linearly from 0 to 1.  In PA, tol = 1 ´ 10-2 (dimensionless).  Thus, the relative permeabilities are ~ 1 for saturations away from residual saturation.

Capillary pressure Pc for both the van Genuchten-Parker and Brooks-Corey models becomes unbounded as brine saturation Sb approaches the residual brine saturation, Sbr.  To avoid unbounded values, Pc is capped at 1 ´ 108 Pa in selected regions (Table PA-4).

Gas density is computed using the RKS equation of state, with the gas assumed to be pure H2.  For a pure gas, the RKS equation of state has the form (Walas 1985, pp. 43−54)

                                                                                                    (PA.51)

where

             R  =  gas constant = 8.31451 Joules (J) mole (mol) -1 K -1

             T  =  temperature (K) = 300.15 K (= 30 °C; 81 °F)

             V  =  molar volume (m3 mol -1)

             a  =  0.42747/Pcrit

             b  =  0.08664 RTcrit  /Pcrit

Table PA-4.  Models for Relative Permeability and Capillary Pressure in Two-Phase Flow

Material

Relative Permeabilitya
(RELP_MOD)

Capillary Pressureb
(CAP_MOD)

Material

Relative Permeabilitya
(RELP_MOD)

Capillary Pressureb
(CAP_MOD)

S_HALITE

4

2

WAS_AREA

12

1

DRZ_0

4

1

DRZ_1

4

1

S_MB139

ANHBCVGPc

2

DRZ_PCS

Sampled

1

S_ANH_AB

ANHBCVGPc

2

CONC_PCS

4

2

S_MB138

ANHBCVGPc

2

UNNAMED

4

1

CAVITY_1

11

1

TAMARISK

4

1

CAVITY_2

11

1

FORTYNIN

4

1

CAVITY_3

11

1

DRF_PCS

12

1

CAVITY_4

11

1

REPOSIT

12

1

IMPERM_Z 

4

1

CONC_MON

4

2

CASTILER

4

2

SHFTU

4

1

OPS_AREA

11

1

SHFTL_T1

4

1

EXP_AREA

11

1

SHFTL_T2

4

1

CULEBRA

4

2

CONC_PLG

4

1

MAGENTA

4

2

BH_OPEN

5

1

DEWYLAKE

4

2

BH_SAND

4

1

SANTAROS

4

1

BH_CREEP

4

1

a   Relative permeability model, where 4 = Brooks-Corey model given by Equation (PA.38) through Equation (PA.40), 5 = linear model given by Equation (PA.47), 11 = linear model given by Equation (PA.48) through Equation (PA.50), 12 = modified Brooks-Corey model to account for cutoff saturation (Nemer 2007c),  and ANHBCVGP ~ use of Brooks-Corey or van Genuchten-Parker model treated as a subjective uncertainty.

b    Capillary pressure model, where 1 = capillary pressure is unbounded, 2 = Pc bounded above by 1 ´ 108 Pa as Sb approaches Sbr.

c    See ANHBCVGP in Table PA-19.

 

            a  = 

                 »  for H2 (Graboski and Daubert 1979)

          Tcrit  =  critical temperature (K)

          Pcrit  =  critical pressure (Pa)

            Tr  =  T / Tcrit = reduced temperature

            w  =  acentric factor

                 =  0 for H2 (Graboski and Daubert 1979)

In order to account for quantum effects in H2, effective critical temperature and pressure values of Tcrit  = 43.6 K and Pcrit = 2.047 ´ 106 Pa are used instead of the true values for these properties (Prausnitz 1969).  Equation (PA.51) is solved for molar volume V.  The gas density rg then is given by

                                                                                                                  (PA.52)

where Mw,H 2 is the molecular weight of H2 (i.e., 2.01588 ´ 10 -3 kg/mol; see Weast 1969, p. B-26).

Brine density rb is defined by Equation (PA.29), with rb0= 1230.0 kg/m3 at a pressure of Pb0 = 1.0132 ´ 105  Pa and cb = 2.5 ´ 10 -10 Pa -1 (Roberts 1996).  Porosity, f, is used as defined by Equation (PA.30) with two exceptions:  in the repository (see Section PA-4.2.3) and in the DRZ and MBs subsequent to fracturing (see Section PA-4.2.4).  The values of f0 and c f used in conjunction with Equation (PA.30) are listed in Table PA-3.  The reference pressure Pb0 in Equation (PA.30) is spatially variable and corresponds to the initial pressures Pb( x, y, −5) (here, −5 means at time equal to −5 years; see Section PA-4.2.2). The gas and brine viscosities ml, g, b in Equation (PA.24) and Equation (PA.25) were assumed to have values of mg = 8.93 ´ 10 -6 Pa s (H2:VISCO; see Vargaftik 1975) and mb = 2.1 ´ 10 -3 Pa s (BRINESAL:VISCO; see McTigue 1993).

The terms qg, qrg, qb, and qrb in Equation (PA.24) and Equation (PA.25) relate to well injection or removal (i.e., qg, qb) and reaction, production, or consumption (i.e., qrg, qrb) of gas and brine, with positive signs corresponding to injection or production and negative signs corresponding to removal or consumption.  In the long-term Salado flow calculations, no injection or removal of gas or brine is calculated using qg and qb.  Thus, qg and qb are equal to zero.  That is, after an intrusion, the borehole is treated as a porous media, rather than a point source or sink of brine and gas.  Furthermore, the mass and pressure lost to a DBR during the intrusion is conservatively ignored in the BRAGFLO calculations.  In the DBR calculations discussed in Section PA-4.7, qg and qb are used to describe injection and production wells in the DBR grid.

In PA, no gas consumption occurs through the term qrg (see below), and gas production has the potential to occur (due to corrosion of steel or microbial degradation of CPR materials) only in the waste disposal regions of the repository (i.e., Waste Panel, South RoR, and North RoR in Figure PA-15).  Thus,

            qrg ³ 0  in waste disposal regions of Figure PA-15

= 0  elsewhere                                                                                                (PA.53)

Gas consumption occurs due to the reaction of CO2 with MgO in the waste panels, and potentially from the sulfidation of steel.  This gas consumption is not modeled using qrg, but is accounted for by reducing the gas generation rate qrg, as discussed in Section PA-4.2.5.  Finally, no brine production occurs, and brine consumption has the potential to occur (due to the consumption of brine during steel corrosion) only in the waste disposal regions of the repository.  Thus,

            qrb   £ 0  in waste disposal regions of Figure PA-15

= 0  elsewhere                                                                                                (PA.54)

More detail on the definition of qrg and qrb is provided in Section PA-4.2.5.

In each two-phase flow simulation, a short period of time representing disposal operations is simulated.  This period of time is called the start-up period, and covers 5 years from t = -5 years to 0 years, corresponding to the amount of time a typical panel is expected to be open during disposal operations.  All grid locations require initial brine pressure and gas saturation at the beginning of the simulation (t = -5 years).

The Rustler and overlying units (except in the shaft) are modeled as horizontal with spatially constant initial pressure in each layer (see Figure PA-15).  Table PA-5 lists the initial brine pressure, Pb, and gas saturation, Sg, for the Rustler.

The Salado (Mesh Rows 3–24 in Figure PA-15) is assumed to dip uniformly q = 1 degree downward from north to south (right to left in Figure PA-15).  Except in the repository excavations and the shaft, brine is initially assumed (i.e., at -5 years) to be in hydrostatic equilibrium relative to an uncertain initial pressure Pb,ref (SALPRES; see Table PA-19) at a reference point located at shaft center at the elevation of the midpoint of MB 139, which is the center of Cell 1266 in Figure PA-17.  This gives rise to the condition

                                                                          (PA.55)

                                           (PA.56)

                               (PA.57)

                                                          (PA.58)

                                                                                                          (PA.59)

Table PA-5.  Initial Conditions in the Rustler

Name

Mesh Row
(Figure PA-15)

Pb (x, y, -5), Pa

Sg (x, y, -5)

Santa Rosa c

33

1.013250 ´ 105

1 - Sb = 0.916
(Sb = SANTAROS:SAT_IBRN)a

Santa Rosa c

32

1.013250 ´ 105

1 - Sb  = 0.916
(Sb = SANTAROS:SAT_IBRN)a

Dewey Lakec

31

1.013250 ´ 105

1 - Sb = 0.916
(Sb = SANTAROS:SAT_USAT)a

Dewey Lakec

30

7.355092 ´ 105

1 - Sb = 0.916
(Sb = SANTAROS:SAT_USAT)a

Forty-ninerc

29

1.47328 ´ 106

0b

Magenta

28

9.465 ´ 105
(MAGENTA:PRESSURE)

0b

Tamariskc

27

1.82709 ´ 106

0b

Culebra

26

9.141 ´ 105

(CULEBRA:PRESSURE)

0b

Los Medaños c

25

2.28346 ´ 106

0b

a    The names in parenthesis are parameters in the WIPP PA Parameter Database.

b    The Rustler is assumed to be fully saturated.  This initial condition is set in the program ICSET.  See Nemer and Clayton (2008, Section 3.2).

c    These pressures are calculated in the ALGEBRA1 step analogously to Equation (PA.55) below, using the brine density of 1220 kg/m3.  See subsequent discussion taking θ = 0 and the reference point (xref, yref) at the top of the Dewey Lake.  See the ALGEBRA input file ALG1_BF_CRA09.INP in library LIBCRA09_BF, class CRA09-1 on the WIPP PA cluster for details.  See Nemer and Clayton (2008, Section 4.1.7) for details on the ALGEBRA1 step.

 

where

     h(x, y)      is defined in Equation (PA.33)

          rb0  =  1220 kg/m3 (BRINESAL:DNSFLUID)

            cb  =  3.1 ´ 10 -10 Pa -1 (BRINESAL:COMPRES)

             g =  9.80665 m/s2

        Pb,ref  =  1.01325 ´ 105 Pa (BRINESAL:REF_PRES)

          Pb0  =  sampled far-field pressure in the undisturbed halite (S_HALITE:PRESSURE)

In the Salado, initial gas saturation Sg( x, y, -5) = 0 (see Nemer and Clayton 2008, Section 4.1.6).

The Castile (Mesh Rows 1 and 2) is modeled as horizontal, and initial brine pressure is spatially constant within each layer (no dip), except that the brine reservoir is treated as a different material from rest of Castile and has a different initial pressure which is a sampled parameter.  Specifically, outside the brine reservoir, pressure is calculated using Equation (PA.55) above with no dip (q = 0) in the ALGEBRA1 step.  Within the reservoir, Pb( x, y, -5) = BPINTPRS, the uncertain initial pressure in the reservoir (see Table PA-19).  Initial gas saturation Sg( x, y, -5) = 0.

Within the shaft (areas Upper Shaft, Lower Shaft, and CONC_MON) and panel closures (areas CONC_PCS and DRF_PCS), Pb( x, y, -5) = 1.01325 ´ 105 Pa and Sg( x, y, -5) = 0.  Within the excavated areas (Waste Panel, South RoR, and North RoR, Ops and Exp), Pb( x, y, -5) = 1.01325 ´ 105 Pa and Sg( x, y, -5) = 0.

At the end of the initial five-year start-up period and the beginning of the regulatory period (t = 0 years), brine pressure and gas saturation are reset in the shaft, panel closures, and excavated areas.  In the shaft (areas Upper Shaft, Lower Shaft, and CONC_MON) and panel closures (areas CONC_PCS and DRF_PCS), Pb( x, y, 0) = 1.01325 ´ 105 Pa and Sg( x, y, 0) = 1 ´ 10 -7 (see CONC_MON:SAT_IBRN).  In the waste disposal regions (areas Waste Panel, South RoR, and North RoR), Pb( x, y, 0) = 1.01325 ´ 105 Pa and Sg( x, y, 05) = 0.985 (see WAS_AREA:SAT_IBRN).  In the other excavated areas, Pb( x, y, 0) = 1.01325 ´ 105 Paand Sg( xy, 0) = 1.0.

Salt creep occurs naturally in the Salado halite in response to deviatoric stress.  Inward creep of rock is generally referred to as creep closure.  Creep closure of excavated regions begins immediately from excavation-induced deviatoric stress.  If the rooms were empty, closure would proceed to the point where the void volume created by the excavation would be eliminated as the surrounding formation returned to a uniform stress state.  In the waste disposal region, inward creep of salt causes consolidation of the waste, and this waste consolidation continues until loading in the surrounding rock is uniform, at which point salt creep and waste consolidation ceases.  The amount of waste consolidation that occurs and the time it takes to consolidate are governed by the waste properties (e.g., waste strength, modulus, etc.), the surrounding rock properties, the dimensions and location of the room, and relative quantities of brine and gas present.

The porosity of the waste disposal regions and neighboring access drifts (i.e., Waste Panel, South RoR, North RoR, and DRF_PCS in Figure PA-15) is assumed to change through time due to creep closure of the halite surrounding the excavations.  The equations on which BRAGFLO is based do not incorporate this type of deformation.  Therefore, the changes in repository porosity due to halite deformation are modeled in a separate analysis with the geomechanical program SANTOS, which implements a quasi-static, large-deformation, finite-element procedure (Stone 1997).  Interpolation procedures are then used with the SANTOS results to define porosity (f) within the repository as a function of time, pressure, and gas generation rate.

For more information on the generation of the porosity surface for BRAGFLO in PA, see Appendix PORSURF-2009.

Fracturing within the anhydrite MBs (i.e., regions MB 138, Anhydrite AB, and MB 139 in Figure PA-15) and in the DRZ (region DRZ in Figure PA-15) is assumed to occur at pressures slightly above lithostatic pressure, and is implemented through a pressure-dependent compressibility cr( Pb) (Mendenhall and Gerstle 1995).  Specifically, MB fracturing begins at a brine pressure of

                                                                                                                (PA.60)

where Pbi and Pb0 are spatially dependent (i.e., Pb0 = P(x, y, 0) as in Section PA-4.2.2) and DPi = 2 ´ 105 Pa (see S_MB138:PI_DELTA in Fox 2008, Table 30)

Fracturing ceases at a pressure of

                                                                                                               (PA.61)

and a fully fractured porosity of

                                                                                                   (PA.62)

where DPa = 3.8 ´ 106 Pa (see S_MB138:PF_DELTA in Fox 2008, Table 30), f0 is spatially dependent (Table PA-3), and Dfa = 0.04, 0.24, and 0.04 for anhydrite materials S_MB138, S_ANH_AB, and S_MB139, respectively (see S_MB138:DPHIMAX in Fox 2008, Table 30).

Once fractured, compressibility cr becomes a linear function

                                                                             (PA.63)

of brine pressure for Pbi £ Pb £ Pba, with cra defined so that the solution f of

                                                (PA.64)

satisfies f (Pba) = fa; specifically, cra is given by

                                                       (PA.65)

The permeability kf( Pb) of fractured material at brine pressure Pb is related to the permeability of unfractured material at brine pressure Pbi by

                                                                                                     (PA.66)

where k is the permeability of unfractured material (i.e., at Pbi) and n is defined so that kf( Pba) = 1 ´ 10 -9 m2 (i.e., n is a function of k, which is an uncertain input to the analysis; see ANHPRM in Table PA-19).  When fracturing occurs, kf( Pb) is used instead of k in the definition of the permeability for the fractured areas of the anhydrite MBs.

Fracturing is also modeled in the DRZ region in Figure PA-15.  The fracture model implementation is the same as for the anhydrite materials.  In this case, fracturing would be in halite rather than anhydrite, but because of the limited extent of the DRZ and the proximity of the nearby interbeds, this representation was deemed acceptable by the Salado Flow Peer Review panel (Caporuscio, Gibbons, and Oswald 2003).

Gas production is assumed to result from anoxic corrosion of steel and the microbial degradation of CPR materials.  Thus, the gas generation rate qrg in Equation (PA.24) is of the form

                                                                                                            (PA.67)

where qrgc is the rate of gas production per unit volume of waste (kg/m3/s) due to anoxic corrosion of Fe-base metals, and qrgm is the rate of gas production per unit volume of waste (kg/m3/s) due to microbial degradation of CPR materials.  Furthermore, qrb in Equation (PA.25) is used to describe the consumption of brine during the corrosion process.

Gas generation takes place only within the waste disposal regions (i.e., Waste Panel, South RoR, and North RoR in Figure PA-15) and all the generated gas is assumed to have the same properties as H2 (see discussion in Appendix MASS-2009, Section MASS-3.2).  In PA, the consumable materials are assumed to be homogeneously distributed throughout the waste disposal regions (i.e., the concentration of Fe-base metals and CPR materials in the waste is not spatially dependent).  A separate analysis examined the potential effects on PA results of spatially varying Fe-base metal and CPR material concentrations, and concluded that PA results are not affected by representing these materials with spatially varying concentrations (see Appendix MASS-2009, Section MASS-21.0).

The rates qrgc, qrb, and qrgm (kg/m3/s) are defined by

gas generation by corrosion

                                                       (PA.68)

brine consumption by corrosion

                                                                       (PA.69)

and microbial gas generation

                                                       (PA.70)

where

                          Ds    = surface area concentration of steel in the repository (m2  surface area steel/ m3 disposal volume)

                          Dc  =  mass concentration of cellulosics in the repository (kg biodegradable material/m3 disposal volume)

                       =  molecular weight of H2 (kg H2/mol H2)

                    =  molecular weight of water (H2O) (kg H2O/mol H2O)

                          Rci  =  corrosion rate under inundated conditions (m/s)

                         Rch  =  corrosion rate under humid conditions (m/s)

                         Rmi  =  rate of cellulose biodegradation under inundated conditions (mol C6H10O5/kg C6H10O5/s)

                         Rmh  =  rate of cellulose biodegradation under humid conditions (mol C6H10O5/kg C6H10O5/s)

                        Sb,eff  =  effective brine saturation due to capillary action in the waste materials (see Equation (PA.90) in Section PA-4.2.6)

                          = 

            =  stoichiometric coefficient for gas generation due to corrosion of steel, i.e., moles of H2 produced by the corrosion of 1 mole of Fe (mol H2/mol Fe)

        =  stoichiometric coefficient for brine consumption due to corrosion of steel, i.e., moles of H2O consumed per mole of H2 generated by corrosion (mol H2O/mol H2)

               =  average stoichiometric factor for microbial degradation of cellulose, i.e., the moles of H2 generated per mole of carbon consumed by microbial action (mol H2/mol C)

                         rFe  =  molar density of steel (mol/m3)

                          Bfc  =  parameter (WAS_AREA:BIOGENFC; discussed in detail later in this section) uniformly sampled from 0 to 1, used to account for the uncertainty in whether microbial gas generation could be realized in the WIPP at experimentally measured rates.

The products Rci  Ds  rFe Xc , Rch  Ds  rFe Xc, Rmi  Dc  y, and Rmh in Equation (PA.68) and Equation (PA.70) define constant rates of gas generation (mol/m3/s) that continue until the associated substrate (i.e., steel or cellulose) is exhausted (i.e., zero order kinetics are assumed).  The terms Sb,eff and  in Equation (PA.68) and Equation (PA.70), which are functions of location and time, correct for the amount of substrate exposed to inundated and humid conditions, respectively.  All the corrosion and microbial action is assumed to cease when no brine is present, which is the reason that 0 replaces Sg = 1 in the definition of .  In PA, Rch = 0 and Rci, Rmh, and Rmi are defined by uncertain variables (see WGRCOR, WGRMICH, WGRMICI in Table PA-19).  However, Rmh is now sampled based on the sampled value of Rmi: see Nemer and Clayton (2008, Section 5.1.3).  Further, MH 2 = 2.02 ´10 -3 kg/mol (Lide 1991, pp. 1-7, 1-8), MH 2O = 1.80 ´ 10 -2 kg/mol (Lide 1991, pp. 1-7, 1-8), rFe = 1.41 ´ 105 mol/m3 (Telander and Westerman 1993), and Ds, Dc, Xc(H2O|H2), Xc(H2|Fe), and y(H2| C) are discussed below.

The concentration Ds in Equation (PA.68) is defined by

                                                                                                               (PA.71)

where

         Ad  =  surface area of steel associated with a waste disposal drum (m2/drum)

        VR  =  initial volume of a single room in the repository (m3)

         nd  =  ideal number of waste drums that can be close-packed into a single room

In PA, Ad = 6 m2/drum (REFCON:ASDRUM), VR = 3,644 m3 (REFCON:VROOM), and nd = 6804 drums (REFCON:DRROOM).

The biodegradable materials to be disposed at the WIPP consist of cellulosic materials, plastics, and rubbers. Cellulosics have been demonstrated experimentally to be the most biodegradable of these materials (Francis, Gillow, and Giles 1997).  The occurrence of significant microbial gas generation in the repository will depend on whether (1) microbes capable of consuming the emplaced organic materials will be present and active, (2) sufficient electron acceptors will be present and available, and (3) enough nutrients will be present and available.

During the CRA-2004 PABC, the EPA (Cotsworth 2005) indicated that the probability that microbial gas generation could occur (WMICDFLG) should be set equal to 1 in PA calculations.  In the CRA-2004, the probability that microbial gas generation could occur was assigned a value of 0.5.  To comply with the EPA’s letter, in the CRA-2004 PABC and the CRA-2009 PA the parameter WMICDFLG was changed so that the probability that microbial gas generation could occur was set to 1 while preserving the previous probability distribution on whether CPR could be degraded.  This is summarized in Table PA-6, and is discussed further in Nemer and Stein

Table PA-6.    Probabilities for Biodegradation of Different Organic Materials (WAS_AREA:PROBDEG) in the CRA-2009 PA and the CRA-2004 PA

WAS_AREA:PROBDEG

Meaning

Probability CRA-2004

Probability CRA-2009

0

No microbial degradation can occur

0.5

0.0

1

Biodegradation of only cellulose can occur

0.25

0.75

2

Biodegradation of CPR materials can occur

0.25

0.25

 

(2005, Section 5.4).  Because there are significant uncertainties in whether the experimentally observed gas-generation rates could be realized in the WIPP repository, during the CRA-2004 PABC the EPA agreed to allow the DOE to multiply the sampled microbial rates by a parameter (WAS_AREA:BIOGENFC) uniformly sampled from 0 to 1.  This is discussed further in Nemer, Stein, and Zelinski (2005, Section 4.2.2).

In cases where biodegradation of rubbers and plastics occur, rubbers and plastics are converted to an equivalent quantity of cellulosics based on their carbon equivalence (Wang and Brush 1996a). This produces the density calculation

for biodegradation of cellulosics only

(PA.72)

for biodegradation of CPR materials

 

where mcel is the mass of cellulosics (kg), mr is the mass of rubbers (kg), and mp is the mass of plastics (kg).

In the CRA-2009 PA, the emplacement materials (cellulose and plastic) have been added to mcel and mp.  This is discussed further in Nemer and Clayton (2008, Section 5.1.1).  Density values for CPR materials can be found in Fox (2008, Table 34).

The most plausible corrosion reactions after closure of the WIPP are believed to be (Wang and Brush 1996a)

                                                        Fe + 2H2O = Fe(OH)2  + H2                                     (PA.73)

and

                                                         3Fe + 4H2O = Fe3O4 + 4H2                                     (PA.74)

When normalized to 1 mole of Fe and linearly weighted by the factors x and , the two preceding reactions become

                                     (PA.75)

where x and  are the fractions of Fe consumed in the reactions in Equation (PA.73) and Equation (PA.74), respectively.  Although magnetite (Fe3O4) has been observed to form on Fe as a corrosion product in low-Mg anoxic brines at elevated temperatures (Telander and Westerman 1997) and in oxic brine (Haberman and Frydrych 1988), there is no evidence that it will form at WIPP repository temperatures.  If Fe3O4 were to form, H2 would be produced (on a molar basis) in excess of the amount of Fe consumed. However, anoxic corrosion experiments (Telander and Westerman 1993) did not indicate the production of H2 in excess of the amount of Fe consumed.  Therefore, the stoichiometric factor x in Reaction (PA.75) is set to 1.0 (i.e., x = 1), which implies that Reaction (PA.73) represents corrosion. Thus, the stoichiometric factor for corrosion is

                                                                            (PA.76)

which implies that one mole of H2 is produced for each mole of Fe consumed, and the stoichiometric factor for brine consumption is

                                                                     (PA.77)

which implies that two moles of H2O are consumed for each mole of H2 produced.

The most plausible biodegradation reactions after closure of the WIPP are believed to be (Wang and Brush 1996a)

denitrification             C6H1OO5 + 4.8H+ + 4.8NO3 - = 7.4H2O + 6CO2 + 2.4N2                (PA.78)

sulfate reduction              C6H1OO5 + 6H+ + 3SO42 - = 5H2O + 6CO2 + 3H2S                    (PA.79)

methanogenesis                            C6H1OO5 + H2O = 3CH4 + 3CO2                                  (PA.80)

Accumulation of CO2 produced by the above reactions could decrease pH and add carbonate to increase An solubility in the repository (Wang and Brush 1996b).

However, in the CRA-2004 PABC, the EPA (Cotsworth 2005) directed the DOE to remove methanogenesis (Equation (PA.80)) from PA.  The EPA cited the presence of calcium sulfate as gypsum and anhydrite in the bedded salt surrounding the repository as possible sources of sulfate.  These sources of sulfate would, if accessible, promote sulfate reduction (Equation (PA.82)), which is energetically and kinetically favored over methanogenesis.  In response, the DOE removed methanogenesis from PA.  Additionally, the DOE removed Fe sulfidation from PA as a gas-consuming reaction because sulfidation of steel produces one mole of H2 for every mole of H2S consumed,

                                                              Fe + H2S = FeS + H2                                          (PA.81)

This and the removal of methanogenesis are discussed fully in Nemer and Zelinski (2005).

To provide added assurance of WIPP performance, a sufficient amount of MgO is added to the repository to remove CO2 (Bynum et al. 1997).  MgO in polypropylene “supersacks” is emplaced on top of the three-layer waste stacks to create conditions that reduce An solubilities in the repository (see Appendix MgO-2009 and Appendix SOTERM-2009, Section SOTERM-2.3).  If brine flows into the repository, MgO will react with water in brine and in the gaseous phase to produce brucite (Mg[OH]2).  MgO will react with essentially all of the CO2 that could be produced by complete microbial consumption of the CPR materials in the waste, and will create hydromagnesite with the composition Mg5(CO3)4(OH)2·4H2O (Appendix MgO-2009; Appendix SOTERM-2009, Section SOTERM-2.0).  The most important MgO hydration and carbonation reactions that will occur in the WIPP are

                                                MgO + H2O(aq and/or g) → Mg(OH)2                            (PA.82)

                                5Mg(OH)2 + 4CO2(aq and/or g) → Mg5(CO3)4(OH)2×4H2O            (PA.83)

In these equations, “aq and/or g” indicates that the H2O or CO2 that reacts with MgO and/or brucite could be present in the aqueous phase (brine) and/or the gaseous phase.  The removal of CO2 by MgO is not explicitly modeled as a separate reaction in PA.  Rather, the effect of CO2 consumption is accounted for by modifying the stoichiometry of Reaction (PA.78), Reaction (PA.79), and Reaction (PA.80) to remove the CO2 from the mass of gas produced by microbial action.

The average stoichiometry of Reaction (PA.78) , Reaction (PA.79), and Reaction (PA.80), is

                                                      (PA.84)

where the average stoichiometric factor y in Reaction (PA.84) represents the number of moles of gas produced and retained in the repository from each mole of carbon consumed.  This factor y depends on the extent of the individual biodegradation pathways in Reaction (PA.84), and the consumption of CO2 by MgO.  A range of values for y is estimated by considering the maximum mass of gas that can be produced from consumption of cellulosics (Mcel) and Fe-base metals (MFe), and is derived as follows (Wang and Brush 1996b).  Estimates of the maximum quantities Mcel and MFe (mol) of cellulosics (i.e., C6H10O5) and steels that can be potentially consumed in 10,000 years are given by

                                                                  (PA.85)

                                                                (PA.86)

where mcel and mFe are the masses (kg) of cellulosics (see Equation (PA.72) for definition) and steels initially present in the repository.  The mass of cellulosics that can be consumed is determined by the uncertain parameter WMICDFLG (see Table PA-19).  The mass of steels mFe = 5.16 ´ 107 kg; this value is calculated as

                                                              (PA.87)

where VCH and VRH are the volumes of CH-TRU and RH-TRU waste, rWCH  and rWRH are the Fe densities in CH-TRU and RH-TRU waste, and rCCH and rCRH  are the Fe densities of the containers of CH-TRU and RH-TRU waste (see Fox 2008, Table 34).  The terms 6000 mcel/162 and 1000 mFe/56 in Equation (PA.85) and Equation (PA.86) equal the inventories in moles of cellulosics and steel, respectively.  The terms 3.2 ´ 1011  Rm mcel and 4.4 ´ 1016  Rci Ad nd equal the maximum amounts of cellulosics and steel that could be consumed over 10,000 years.  In Equation (PA.85), Rm = max{Rmh, Rmi}, where Rmh and Rmi are defined by uncertain variables (see WGRMICH and WGRMICI in Table PA-19, respectively), and 3.2 ´ 1011 = (3.15569 ´ 107 s/yr) (104 yr).  In Equation (PA.86), Ad nd is the total surface area of all drums (m2) and the factor 4.4 ´ 1016 = (3.15569 ´ 107 s/yr) (104 yr) (1.41 ´ 105 mol/ m3), where rFe = 1.41 ´ 105 mol/m3 (see Equation (PA.72)) (Telander and Westerman 1993), converts the corrosion rate from m/s to mol/m2/s.

A range of possible values for the average stoichiometric factor y in Reaction (PA.84) can be obtained by considering individual biodegradation pathways involving Mcel and accounting for the removal of CO2 by the MgO. 

In the absence of methanogenesis and steel sulfidation, y from Equation (PA.84) becomes

                                                                           (PA.88)

                                                                                           (PA.89)

where  is the quantity of NO3 - (mols) initially present in the repository.  Specifically, = 4.31 ´ 107 mol (Fox 2008, Table 39).

Capillary action (wicking) is the ability of a material to carry a fluid by capillary forces above the level it would normally seek in response to gravity.  In the current analysis, this phenomena is accounted for by defining an effective saturation given by


                                                                                                                                         (
PA.90)

where

      Sb,eff  =  effective brine saturation

         Sb  =  brine saturation

      Swick  =  wicking saturation

      Smin    =  minimum brine saturation at which code can run in the waste-filled areas

          α  =  smoothing parameter = -1000

The effective saturation given by Equation (PA.90) differs from that in the CRA-2004 PA in that Sb,eff  now approaches zero as Sb approaches a small value Smin.  In simulations where Fe corrosion dried out the repository, the time required to complete the simulation could be quite long.  In order to speed up the code and increase robustness, the parameter Smin was added.  For PA, Smin = 0.015, which was small enough to not affect the results, while greatly reducing run time.  This is explained fully in Nemer and Clayton (2008, Section 5.2.2).

The effective saturation is used on a grid block basis within all waste regions (Waste Panel, South RoR, and North RoR in Figure PA-15).  The wicking saturation, Swick, is treated as an uncertain variable (see WASTWICK in Table PA-19).  The effective brine saturation Sb,eff is currently used only to calculate the corrosion of steel (Equation (PA.68)) and the microbial degradation of cellulose (Equation (PA.70)), and does not directly affect the two-phase flow calculations indicated.

The WIPP excavation includes four shafts that connect the repository region to the surface: the air intake shaft, salt handling shaft, waste handling shaft, and exhaust shaft. In PA, these four shafts are modeled as a single shaft.  The rationale for this modeling treatment is set forth in SNL (Sandia National Laboratories 1992, Volume 5, Section 2.3).

The shaft seal model included in the PA grid (Column 43 in Figure PA-15) is the simplified shaft model used in the CRA-2004 PA and the CRA-2004 PABC.  The simplified shaft seal model used in PA is described by Stein and Zelinski (2003) and is briefly discussed below; this model was approved by the Salado Flow Peer Review Panel (Caporuscio, Gibbons, and Oswald 2003).

The planned design of the shaft seals involves numerous materials, including earth, crushed salt, clay, asphalt, and Salado Mass Concrete (SMC) (see the CCA, Appendix SEAL).  The design is intended to control both short-term and long-term fluid flow through the Salado portion of the shafts.  For the CCA PA, each material in the shaft seal was represented in the BRAGFLO grid.  Analysis of the flow results from the CCA PA and the subsequent CCA PAVT (Sandia National Laboratories 1997 and U.S. Department of Energy 1997) indicated that no significant flows of brine or gas occurred in the shaft during the 10,000-year regulatory period.  As a result of these analyses, a simplified shaft seal model was developed for the CRA-2004 PA.

A conceptual representation of the simplified shaft seal system used in PA is shown in Figure PA-18.  The simplified model divides the shaft into three sections:  an upper section (shaft seal above the Salado), a lower section (within the Salado), and a concrete monolith section within the repository horizon.  A detailed discussion of how the material properties were assigned for the simplified shaft seal model is included in James and Stein (2003).  The permeability value used to represent the upper and lower sections is defined as the harmonic mean of the component materials’ permeability in the detailed shaft seal model (including permeability adjustments made for the DRZ assumed to surround the lower shaft seal section within the Salado).  Porosity is defined as the thickness-weighted mean porosity of the component materials.  Other material properties are described in James and Stein (2003).

The lower section of the shaft experiences a change in material properties at 200 years.  This change simulates the consolidation of seal materials within the Salado and significantly decreases permeability.  This time was chosen as a conservative overestimate of the amount of time expected for this section of the shaft to become consolidated.  The concrete monolith section of the shaft is unchanged from the CCA PA and is represented as being highly permeable for 10,000 years to ensure that fluids can access the north end (operations and experimental areas) in the model.  In three thin regions at the stratigraphic position of the anhydrite MBs, the shaft seal is modeled as MB material (Figure PA-18).  This model feature is included so that fluids flowing in the DRZ and MB fractures can access the interbeds to the north of the repository “around” the shaft seals.  Because these layers are so thin, they have virtually no effect on the effective permeability of the shaft seal itself.

The simplified shaft model was tested in the AP-106 analysis (Stein and Zelinski 2003), which supports the Salado Flow Peer Review (Caporuscio, Gibbons, and Oswald 2003).  The results of the AP-106 analysis demonstrate that vertical brine flow through the simplified shaft model is comparable to brine flows seen through the detailed shaft model used in the CCA PA and subsequent CCA PAVT calculations.

PA includes panel closures models that represent the Option D panel closure design (see the CRA-2004, Chapter 6.0, Section 6.4.3), which are unchanged from the CRA-2004.  Option D closures (Figure PA-19) are designed to allow minimal fluid flow between panels.  PA explicitly represents selected Option D panel closures in the computational grid using a model approved by the Salado Flow Peer Review Panel (Caporuscio, Gibbons, and Oswald 2003).  The Option D panel closure design has several components: an SMC monolith, which extends into the DRZ in all directions; an empty drift section; and a block and mortar explosion wall (Figure PA-19). Each set of panel closures are represented in the BRAGFLO grid by 4 materials in 13 grid cells (Figure PA-20):

Figure PA-18.    Schematic View of the Simplified Shaft Model (numbers on right indicate dimensions in meters)

Figure PA-19.  Schematic Side View of Option D Panel Closure

·       Six cells of panel closure concrete (area CONC_PCS, material CONC_PCS)

·       One cell above and one cell below the concrete material consisting of MB anhydrite (areas MB 139 and Anhydrite AB, materials S_MB139 and S_ANH_AB, respectively)

·       Two cells of healed DRZ above Anhydrite AB and the panel closure system (PCS) (area DRZ_PCS, material DRZ_PCS)

·       Three cells of empty drift and explosion wall (area DRF_PCS, material DRF_PCS)

Figure PA-20.  Representation of Option D Panel Closures in the BRAGFLO Grid

Properties for the materials making up the panel closure system are listed in Table PA-3.

The Option D panel closure design requires the use of a salt-saturated concrete, identified as SMC, as specified for the shaft seal system.  The design of the shaft seal system and the properties of SMC are described in the CCA, Appendix SEAL.  The BRAGFLO grid incorporates the material CONC_PCS, which is assigned the material properties of undegraded SMC and is used to represent the concrete portion of the Option D panel closure system (Figure PA-15).  A double-thickness concrete segment represents the northernmost set of panel closures (between the north RoR and the operations area).  This feature represents the two sets of panel closures in series that will be emplaced between the waste-filled repository and the shaft.

In the BRAGFLO grid, regions where the Option D panel closures intersect the MBs are represented as blocks of MB material (Figure PA-15).  This representation is warranted for two reasons:

1.   The MB material has a very similar permeability distribution (10 -21 to 10 -17.1 m2) to the concrete portion of the Option D panel closures (10 -20.699  to 10 -17 m2), and thus, assigning this material as anhydrite MB in the model has essentially the same effect as calling it concrete, as long as pressures are below the fracture initiation pressure.

2.   In the case of high pressures, it is expected that fracturing may occur in the anhydrite MBs and flow could go “around” the panel closures and out of the two-dimensional plane considered in the model grid.  In this case, the flow would be through the MB material, which incorporates a fracture model, as described above.

After constructing the concrete portion of the panel closure, the salt surrounding the monolith will be subjected to compressive stresses, which will facilitate the rapid healing of disturbed halite.  The rounded configuration of the monolith creates a situation very favorable for concrete durability:  high compressive stresses and low stress differences.  In turn, the compressive stresses developed within the salt will quickly heal any damage caused by construction excavation, thereby eliminating the DRZ along this portion of the panel closure.  The permeability of the salt immediately above and below the rigid concrete monolith component of Option D will approach the intrinsic permeability of the undisturbed Salado halite.

To represent the DRZ above the monoliths, PA uses the material DRZ_PCS in the BRAGFLO grid (Figure PA-15).  The values assigned to DRZ_PCS are the same as those used for the DRZ above the excavated areas (material DRZ_1, see Table PA-3), except for the properties PRMX_LOG, PRMY_LOG, and PRMZ_LOG, the logarithms of permeability in the x, y, and z directions, respectively.  These permeability values are assigned the same distributions used for the material CONC_PCS.  In this instance, the values are based on the nature of the model setup, not directly on experimental data (although the general range of the distribution agrees with experimental observations of healed salt).  The use of these permeabilities ensures that any fluid flow is equally probable through or around the Option D panel closures, and represents the range of uncertainty that exists in the performance of the panel closure system.

The DRF_PCS is the material representing the empty drift and explosion wall.  For simplicity, this material is assumed to have hydrologic properties equivalent to the material representing the waste panel and is used for the three sets of panel closures represented in the grid (Figure PA-15).  The creep closure model is applied to this material to be consistent with the neighboring materials.  The assignment of a high permeability to this region for the PA calculations, which contains the explosion wall, is justified because the explosion wall is not designed to withstand the stresses imposed by creep closure beyond the operational period and will be highly permeable following rapid room closure.

The major disruptive event in PA is the penetration of the repository by a drilling intrusion.  In the undisturbed scenario (see Section PA-6.7.1), these blocks have the material properties of the neighboring stratigraphic or excavated modeling unit, and there is no designation in the borehole grid except for the reduced lateral dimensions of this particular column of grid blocks.

In the scenarios simulating drilling disturbance, these cells start out with the same material properties as in the undisturbed scenario, but at the time of intrusion, the borehole grid blocks are reassigned to borehole material properties.  The drilling intrusion is modeled by modifying the permeability of the grid blocks in Column 26 of Figure PA-15 (values listed in Table PA-7).  Furthermore, the drilling intrusion is assumed to produce a borehole with a diameter of 12.25 in. (0.31 m) (Vaughn 1996, Howard 1996), borehole fill is assumed to be incompressible, capillary effects are ignored, residual gas and brine saturations are set to zero, and porosity is set to 0.32 (see materials CONC_PLG, BH_OPEN, BH_SAND, and BH_CREEP in Table PA-3).  When a borehole that penetrates pressurized brine in the Castile is simulated (i.e., an E1 intrusion), the permeability modifications indicated in Table PA-7 extend from the ground surface (i.e., Grid Cell 2155 in Figure PA-17) to the base of the pressurized brine (i.e., Grid Cell 2225 in Figure PA-17).  When a borehole that does not penetrate pressurized brine in the Castile is under consideration (i.e., an E2 intrusion), the permeability modifications indicated in Table PA-7 stop at the bottom of the lower DRZ (i.e., Grid Cell 1111 in Figure PA-17).

High-pressure Castile brine was encountered in several WIPP-area boreholes, including the WIPP-12 borehole within the controlled area and the ERDA-6 borehole northeast of the site.  Consequently, the conceptual model for the Castile includes the possibility that brine reservoirs underlie the repository.  The E1 and E1E2 scenarios include borehole penetration of both the repository and a brine reservoir in the Castile.

Unless a borehole penetrates both the repository and a brine reservoir in the Castile, the Castile is conceptually unimportant to PA because of its expected low permeability.  Two regions are specified in the disposal system geometry of the Castile horizon:  the Castile (Rows 1 and 2 in

Table PA-7.  Permeabilities for Drilling Intrusions Through the Repository

Time After Intrusion

Assigned Permeabilities

0–200 years

Concrete plugs are assumed to be emplaced at the Santa Rosa (i.e., a surface plug with a length of 15.76 m; corresponds to Grid Cells 2113, 2155 in Figure PA-17) and the Los Medaños Member of the Rustler (i.e., a plug at the top of the Salado with a length of 36 m; corresponds to Grid Cell 1644 in Figure PA-17).  Concrete plugs are assumed to have a permeability log-uniformly sampled between 10-19 m2 to 10-17m2 (see material CONC_PLG in Table PA-4).  The open portions of the borehole are assumed to have a permeability of 1 ´ 10-9 m2.

200–1200 years

Concrete plugs are assumed to fail after 200 years (U.S. Department of Energy 1995).  An entire borehole is assigned a permeability typical of silty sand log-uniformly sampled between 10-16.3 m2 and 10-11 m2 (see parameter BHPRM in Table PA-19 and material BH_SAND in Table PA-4).

> 1200 years

Permeability of borehole reduced by one order of magnitude in the Salado beneath the repository due to creep closure of borehole (Thompson et al. 1996) (i.e., k = 10x/10, x = BHPRM, in Grid Cells 2225, 1576, 26, 94, 162, 230 of Figure PA-17).  No changes are made within and above the lower DRZ (see material BH_CREEP in Table PA-4).

 

Figure PA-15) and a reservoir (Row 1, Columns 23 to 45 in Figure PA-15).  The Castile region has an extremely low permeability, which prevents it from participating in fluid flow processes.

It is unknown whether a brine reservoir exists below the repository.  As a result, the conceptual model for the brine reservoirs is somewhat different from those for known major properties of the natural barrier system, such as stratigraphy.  The principal difference is that a reasonable treatment of the uncertainty of the existence of a brine reservoir requires assumptions about the spatial distribution of such reservoir and the probability of intersection (see Attachment MASS-2009, Section MASS.18.1).  A range of probabilities for a borehole hitting a brine reservoir is used (see Section PA-3.6).

In addition to the stochastic uncertainty in the location and hence in the probability of intersecting reservoirs, there is also uncertainty in the properties of reservoirs.  The manner in which brine reservoirs would behave if penetrated is captured by parameter ranges and is incorporated in the BRAGFLO calculations of disposal system performance.  The conceptual model for the behavior of such a brine reservoir is discussed below.  The properties specified for brine reservoirs are pressure, permeability, compressibility, and porosity, and are sampled from parameter ranges (see Table PA-19).

Where they exist, Castile brine reservoirs in the northern Delaware Basin are believed to be fractured systems, with high-angle fractures spaced widely enough that a borehole can penetrate through a volume of rock containing a brine reservoir without intersecting any fractures, and therefore not producing brine.  Castile brine reservoirs occur in the upper portion of the Castile (Popielak et al. 1983).  Appreciable volumes of brine have been produced from several reservoirs in the Delaware Basin, but there is little direct information on the areal extent of the reservoirs or the existence of the interconnection between them.  Data from WIPP-12 and ERDA-6 indicate that fractures have a variety of apertures and permeabilities, and they deplete at different rates.  Brine occurrences in the Castile behave as reservoirs; that is, they are bounded systems.

Determining gas and brine flow in the vicinity of the repository requires solving the two nonlinear PDEs in Equation (PA.24) through Equation (PA.30) on the computational domain in Figure PA-15, along with evaluating appropriate auxiliary conditions.  The actual unknown functions in this solution are Pb and Sg, although the constraint conditions also give rise to values for Pg and Sb.  As two dimensions in space and one dimension in time are in use, Pb, Pg, Sb, and Sg are functions of the form Pb( x, y, t), Pg( x, y, t), Sb( x, y, t), and Sg( x, y, t).

Solving Equation (PA.24) through Equation (PA.30) requires both initial value and boundary value conditions for Pb and Sg.  The initial value conditions for Pb and Sg are given in Section PA-4.2.2.  As indicated there, the calculation starts at time t = −5 years, with a possible resetting of values at t = 0 years, which corresponds to final waste emplacement and sealing of the repository.  The boundary conditions are such that no brine or gas moves across the exterior grid boundary (Table PA-8).  This Neumann-type boundary condition is maintained for all time.  Further, BRAGFLO allows the user to maintain a specified pressure and/or saturation at any grid

Table PA-8.  Boundary Value Conditions for Pg and Pb

Boundaries below (Row 1, y = 0 m) and above (Row 33, y = 1039 m) system for 0 £ x £ 46630 m (Columns 1-68) and -5 yr £ t.  Below, j refers to the unit normal vector in the positive y direction.

×j = 0 Pa / m

No gas flow condition

×j = 0 Pa / m

No brine flow condition

Boundaries at left (Column 1, x = 0 m) and right (Column 68, x = 46630 m) of system for 0 £ y £ 1039 m (Rows 1-33) and -5 yr £ t.  Below, i refers to the unit normal vector in the positive x direction.

 ×i = 0 Pa / m

No gas flow condition

×i = 0 Pa / m

No brine flow condition

 

block.  This is not a boundary condition and is not required to close the problem.  This feature is used to specify Dirichlet-type conditions at the surface grid blocks (Columns 1-68, Row 33, Figure PA-15) and at the far-field locations in the Culebra and Magenta (Columns 1 and 68, Row 26, and Columns 1 and 68, Row 28, Figure PA-15).  These auxiliary conditions are summarized in Table PA-9.

A fully implicit finite-difference procedure is used to solve Equation (PA.24) through Equation (PA.30).  The associated discretization of the gas mass balance equation is given by


Table PA-9.  Auxiliary Dirichlet Conditions for Sg and Pb

Surface Grid Blocks

 = 0.08363

Columns  1–42, 44–68, Row 33, -5 yr £ t

Saturation is not forced at the shaft cell on the surface because its saturation is reset to 1.0 at t = 0 yr.

 = 1.01 ´ 105 Pa

Columns  1–68, row 33, –5 yr £ t

Culebra and Magenta Far Field

 = 9.14 ´ 105 Pa

i = 1 and 68, j = 26, –5 yr £ t (Culebra)

 = 9.47 ´ 105 Pa

i = 1 and 68, j = 28, –5 yr £ t (Magenta)

 

                                            (PA.91)

where F represents the phase potentials given by

and

the subscripts are defined by

              i  =  x-direction grid index

              j  =  y-direction grid index

      =  x-direction grid block interface

     =  y-direction grid block interface

             xi  =  grid block center in the x-coordinate direction (m)

             yj   =  grid block center in the y-coordinate direction (m)

           =  grid block length in the x-coordinate direction (m)

          =  grid block length in the y-coordinate direction (m)

the superscripts are defined by

             n  =  index in the time discretization, known solution time level

         n+1  =  index in the time discretization, unknown solution time level

and the interblock densities are defined by

           

           

           

           

The interface values of krg in Equation (PA.91) are evaluated using upstream weighted values (i.e., the relative permeabilities at each grid block interface are defined to be the relative permeabilities at the center of the adjacent grid block with the highest potential).  Further, interface values for arg kx/ mg and arg ky/ mg are obtained by harmonic averaging of adjacent grid block values for these expressions.  Currently all materials are isotropic, i.e. kx = ky = kz.

The discretization of the brine mass balance equation is obtained by replacing the subscript for gas, g, by the subscript for brine, b.  As a reminder, Pg and Sb are replaced in the numerical implementation with the substitutions indicated by Equation (PA.27) and Equation (PA.26), respectively.  Wells are not used in the conceptual model for long-term Salado flow calculations, but they are used for DBR calculations.  Thus, for long-term Salado flow calculations, the terms qg and qb are zero.  For long-term Salado flow calculations, the wellbore is not treated by a well model, but rather is explicitly modeled within the grid as a distinct material region (i.e., Upper Borehole and Lower Borehole in Figure PA-15).

The resultant coupled system of nonlinear brine and gas mass balance equations is integrated in time using the Newton-Raphson method with upstream weighting of the relative permeabilities, as previously indicated.  The primary unknowns at each computational cell center are brine pressure and gas saturation.

The Darcy velocity vectors vg( x, y, t) and vb( x, y, t) for gas and brine flow (m3/m2/s = m/s) are defined by the expressions

                                                                           (PA.92)

and

                                                                             (PA.93)

Values for vg and vb are obtained and saved as the numerical solution of Equation (PA.24) through Equation (PA.30) is carried out.  Cumulative flows of gas, Cg( t, B), and brine, Cb( t, B), from time 0 to time t across an arbitrary boundary B in the domain of (Figure PA-15) is then given by

                                                             (PA.94)

for l = g, b, where a(x, y) is the geometry factor defined in Equation (PA.32), n(x, y) is an outward-pointing unit normal vector, and  denotes a line integral.  As an example, B could correspond to the boundary of the waste disposal regions in Figure PA-15.  The integrals defining Cg( t, B) and Cb( t, B) are evaluated using the Darcy velocities defined by Equation (PA.92) and Equation (PA.93).  Due to the dependence of gas volume on pressure, Cg( t, B) is typically calculated in moles or cubic meters at standard temperature and pressure, which requires an appropriate change of units for vg in Equation (PA.95).

Additional information on BRAGFLO and its use in the CRA-2009 PA can be found in the BRAGFLO user’s manual (Nemer 2007c), the BRAGFLO design document (Nemer 2007b) and the analysis package for the Salado flow calculations in the CRA-2009 PA (Nemer and Clayton 2008).

The NUTS code is used to model radionuclide transport in the Salado.  NUTS models radionuclide transport within all regions for which BRAGFLO computes brine and gas flow, and for each realization uses as input the corresponding BRAGFLO velocity field, pressures, porosities, saturations, and other model parameters, including, for example, the geometrical grid, residual saturation, material map, and compressibility.  Of the radionuclides that are transported vertically due to an intrusion or up shaft, it is assumed that the lateral radionuclide transport is in the most transmissive unit, the Culebra.  Therefore, the radionuclide transport through the Dewey Lake to the accessible environment (fDL in Equation [PA.33]) and to the land surface due to long-term flow (fS in Equation [PA.33]) are set to zero for consistency.

The PA uses NUTS in two different modes.  First, the code is used in a computationally fast screening mode to identify those BRAGFLO realizations for which it is unnecessary to do full transport calculations because the amount of contaminated brine that reaches the Culebra or the LWB within the Salado is insufficient to significantly contribute to the total integrated release of radionuclides from the disposal system.  For the remaining realizations, which have the possibility of consequential release, a more computationally intensive calculation of each radionuclide’s full transport is performed (see Section PA-6.7.2).

This section describes the model used to compute radionuclide transport in the Salado for E0, E1, and E2 scenarios (defined in Section PA-2.3.2).  The model for transport in the E1E2 scenario, which is computed using the PANEL code, is described in Section PA-4.4.

NUTS models radionuclide transport by advection (see Attachment MASS-2009, Section MASS-13.5).  NUTS disregards sorptive and other retarding effects throughout the entire flow region.  Physically, some degree of retardation must occur at locations within the repository and the geologic media; it is therefore conservative to ignore retardation processes.  NUTS also ignores reaction-rate aspects of dissolution and colloid formation processes, and mobilization is assumed to occur instantaneously.  Neither molecular nor mechanical dispersion is modeled in NUTS.  These processes are assumed to be insignificant compared to advection, as discussed further in Attachment MASS-2009, Section MASS-13.5.

Colloidal actinides are subject to retardation by chemical interaction between colloids and solid surfaces and by clogging of small pore throats (i.e., by sieving).  There will be some interaction of colloids with solid surfaces in the anhydrite interbeds.  Given the low permeability of intact interbeds, it is likely that pore apertures will be small and some sieving will occur.  However, colloidal particles, if not retarded, are transported slightly more rapidly than the average velocity of the bulk liquid flow.  Because the effects on transport of slightly increased average pore velocity and retarded interactions with solid surfaces and sieving offset one another, the DOE assumes residual effects of these opposing processes will be either small or beneficial, and does not incorporate them when modeling An transport in the Salado interbeds.

If brine in the repository moves into interbeds, it is likely that mineral precipitation reactions will occur.  Precipitated minerals may contain actinides as trace constituents.  Furthermore, colloidal-sized precipitates will behave like mineral-fragment colloids, which are destabilized by brines, quickly agglomerating and settling by gravity.  The beneficial effects of precipitation and coprecipitation are neglected in PA.

Fractures, channeling, and viscous fingering may also impact transport in Salado interbeds, which contain natural fractures.  Because of the low permeability of unfractured anhydrite, most fluid flow in interbeds will occur in fractures.  Even though some properties of naturally fractured interbeds are characterized by in situ tests, other uncertainty exists in the characteristics of the fracture network that may be created with high gas pressure in the repository.  The PA modeling system accounts for the possible effects on porosity and permeability of fracturing by using a fracturing model (see Section PA-4.2.4).  The processes and effects associated with fracture dilation or fracture propagation not already captured by the PA fracture model are negligible (see the CCA, Appendix MASS, Section MASS.13.3 and Appendix MASS Attachment 13.2).  Of those processes not already incorporated, channeling has the greatest potential effect.

Channeling is the movement of fluid through the larger-aperture sections of a fracture network with locally high permeabilities.  It could locally enhance An transport.  However, it is assumed that the effects of channeled flow in existing or altered fractures will be negligible for the length and time scales associated with the disposal system.  The DOE believes this assumption is reasonable because processes are likely to occur that limit the effectiveness of channels or the dispersion of actinides in them.  First, if gas is present in the fracture network, it will be present as a nonwetting phase and will occupy the portions of the fracture network with relatively large apertures, where the highest local permeabilities will exist.  The presence of gas thus removes the most rapid transport pathways from the contaminated brine and decreases the impact of channeling.  Second, brine penetrating the Salado from the repository is likely to be completely miscible with in situ brine.  Because of miscibility, diffusion or other local mixing processes will probably broaden fingers (reduce concentration gradients) until the propagating fingers are indistinguishable from the advancing front.

Gas will likely penetrate the liquid-saturated interbeds as a fingered front, rather than a uniform front.  Fingers form when there is a difference in viscosity between the invading fluid (gas) and the resident fluid (liquid brine), and because of channeling effects.  This process does not affect An transport, however, because actinides of interest are transported only in the liquid phase, which will not displace gas in the relatively high-permeability regions due to capillary effects.

The following system of PDEs is used to model radionuclide transport in the Salado:

               -Ñ×       (PA.95)

                                                                            (PA.96)

for l = 1, 2, …, nR, where

      =  Darcy velocity vector (m3/m2/s = m/s) for brine (supplied by BRAGFLO from solution of Equation (PA.93))

     Cbl  =  concentration (kg/m3) of radionuclide l in brine

     Csl  =  concentration (kg/m3) of radionuclide l in solid phase (i.e., not in brine), with concentration defined with respect to total (i.e., bulk) formation volume (only used in repository; see Figure PA-15)

      Sl  =  linkage term (kg/m3/s) due to dissolution/precipitation between radionuclide l in brine and in solid phase (see Equation (PA.97))

       f  =  porosity (supplied by BRAGFLO from solution of Equation (PA.24) through Equation (PA.30))

      Sb  =  brine saturation (supplied by BRAGFLO from solution of Equation (PA.24) through Equation (PA.30))

      ll  =  decay constant (s -1) for radionuclide l

   P(l)  =  {p:  radionuclide p is a parent of radionuclide l}

      nR  =  number of radionuclides,

and a is the dimension-dependent geometry factor in Equation (PA.32).  PA uses a two-dimensional representation for fluid flow and radionuclide transport in the vicinity of the repository, with a defined by the element depths in Figure PA-15.  Although omitted for brevity, the terms a, vb, Cbl, Csl, Sl, Sb, and f are functions a(x, y), vb( x, y, t), Cbl( x, y, t), Csl( x, y, t), Sl( x, y, t), Sb( x, y, t), and f(x, y, t) of time t and the spatial variables x and y.  Equation (PA.95) and Equation (PA.96) are defined and solved on the same computational grid (Figure PA-15) used by BRAGFLO for the solution of Equation (PA.24) through Equation (PA.30).

Radionuclides are assumed to be present in both brine (Equation (PA.95)) and in an immobile solid phase (Equation (PA.96)), although radionuclide transport takes place only by brine flow (Equation (PA.95)).  The maximum radionuclide concentration is assumed to equilibrate instantly for each element (Section PA-4.3.2).  Then each individual radionuclide equilibrates between the brine and solid phases based on the maximum concentration of its associated element and the mole fractions of other isotopes included in the calculation.  The linkage between the brine and solid phases in Equation (PA.95) and Equation (PA.96) is accomplished by the term Sl, where

         
                                                                                                                                         (
PA.97)

where

         =  maximum concentration (kg/m3) of element El(l) in oxidation state Ox(l) in brine type Br(t), where El(l) denotes the element of which radionuclide l is an isotope, Ox(l) denotes the oxidation state in which element El(l) is present, and Br(t) denotes the type of brine present in the repository at time t (see Section PA-4.3.2 for definition of ).

                                                       =  concentration (kg/m3) of element El(l) in brine (p = b) or solid (p = s), which is equal to the sum of concentrations of radionuclides that are isotopes of same element as radionuclide l, where k Î El(l) only if k is an isotope of element El(l):

                                                                                                   (PA.98)

                            =  difference (kg/m3) between maximum concentration of element El(l) in brine and existing concentration of element El(l) in brine

                                                (PA.99)

                                     MFpl = mole fraction of radionuclide l in phase p, where p = b (brine) or
p = s (solid)

                                                                              (PA.100)

                                      CMl = conversion factor (mol/kg) from kilograms to moles for radionuclide l

                                                      = Dirac delta function (s -1) (d(t - t) = 0 if t ¹ t and )

In the CCA and the CRA-2004, the function ST had an additional parameter, Mi, representing the presence or absence of microbial activity. However, during preparations for the CRA-2004 PABC, the EPA directed the DOE to change the probability of microbial activity to 1.0 (Cotsworth 2005). The DOE therefore revised the probability distribution so that all vectors would have the possibility for microbial degradation of cellulose, while retaining the percentage of vectors showing microbial degradation of rubbers and plastic at 0.25 (Nemer 2005). This, in effect, sets Mi equal to “Yes” for all vectors.  Mi is therefore suppressed as an argument of ST in the remainder of this section.

The terms Sl, Sb, Cp,El(l), MFpl , and f are functions of time t and the spatial variables x and y, although the dependencies are omitted for brevity.  The Dirac delta function, d(tt), appears in Equation (PA.97) to indicate that the adjustments to concentration are implemented instantaneously within the numerical solution of Equation (PA.95) and Equation (PA.96) whenever a concentration imbalance is observed.

The velocity vectorvb in Equation (PA.95) and Equation (PA.96) is defined in Equation (PA.93) and is obtained from the numerical solution of Equation (PA.24) through Equation (PA.30).  If B denotes an arbitrary boundary (e.g., the LWB) in the domain of Equation (PA.95) and Equation (PA.96) (as shown in Figure PA-15), the cumulative transport of Cl( t, B) of radionuclide l from time 0 to time t across B is given by

                                            (PA.101)

where n(x, y) is an outward-pointing unit normal vector and  denotes a line integral over B.

Equation (PA.95) and Equation (PA.96) models advective radionuclide transport due to the velocity vector vb.  Although the effects of solubility limits are considered, no chemical or physical retardation is included in the model.  Molecular diffusion is not included in the model, because the radionuclides under consideration have molecular diffusion coefficients on the order of 10 -10 m2/s, and thus can be expected to move only approximately 10 m over 10,000 years due to molecular diffusion.  Mechanical dispersion is also omitted in NUTS, as the uniform initial radionuclide concentrations assumed within the repository and the use of time-integrated releases in assessing compliance with section 191.13 mean that mechanical dispersion will have negligible impact on overall performance.

A maximum concentration ST( Br, Ox, El) (mol/liter [L]) is calculated for each brine type (Br Î {Salado, Castile}), oxidation state (Ox Î {III, IV, V, VI}), and element (El Î {Am, Pu, U, Th}).  The maximum concentration is given by

                                                          (PA.102)

where SD( Br, Ox) is the dissolved solubility (mol/L) and SC( Br, Ox, El) is the concentration (mol/L) of the element sorbed to colloids.

The dissolved solubility SD( Br, Ox) is given by

                                                                     (PA.103)

where

    SFMT( Br, Ox) = dissolved solubility (mol/L) calculated by Fracture-Matrix Transport (FMT) model (WIPP Performance Assessment 1998a) for brine type Br and oxidation state Ox

                  = logarithm (base 10) of uncertainty factor for solubilities calculated by FMT expressed as a function of oxidation state Ox

Table PA-10 lists the calculated values of SFMT( Br, Ox); details of the calculation are provided in Appendix SOTERM-2009 and Brush (2005).  The uncertainty factors UF(Ox) are determined by the uncertain parameters WSOLVAR3 for the III oxidation state and WSOLVAR4 for the IV oxidation state; definitions of these parameters are provided in Table PA-19, and further details regarding the calculation of the solubilities is given in Appendix SOTERM-2009, Section SOTERM-5.1.3.

Table PA-10.     Calculated Values for Dissolved Solubility (see Appendix SOTERM-2009, Table SOTERM-16)

Brine

Oxidation State

III

IV

V

VIa

Salado

3.87 ´ 10-7

5.64 ´ 10-8

3.55 ´ 10-6

1.0 ´ 10-3

Castile

2.88 ´ 10-7

6.79 ´ 10-8

8.24 ´ 10-7

1.0 ´ 10-3

a  Values for the VI solubilities were mandated by the EPA (U.S. Environmental Protection Agency 2005).

 

The concentration (mol/L) of the element sorbed to colloids SC (Br, Ox, El) is given by

                                                                                          (PA.104)

where

             = solubility (concentration expressed in mol/L) in brine type Br of element El in oxidation state Ox resulting from humic colloid formation

                                      = 

         = scale factor used as a multiplier on SD( Br, Ox) in definition of SHum( Br, Ox, El) (see Table PA-11)

                          UBHum  = upper bound on solubility (concentration expressed in mol/L) of individual An elements resulting from humic colloid formation

                                      =  1.1 ´ 10-5 mol/L

               = solubility (concentration expressed in mol/L) in brine type Br of element El in oxidation state Ox resulting from microbial colloid formation

                                      =

                     = scale factor used as multiplier on SD( Br, Ox) in definition of SMic( Br, Ox, El) (see Table PA-12)

                    = upper bound on solubility (concentration expressed in mol/L) of element El in oxidation state Ox resulting from microbial colloid formation (see Table PA-12)

Table PA-11.     Scale Factor SFHum( Br, Ox, El) Used in Definition of SHum( Br, Ox, El) (see Appendix SOTERM-2009, Table SOTERM-21)

Brine

Element (Oxidation State)

Am(III)

Pu(III)

Pu(IV)

U(IV)

Th(IV)

Np(IV)

Np(V)

U(VI)

Salado

0.19

0.19

6.3

6.3

6.3

6.3

9.1 ´ 10-4

0.12

Castile

WPHUMOX3a

WPHUMOX3a

6.3

6.3

6.3

6.3

7.4 ´ 10-3

0.51

a  See Table PA-19.

 

Table PA-12.     Scale Factor SFMic( Ox, El) and Upper Bound UBMic( Ox, El) (mol/L) Used in Definition of SMic( Br, Ox, El) (see Appendix SOTERM-2009, Table SOTERM-21)

Parameter

Element (Oxidation State)

Am(III)

Pu(III)

Pu(IV)

U(IV)

Th(IV)

Np(IV)

Np(V)

U(VI)

SFMic(Ox, El)

3.6

0.3

0.3

0.0021

3.1

12.0

12.0

0.0021

UBMic(Ox, El)

1

6.8 ´ 10-5

6.8 ´ 10-5

0.0021

0.0019

0.0027

0.0027

0.0021

 

                                   = solubility (concentration in mol/L) of element El resulting from An intrinsic colloid formation

                                      = 

                               SMn  = solubility (concentration in mol/L) of individual An element resulting from mineral fragment colloid formation

                                      =  2.6 ´ 10 -8 mol/L

Since the solution of Equation (PA.95) and Equation (PA.96) for this many radionuclides and decay chains would be very time-consuming, the number of radionuclides for direct inclusion in the analysis was reduced using the algorithm shown in Figure PA-21.  The CRA-2009 PA uses the same reduction algorithm as the CCA PA (see the CCA, Appendix WCA); the algorithm was found to be acceptable in the CCA review (U.S. Environmental Protection Agency 1998a, Section 4.6.1.1).

Figure PA-21. Selecting Radionuclides for the Release Pathways Conceptualized by PA


Using Figure PA-21, the number of radionuclides considered was reduced from 135 to a total of 29 included in the decay calculations carried out by the PANEL code (Garner and Leigh 2005). These radionuclides belong to the following decay chains:

                                                        (PA.105)

                                                                            (PA.106)

                            (PA.107)

                                           (PA.108)

Radionuclides considered in the decay calculations that do not belong to one of the decay chains listed above are 147Pm, 137Cs, and 90Sr.  In addition, some intermediates with extremely short half-lives, such as 240U, were omitted from the decay chains.

Further simplification of the decay chains is possible based on the total inventories.  Releases of radionuclides whose inventories total less than one EPA unit are essentially insignificant, as any release that transports essentially all of a given species outside the LWB will be dominated by the releases of other species with much larger inventories.  In addition, 137Cs and 90Sr can be omitted because their concentrations drop to below 1 EPA unit within 150 years, which makes it improbable that a significant release of these radionuclides will occur.  Isotopes such as 241Pu, whose decay could affect the inventory of measurable isotopes, were reinstated.

After the reduction of radionuclides outlined inFigure PA-21 and the above paragraph, the following 10 radionuclides remained from the decay chains shown above:

                                                                              (PA.109)

                                                                                                                            (PA.110)

                                                                                                                            (PA.111)

                                                                            (PA.112)

238Pu does not significantly affect transport calculations because of its short half-life (87.8 years).  The remaining nine radionuclides were then further reduced by combining those with similar decay and transport properties.  In particular, 234U, 230Th, and 239Pu were used as surrogates for the groups {234U, 233U}, {230Th, 229Th}, and {242Pu, 240Pu, 239Pu}, with the initial inventories of 234U, 230Th, and 239Pu being increased to account for the additional radionuclide(s) in each group.

In increasing the initial inventories, the individual radionuclides were combined on either a mole or curie basis (i.e., moles added and then converted back to curies, or curies added directly).  In each case, the method that maximized the combined inventory was used; thus, 233U was added to 234U, 240Pu to 239Pu, and 229Th to 230Th by curies, while 242Pu was added to 239Pu by moles. In addition, 241Pu was added to 241Am by moles because 241Pu has a half-life of 14 years and will quickly decay to 241Am, and neglect of this ingrowth would underestimate the 241Am inventory by about 4% (Table PA-13).  The outcome of this process was the following set of five radionuclides in three simplified decay chains:

                                                                (PA.113)

which were then used with Equation (PA.95) and Equation (PA.96) for transport in the vicinity of the repository.  As 238Pu does not significantly affect transport calculations, only the remaining four radionuclides are used in the Culebra transport calculations (Section PA-4.9).  These radionuclides account for 99% of the EPA units in the waste after 2,000 years (Garner and Leigh 2005), and hence will dominate any releases by transport.

Table PA-13.  Combination of Radionuclides for Transport

Combination

Isotope Initial Values

Combination Procedure

Combined Inventory

233U ® 234U

1.23 ´ 103 Ci 233U
3.44 ´ 102 Ci 234U

1.23 ´ 103 Ci 233U
    ® 1.23 ´ 103 Ci 234U

1.57 ´ 103 Ci 234U

242Pu ® 239Pu

 

 

 

 

240Pu ® 239Pu

1.27 ´ 101 Ci 242Pu




9.55 ´ 104 Ci 240Pu
5.82 ´ 105 Ci 239Pu

1.27 ´ 101 Ci 242Pu
    = 1.32 ´ 101 moles 242Pu
    ® 1.32 ´ 101 moles 239Pu
    = 1.98 ´ 102 Ci 239Pu

9.55 ´ 104 Ci 240Pu
     ® 9.55 ´ 104 Ci 239Pu

6.78 ´ 105 Ci 239Pu

229Th ® 230Th

5.21 ´ 100 Ci 229Th
1.80 ´ 10-1 Ci 230Th

5.21 ´ 100 Ci 229Th
    ® 5.21 ´ 100 Ci 230Th

5.39 ´ 100 Ci 230Th

241Pu ® 241Am

4.48 ´ 105 Ci 241Pu

5.18 ´ 105 Ci 241Am

5.38 ´ 105 Ci 241Pu
    = 1.79 ´ 101 moles 241Pu
    ® 1.79 ´ 101 moles 241Am
    = 1.50 ´ 104 Ci 241Am

5.33 ´ 105 Ci 241Am

 

All BRAGFLO realizations are first evaluated using NUTS in a screening mode to identify those realizations for which a significant release of radionuclides to the LWB cannot occur.  The screening simulations consider an infinitely soluble, nondecaying, nondispersive, and nonsorbing species as a tracer element.  The tracer is given a unit concentration in all waste disposal areas of 1 kg/m3.  If the amount of tracer that reaches the selected boundaries (the top of the Salado and the LWB within the Salado) does not exceed a cumulative mass of 10–7 kg within 10,000 years, it is assumed there is no consequential release to these boundaries.  If the cumulative mass outside the boundaries within 10,000 years exceeds 10–7 kg, a complete transport analysis is conducted.  The value of 10 -7 kg is selected because, regardless of the isotopic composition of the release, it corresponds to a normalized release less than 10–6 EPA units, the smallest release displayed in CCDF construction (Stockman 1996).  The largest normalized release would be 9.98 × 10–7 EPA units, corresponding to 10 -7 kg of 241Am if the release was entirely 241Am.

For BRAGFLO realizations with greater than 10 -7 kg reaching the boundaries in the tracer calculations, NUTS models the transport of five different radionuclide species (241Am, 239Pu, 238Pu, 234U, and 230Th).  These radionuclides represent a larger number of radionuclides:  as discussed in Section PA-4.3.3, radionuclides were grouped together based on similarities, such as isotopes of the same element and those with similar half-lives, to simplify the calculations.  For transport purposes, solubilities are lumped to represent both dissolved and colloidal forms.  These groupings simplify and expedite calculations.

Equation (PA.95) and Equation (PA.96) are numerically solved by the NUTS program (WIPP Performance Assessment 1997a) on the same computational grid (Figure PA-15) used by BRAGFLO for the solution of Equation (PA.24), Equation (PA.25), Equation (PA.26), Equation (PA.27), Equation (PA.28), Equation (PA.29), and Equation (PA.30).  In the solution procedure, Equation (PA.95) and Equation (PA.96) are numerically solved with Sl = 0 for each time step, with the instantaneous updating of concentrations indicated in Equation (PA.97) and the appropriate modification to Csl in Equation (PA.96) taking place after the time step.  The solution is carried out for the five radionuclides indicated in Equation (PA.113).

The initial value and boundary value conditions used with Equation (PA.95) and Equation (PA.96) are given in Table PA-14.  At time t = 0 (corresponding to the year 2033), the total inventory of each radionuclide is assumed to be in brine; the solubility constraints associated with Equation (PA.97) then immediately adjust the values for Cbl( x, y, t) and Csl( x, y, t) for consistency with the constraints imposed by ST( Br, Ox, El) and available radionuclide inventory.

The nR PDEs in Equation (PA.95) and Equation (PA.96) are discretized in two dimensions and then developed into a linear system of algebraic equations for numerical implementation.  The following conventions are used in the representation of each discretized equation:

·       The subscript b is dropped from Cbl, so that the unknown function is represented by Cl.

·       A superscript n denotes time tn, with the assumption that the solution Cl is known at time tn and is to be propagated to time tn +1.

·       The grid indices are i in the x-direction and j in the y-direction, and are the same as the BRAGFLO grid indices.

Table PA-14.  Initial and Boundary Conditions for Cbl( x, y, t) and Csl( x, y, t)

Initial Conditions for Cbl(x, y, t) and Csl(x, y, t)

   =   if (x, y) is a point in the repository (i.e., areas Waste Panel, South RoR and North RoR, in Figure PA-15), where Al(0) is the amount (kg) of radionuclide l present at time t = 0 (Table PA-13) and Vb(0) is the amount (m3) of brine in repository at time t = 0 (from solution of Equation (PA.24) through Equation (PA.30) with BRAGFLO) for all (x, y).

                       = 0  otherwise.

= 0  if (x, y) is a point in the repository.

Boundary Conditions for Cbl(x, y, t)

            = , where B is any subset of the outer boundary of the computational grid in Figure PA-15,  is the flux (kg/s) at time t of radionuclide l across B, vb(x, y, t) is the Darcy velocity (m3/m2/s) of brine at (x, y) on B and is obtained from the solution of Equation (PA.24) through Equation (PA.30) by BRAGFLO, n(x, y) denotes an outward-pointing unit normal vector, and  denotes a line integral along B.

 

·       Fractional indices refer to quantities evaluated at grid block interfaces.

·       Each time step by NUTS is equal to 20 BRAGFLO time steps because BRAGFLO stores results (here, vb, f, and Sb) every 20 time steps.

The following finite-difference discretization is used for the lth equation in each grid block (i, j):

                                                                                                                      (PA.114)

where qb is the grid block interfacial brine flow rate (m3/s) and VR is the grid block volume (m3).  The quantity qb is based on  and a in Equation (PA.95) and Equation (PA.96), and the quantity VR is based on grid block dimensions (Figure PA-15) and a.

The interfacial values of concentration in Equation (PA.114) are discretized using the one-point upstream weighting method (Aziz and Settari 1979), which results in

                                                                       (PA.115)

where w derives from the upstream weighting for flow between adjacent grid blocks and is defined by

By collecting similar terms, Equation (PA.115) can be represented by the linear equation

                                     (PA.116)

where

Given the form of Equation (PA.116), the solution of Equation (PA.95) and Equation (PA.96) has now been reduced to the solution of nR ´ nG linear algebraic equations in nR ´ nG unknowns, where nR is the number of equations for each grid block (i.e., the number of radionuclides) and nG is the number of grid blocks into which the spatial domain is discretized (Figure PA-15).

The system of PDEs in Equation (PA.95) and Equation (PA.96) is strongly coupled because of the contribution from parental decay to the equation governing the immediate daughter.  Consequently, a sequential method is used to solve for the radionuclide concentrations by starting at the top of a decay chain and working down from parent to daughter.  This implies that when solving Equation (PA.116) for the lth isotope concentration, all parent concentrations occurring in the right-hand-side term R are known.  The system of equations is then linear in the concentrations of the lth isotope.  As a result, solving Equation (PA.95) and Equation (PA.96) is reduced from the solution of one algebraic equation at each time step with nR ´ nG unknowns to the solution of nR algebraic equations each with nG unknowns at each time step, which can result in a significant computational savings.

The matrix resulting from one-point upstream weighting has the following structural form for a 3 ´ 3 system of grid blocks, and a similar structure for a larger number of grid blocks:

 

1

2

3

4

5

6

7

8

9

1

X

X

0

X

 

 

 

 

 

2

X

X

X

0

X

 

 

 

 

3

0

X

X

0

0

X

 

 

 

4

X

0

0

X

X

0

X

 

 

5

 

X

0

X

X

X

0

X

 

6

 

 

X

0

X

X

0

0

X

7

 

 

 

X

0

0

X

X

0

8

 

 

 

 

X

0

X

X

X

9

 

 

 

 

 

X

0

X

X

 

where X designates possible nonzero matrix entries, and 0 designates zero entries within the banded structure. All entries outside of the banded structure are zero.  Because of this structure, a banded direct elimination solver (Aziz and Settari 1979, Section 8.2.1) is used to solve the linear system for each radionuclide.  The bandwidth is minimized by first indexing equations in the coordinate direction with the minimum number of grid blocks.  The coefficient matrix is stored in this banded structure, and all infill coefficients calculated during the elimination procedure are contained within the band structure.  Therefore, for the matrix system in two dimensions, a pentadiagonal matrix of dimension IBW ´ nG is inverted instead of a full nG ´ nG matrix, where IBW is the bandwidth.

The numerical implementation of Equation (PA.96) enters the solution process through updates to the radionuclide concentrations in Equation (PA.115) between each time step, as indicated in Equation (PA.97).  The numerical solution of Equation (PA.95) and Equation (PA.96) also generates the concentrations required to integrate numerically evaluating the integral that defines Cl( t, B) in Equation (PA.101).

Additional information on NUTS and its use in WIPP PA can be found in the NUTS users manual (WIPP Performance Assessment 1997a) and in the analysis package of Salado transport calculations for the CRA-2009 PA (Ismail and Garner 2008).  Furthermore, additional information on dissolved and colloidal actinides is given in Appendix SOTERM-2009.

This section describes the model used to compute radionuclide transport in the Salado for the E1E2 scenario.  The model for transport in E0, E1, and E2 scenarios is described in Section PA-4.3.

A relatively simple mixed-cell model is used for radionuclide transport in the vicinity of the repository after an E1E2 intrusion, when connecting flow between two drilling intrusions into the same waste panel is assumed to take place.  With this model, the amount of radionuclide l contained in a waste panel is represented by

                                                                              (PA.117)

where

       =  amount (mol) of radionuclide l in waste panel at time t

     =  concentration (mol/m3) of radionuclide l in brine in waste panel at time t (Equation (PA.118) and Equation (PA.119))

        =  rate (m3/s) at which brine flows out of the repository at time t (supplied by BRAGFLO from solution of Equation (PA.93))

and ll and P(l) are defined in conjunction with Equation (PA.95) and Equation (PA.96).

The brine concentration Cbl in Equation (PA.117) is defined by

                           (PA.118)

                               (PA.119)

where

     = mole fraction of radionuclide l in waste panel at time t

                                                                 =                                             (PA.120)

          = volume (m3) of brine in waste panel at time t (supplied by BRAGFLO from solution of Equation (PA.24), Equation (PA.25), Equation (PA.26), Equation (PA.27), Equation (PA.28), Equation (PA.29), and Equation (PA.30))

and ST[ Br, Ox, El] is supplied by Equation (PA.102).

For Equation (PA.118) and Equation (PA.119), ST[ Br, Ox, El] must be expressed in units of mol/L.  In other words, Cbl( t) is defined to be the maximum concentration (ST in Equation (PA.102)) if there is sufficient radionuclide inventory in the waste panel to generate this concentration (Equation (PA.118)); otherwise, Cbl( t) is defined by the concentration that results when all the relevant element in the waste panel is placed in solution (Equation (PA.119)).  The dissolved and colloidal An equilibrate instantly for each element.

Given rb and Cbl, evaluation of the integral

                                                                                            (PA.121)

provides the cumulative release Rl( t) of radionuclide l from the waste panel through time t.

Equation (PA.117) is numerically evaluated by the PANEL model (WIPP Performance Assessment 1998b) using a discretization based on time steps of 50 years or less.  Specifically, Equation (PA.117) is evaluated with the approximation

              
                                                                                                                                       (
PA.122)

where

      = gain in radionuclide l due to the decay of precursor radionuclides between tn

and tn+1 (see Equation (PA.123)),  = .

As the solution progresses, values for Cbl( tn) are updated in consistency with Equation (PA.118) and Equation (PA.119), and the products rb( tn) Cbl( tn) are accumulated to provide an approximation to Rl in Equation (PA.121).

The term Gl( tn, tn+1) in Equation (PA.122) is evaluated with the Bateman equations (Bateman 1910), with PANEL programmed to handle decay chains of up to five (four decay daughters for a given radionuclide).  As a single example, if radionuclide l is the third radionuclide in a decay chain (i.e., l = 3) and the two preceding radionuclides in the decay chain are designated by l = 1 and l = 2, then

                                                                       (PA.123)

in Equation (PA.122).

The preceding model is used in two ways in PA.  First, Equation (PA.121) estimates releases to the Culebra associated with E1E2 intrusion scenarios (see Section PA-6.7.3).  Second, the radionuclide concentrations are computed using the minimum brine volume for a significant release to estimate DBRs (see Section PA-6.8.2.3).  The calculation of the minimum brine volume used in the CRA-2009 PA is found in Stein (2005). The calculated concentrations are the Sl term indicated in Equation (PA.97) which are used in the NUTS calculations discussed in Section PA-4.3.

For E1E2 intrusions, the initial amount Al of radionuclide l is the inventory of the decayed isotope at the time of the E1 intrusion.  PANEL calculates the inventory of each of the 29 radioisotopes throughout the regulatory period.  The initial concentration Cbl of radionuclide l is computed by Equation (PA.117), Equation (PA.118), and Equation (PA.119).  For the DBR calculations, the initial amount Al of radionuclide l is the inventory of the isotope at the time of repository closure.

Additional information on PANEL and its use in the CRA-2009 PA calculations can be found in the PANEL user’s manual (WIPP Performance Assessment 2003b), the analysis package for PANEL calculations (Garner and Leigh 2005), and the analysis package for Salado transport calculations in the CRA-2009 PA (Ismail and Garner 2008).

Cuttings are waste solids contained in the cylindrical volume created by the cutting action of the drill bit passing through the waste, while cavings are additional waste solids eroded from the borehole by the upward-flowing drilling fluid within the borehole.  The releases associated with these processes are computed within the CUTTINGS_S code (WIPP Performance Assessment 2003c).  The mathematical representations used for cuttings and cavings are described in this section.

The uncompacted volume of cuttings removed and transported to the surface in the drilling fluid, Vcut, is given by

                                                                                                (PA.124)

where A is the drill bit area (m2), Hi is the initial (or uncompacted) repository height (3.96 m), and D is the drill-bit diameter (0.31115 m) (Fox 2008, Table 13).  For drilling intrusions through RH-TRU waste, Hi = 0.509 m is used (Fox 2008, Table 45).

The cavings component of the direct surface release is caused by the shearing action of the drilling fluid on the waste as it flows up the borehole annulus.  Like the cuttings release, the cavings release is assumed to be independent of the conditions that exist in the repository during a drilling intrusion.

The final diameter of the borehole depends on the diameter of the drillbit and on the extent to which the actual borehole diameter exceeds the drill-bit diameter.  Although a number of factors affect erosion within a borehole (Chambre Syndicale de la Recherche et de la Production du Petrole et du Gaz Naturel 1982), the most important is the fluid shear stress on the borehole wall (i.e., the shearing force per unit area, N/m2) resulting from circulating drilling fluids (Darley 1969, Walker and Holman 1971).  As a result, PA estimates cavings removal with a model based on the effect of shear stress on the borehole diameter.  In particular, the borehole diameter is assumed to grow until the shear stress on the borehole wall is equal to the shear strength of the waste, which is the limit below which waste erosion ceases.

The final eroded diameter Df (m) of the borehole through the waste determines the total volume V (m3) of uncompacted waste removed to the surface by circulating drilling fluid.  Specifically,

                                                                                         (PA.125)

where Vcav is the volume (m3) of waste removed as cavings.

Most borehole erosion is believed to occur in the vicinity of the drill collar (Figure PA-22) because of decreased flow area and consequent increased mud velocity (Rechard, Iuzzolino, and Sandha 1990, Letters 1a and 1b, App. A).  An important determinant of the extent of this erosion is whether the flow of the drilling fluid in the vicinity of the collar is laminar or turbulent.  PA uses Reynolds numbers to distinguish between the occurrence of laminar flow and turbulent flow.  The Reynolds number is the ratio between inertial and viscous (or shear) forces in a fluid, and can be expressed as (Fox and McDonald 1985)

                                                                                                                 (PA.126)

where Re is the Reynolds number (dimensionless), rf is the fluid density (kg/m3), De is the equivalent diameter (m),  is the fluid speed (m s -1), and h is the fluid viscosity (kg m -1 s -1).

Typically, rf, v, and h are averages over a control volume with an equivalent diameter of De, where rf = 1.21 ´ 103 kg/m3 (Fox 2008, Table 13), v = 0.7089 m s -1 (based on 40 gal/min/in of drill diameter) (Berglund 1992), and De = 2 (R - Ri), as shown in Figure PA-22.  The diameter of the drill collar (i.e., 2Ri in Figure PA-22) is 8.0 in = 0.2032 m (Ismail 2008).  The determination of h is discussed below.  PA assumes that Reynolds numbers less than 2100 are associated with laminar flow, while Reynolds numbers greater than 2100 are associated with turbulent flow (Walker 1976).

Drilling fluids are modeled as non-Newtonian, which means that the viscosity h is a function of the shear rate within the fluid (i.e., the rate at which the fluid velocity changes normal to the flow direction, m/s/m).  PA uses a model proposed by Oldroyd (1958) to estimate the viscosity of drilling fluids.  As discussed in the Drilling Mud and Cement Slurry Rheology Manual (Chambre Syndicale de la Recherche et de la Production du Petrole et du Gaz Naturel 1982), the Oldroyd model leads to the following expression for the Reynolds number associated with the helical flow of a drilling fluid within an annulus:

                                                                                                      (PA.127)

where rf, De, and v are defined as in Equation (PA.126), and h ¥ is the asymptotic value for the derivative of the shear stress (t , kg m -1 s -2) with respect to the shear rate (G, s -1) obtained as the shear rate increases (i.e.,  as ).  PA uses Equation (PA.127) to determine whether drilling fluids in the area of the drill collar are undergoing laminar or turbulent flow.

The Oldroyd model assumes that the shear stress t  is related to the shear rate G through the relationship

                                                                                                      (PA.128)

where h0 is the asymptotic value of the viscosity (kg m -1 s -1) that results as the shear rate G approaches zero, and s1 and s2 are constants (s2).  The expression leads to

                                                                                                               (PA.129)

PA uses values of h0 = 1.834 ´ 10 -2 kg m -1 s -1, s1 = 1.082 ´ 10 -6 s2, and s2 = 5.410 ´ 10 -7 s2 (Berglund 1996), from which viscosity in the limit of infinite shear rate is found to be h ¥ =

Figure PA-22.  Detail of Rotary Drill String Adjacent to Drill Bit

9.17 ´ 10 -3 kg m -1 s -1.  The quantity h ¥ is comparable to the plastic viscosity of the fluid (Chambre Syndicale de la Recherche et de la Production du Petrole et du Gaz Naturel 1982).

As previously indicated, different models are used to determine the eroded diameter Df of a borehole depending on whether flow in the vicinity of the drill collar is laminar or turbulent.  The model for borehole erosion in the presence of laminar flow is described next, and then the model for borehole erosion in the presence of turbulent flow is described.

As shown by Savins and Wallick (1966), the shear stresses associated with the laminar helical flow of a non-Newtonian fluid, as a function of the normalized radius, r, can be expressed as

                                                                             (PA.130)

for Ri/ R £ r £ 1, where Ri and R are the inner and outer radii within which the flow occurs, as indicated in Figure PA-22; t(R,r) is the shear stress (kg m -1 s -2) at a radial distance DR beyond the inner boundary (i.e., at r = (Ri + DR)/R); and the variables C,J, and l depend on R and satisfy conditions Equation (PA.132) through Equation (PA.134) indicated below.  The shear stress at the outer radius R is given by

                                                                                      (PA.131)

As previously indicated, the borehole radius R is assumed to increase as a result of erosional processes until a value of R is reached at which t(R, 1) is equal to the shear strength of the waste.  In PA, the shear strength of the waste is treated as an uncertain parameter (see WTAUFAIL in Table PA-19).  Computationally, determining the eroded borehole diameter R associated with a particular value of the waste shear strength requires repeated evaluation of t (R, 1), as indicated in Equation (PA.131), until a value of R is determined for which t(R, 1) equals the shear strength.

The quantities C, J, and l must satisfy the following three conditions (Savins and Wallick 1966) for Equation (PA.131) to be valid:

                                                                                               (PA.132)

                                                                                                       (PA.133)

                                                       (PA.134)

where h, the drilling fluid viscosity (kg m -1 s -1), is a function of R and r; DW is the drill string angular velocity (rad s -1); and Q is the drilling fluid flow rate (m3 s -1).

The viscosity h in Equation (PA.132) through Equation (PA.134) is introduced into the analysis by assuming that the drilling fluid follows the Oldroyd model for shear stress in Equation (PA.128).  By definition of the viscosity h,

                                                                         t = hG                                                   (PA.135)

and from Equation (PA.128)

                                                                                                          (PA.136)

thus the expression in Equation (PA.130) can be reformulated as

                                                                   (PA.137)

As discussed by Savins and Wallick (1966) and Berglund (1992), the expressions in Equation (PA.132) through Equation (PA.134) and Equation (PA.136) can be numerically evaluated to obtain C, J, and l for use in Equation (PA.130) and Equation (PA.131).  In PA, the drill string angular velocity DW is treated as an uncertain parameter (see DOMEGA in Table PA-19), and

                                                                                                          (PA.138)

where v = 0.7089 m s -1 as used in Equation (PA.126), and h0, s 1, and s 2 are defined as in Equation (PA.128) and Equation (PA.129).

The model for borehole erosion in the presence of turbulent flow is now described.  Unlike the theoretically derived relationship for erosion in the presence of laminar flow, the model for borehole erosion in the presence of turbulent flow is empirical.  In particular, pressure loss for axial flow in an annulus under turbulent flow conditions can be approximated by (Chambre Syndicale de la Recherche et de la Production du Petrole et du Gaz Naturel 1982)

                                                                                                            (PA.139)

where DP is the pressure change (Pa), f is the Fanning friction factor (dimensionless), L is the distance (m) over which pressure change DP occurs, and rf , v, and De are defined in Equation (PA.126).

For turbulent pipe flow, f is empirically related to the Reynolds number Re defined in Equation (PA.126) by (Whittaker 1985)

                                                                                (PA.140)

where D is the inside diameter (m) of the pipe and e is a “roughness term” equal to the average depth (m) of pipe wall irregularities.  In the absence of a similar equation for flow in an annulus, Equation (PA.140) is used in PA to define f for use in Equation (PA.139), with D replaced by the effective diameter De = 2(RRi) and e equal to the average depth of irregularities in the waste-borehole interface.  In the present analysis, e  = 0.025 m (Fox 2008, Table 34), which exceeds the value often selected in calculations involving very rough concrete or riveted steel piping (Streeter 1958).

The pressure change DP in Equation (PA.139) and the corresponding shear stress t  at the walls of the annulus are approximately related by

                                                                               (PA.141)

where  is the cross-sectional area of the annulus (see Figure PA-22) and 2pL(R + Ri) is the total surface area of the annulus.  Rearranging Equation (PA.139) and using the relationship in Equation (PA.135) yields

                                                                                                          (PA.142)

which was used in the 1991 and 1992 WIPP PAs to define the shear stress at the surface of a borehole of radius R.  The radius R enters into Equation (PA.132) through Equation (PA.134) through the use of D = 2(RRi) in the definition of f in Equation (PA.140).  As with laminar flow, the borehole radius R is assumed to increase until a value of t(R) is reached that equals the sample value for the shear strength of the waste (i.e., the uncertain parameter WTAUFAIL in Table PA-19).  Computationally, the eroded borehole diameter is determined by solving Equation (PA.142) for R under the assumption that t(R) equals the assumed shear strength of the waste.

For the CRA-2004 PA, a slight modification to the definition of t  in Equation (PA.142) was made to account for drill string rotation when fluid flow in the vicinity of the drill collars is turbulent (Abdul Khader and Rao 1974; Bilgen, Boulos, and Akgungor 1973).  Specifically, an axial flow velocity correction factor (i.e., a rotation factor), Fr, was introduced into the definition of t.  The correction factor Fr is defined by

                                                                     Fr = v2100 / v                                               (PA.143)

where v2100 is the norm of the flow velocity required for the eroded diameters to be the same for turbulent and laminar flow at a Reynolds number of Re = 2100, and is obtained by solving

                                                                                                           (PA.144)

for v2100 with D in the definition of f in Equation (PA.140) assigned the final diameter value that results for laminar flow at a Reynolds number of Re = 2100 (that is, the D in De = 2(RRi) = D – 2Ri obtained from Equation (PA.127) with Re = 2100).  The modified definition of t  is

                                                                                                         (PA.145)

and results in turbulent and laminar flow with the same eroded diameter at a Reynolds number of 2100, where PA assumes that the transition between turbulent and laminar flow takes place.

The following algorithm was used to determine the final eroded radius Rf of a borehole and incorporates a possible transition from turbulent to laminar fluid flow within a borehole:

Step 1.   Use Equation (PA.127) to determine an initial Reynolds number Re, with R initially set to the drill-bit radius, R0 = 0.31115 m (Fox 2008, Table 13).

Step 2.   If Re < 2100, the flow is laminar and the procedure in Section PA-4.5.2.1 is used to determine Rf.  Because any increase in the borehole diameter will cause the Reynolds number to decrease, the flow will remain laminar and there is no need to consider the possibility of turbulent flow as the borehole diameter increases, with the result that Rf  determined in this step is the final eroded radius of the borehole.

Step 3.   If Re ³ 2100, then the flow is turbulent, and the procedure discussed in Section PA-4.5.2.2 is used to determine Rf .  Once Rf is determined, the associated Reynolds number Re is recalculated using Equation (PA.127) and R = Rf .  If the recalculated Re > 2100, a transition from turbulent to laminar flow cannot take place, and the final eroded radius is Rf determined in this step.  If not, go to Step 4.

Step 4.   If  the Reynolds number Re with the new Rf in Step 3 satisfies the inequality Re £ 2100, a transition from turbulent to laminar flow is assumed to have taken place.  In this case, Rf is recalculated assuming laminar flow, with the outer borehole radius R initially defined to be the radius associated with Re = 2100. In particular, the initial value for R is given by the radius at which the transition from laminar to turbulent flow takes place:

                                                                                                   (PA.146)

              which is obtained from Equation (PA.127) by solving for R with Re = 2100.  A new value for Rf is then calculated with the procedure discussed in Section PA-4.5.2.1 for laminar flow, with this value of Rf replacing the value from Step 3 as the final eroded diameter of the borehole.

Step 5.   Once Rf is known, the amount of waste removed to the surface is determined using Equation (PA.125) with Df = 2Rf .

Additional information on CUTTINGS_S and its use in the CRA-2004 PA to determine cuttings and cavings releases can be found in the CUTTINGS_S user’s manual (WIPP Performance Assessment 2003c) and in the analysis package for cuttings and cavings releases (Ismail 2008).

Spallings are waste solids introduced into a borehole by the movement of waste-generated gas towards the lower-pressure borehole.  In engineering literature, the term “spalling” describes the dynamic fracture of a solid material, such as rock or metal (Antoun et al. 2003).  In WIPP PA, the spallings model describes a series of processes, including tensile failure of solid waste, fluidization of failed material, entrainment into the wellbore flow, and transport up the wellbore to the land surface.  Spallings releases could occur when pressure differences between the repository and the wellbore cause solid stresses in the waste exceeding the waste material strength and gas velocities sufficient to mobilize failed waste material.

The spallings model is described in the following sections.  Presented first are the primary modeling assumptions used to build the conceptual model.  Next, the mathematical model and its numerical implementation in the computer code DRSPALL are described.  Finally, implementation of the spallings model in WIPP PA by means of the code CUTTINGS_S is discussed.

Assumptions underlying the spallings model include the future state of the waste, specifications of drilling equipment, and the driller’s actions at the time of intrusion.  Consistent with the other PA models, the spallings model assumes massive degradation of the emplaced waste through mechanical compaction, corrosion, and biodegradation.  Waste is modeled as a homogeneous, isotropic, weakly consolidated material with uniform particle size and shape.  The rationale for selecting the spallings model material properties is addressed in detail by Hansen et al. (1997) and Hansen, Pfeifle, and Lord (2003).

Drilling equipment specifications, such as bit diameter and drilling mud density, are based on surveys of drillers in the Delaware Basin (Hansen, Pfeifle, and Lord 2003).  Assumptions about the driller’s actions during the intrusion are conservative.  Typically, the drilling mud density is controlled to maintain a slightly “overbalanced” condition so that the mud pressure is always slightly higher than the fluid pressures in the formation.  If the borehole suddenly passes through a high-pressure zone, the well can quickly become “underbalanced,” with a resulting fluid pressure gradient driving formation fluids into the wellbore.  This situation is known as a kick and is of great concern to drillers because a violent kick can lead to a blowout of mud, gas, and oil from the wellbore, leading to equipment damage and worker injury.  Standard drilling practice is to watch diligently for kicks.  The first indicator of a kick is typically an increase in mud return rate, leading to an increase in mud pit volume (Frigaard and Humphries 1997).  Downhole monitors detect whether the kick is air, H2S, or brine.  If the kick fluid is air, the standard procedure is to stop drilling and continue pumping mud in order to circulate the air pocket out.  If the mud return rate continues to grow after drilling has stopped and the driller believes that the kick is sufficiently large to cause damage, the well may be shut in by closing the blowout preventer.  Once shut in, the well pressure may be bled off slowly and mud weight eventually increased and circulated to offset the higher formation pressure before drilling continues.  The spallings model simulates an underbalanced system in which a gas kick is assured, and the kick proceeds with no intervention from the drill operation.  Therefore, drilling and pumping continue during the entire blowout event.

The spallings model calculates transient repository and wellbore fluid flow before, during, and after a drilling intrusion.  To simplify the calculations, both the wellbore and the repository are modeled by one-dimensional geometries.  The wellbore assumes a compressible Newtonian fluid consisting of a mixture of mud, gas, salt, and waste solids; viscosity of the mixture varies with the fraction of waste solids in the flow.  In the repository, flow is viscous, isothermal, compressible single-phase (gas) flow in a porous medium.

The wellbore and repository flows are coupled by a cylinder of porous media before penetration, and by a cavity representing the bottom of the borehole after penetration.  Schematic diagrams of the flow geometry prior to and after penetration are shown in Figure PA-23 and Figure PA-24, respectively.  The drill bit moves downward as a function of time, removing salt or waste material.  After penetration, waste solids freed by drilling, tensile failure, and associated fluidization may enter the wellbore flow stream at the cavity forming the repository-wellbore boundary.

Figure PA-23.  Schematic Diagram of the Flow Geometry Prior to Repository Penetration

Figure PA-24.  Schematic Diagram of the Flow Geometry After Repository Penetration

Flow in the well is modeled as a one-dimensional pipe flow with cross-sectional areas corresponding to the appropriate flow area at a given position in the well, as shown in Figure PA-25 and Figure PA-26.  This model is conceptually similar to that proposed by Podio and Yang (1986) for use in the oil and gas industry.  Drilling mud is added at the wellbore entrance by the pump.  Flow through the drill bit is treated as a choke with cross-sectional area appropriate for the bit nozzle area.  At the annulus output to the surface, the mixture is ejected at a constant atmospheric pressure.  The gravitational body force acts in its appropriate direction based on position before or after the bit.

Figure PA-25.  Effective Wellbore Flow Geometry Before Bit Penetration

Figure PA-26.  Effective Wellbore Flow Geometry After Bit Penetration

Prior to drill bit penetration into the repository, gas from the repository can flow through drilling-damaged salt into the well.  After penetration, the cavity at the bottom of the wellbore couples the wellbore flow and the repository flow models; gas and waste material can exit the repository domain into the cavity.  The cavity radius increases as waste materials are moved into the wellbore.

The system of equations representing flow in the wellbore consists of four equations for mass conservation, one for each phase (salt, waste, mud, and gas); one equation for conservation of total momentum; two equations relating gas and mud density to pressure; the definition of density for the fluid mixture; and one constraint imposed by the fixed volume of the wellbore.  The conservation of mass and momentum is described by

                                                                                        (PA.147)

                                                     (PA.148)

where

          q  =  phase (w for waste, s for salt, m for mud, and g for gas)

        Vq    = volume (m3) of phase q

          V  = total volume (m3)

         rq    = density (kg/m3) of phase q, constant for salt and waste (2,180 and 2,650 kg/m3, respectively) and pressure-dependent for gas and mud (see Equation (PA.149) and Equation (PA.150))

          r  =  density of fluid mixture (kg/m3) determined by Equation (PA.151)

          u  = velocity (m/s) of fluid mixture in wellbore

           t  = time (s)

          z  = distance (m) from inlet at top of well

         Sq    = rate of mass (kg/s) in phase q entering and exiting wellbore domain at position z (Equation (PA.162))

     Smom    = rate of momentum (kg m/s2) entering and exiting wellbore domain at position z (Equation (PA.165))

          P  = pressure (Pa) at position z

          g  = standard gravity (9.8067 kg/m/s2)

          F  = friction loss using pipe flow model (kg/m2/s2) determined by Equation (PA.153)

Gas is treated as isothermal and ideal, so the pressure and density are related by Boyle’s law:

                                                                                                                 (PA.149)

where rg,0  is the density of H2 gas at atmospheric pressure and 298 K (8.24182 ´ 10-2 kg/m3).

The mud is assumed to be a compressible fluid, so

                                                                                      (PA.150)

where rm,0 is the density of the mud at atmospheric pressure (1,210 kg/m3) and cm is the compressibility of the mud (3.1 ´ 10-10 Pa-1).

The density of the fluid mixture is determined from the densities and volumes occupied by the phases:

                                                                               (PA.151)

The volume of each phase is constrained by the fixed total volume of the wellbore:

                                                                                                   (PA.152)

The friction loss is a standard formulation for pipe flow (Fox and McDonald 1985), where the head loss per unit length is given as

                                                                                                           (PA.153)

where the hydraulic diameter dh is given by

                                                                                                          (PA.154)

In PA, Do = 0.31115 m throughout the domain.  From the bit to the top of the collar, Di = 0.2032 m; above the collar, Di = 0.1143 m.  The area A is calculated as the area of the annulus between the outer and inner radii:

                                                                                                            (PA.155)

Thus, dh = 0.108 m from the bit to the top of the collar, and dh = 0.197 m above the collar.

The Darcy friction factor f in Equation (PA.153) is determined by the method of Colebrook (Fox and MacDonald 1985).  In the laminar regime, which is assumed to be characterized by Reynolds numbers below 2100 (Walker 1976),

                                                                                                                          (PA.156)

and in the turbulent regime (Re > 2100)

                                                                           (PA.157)

where  is the Reynolds number of the mixture, and h is the viscosity calculated in Equation (PA.158), below.  As the wellbore mixture becomes particle-laden, the viscosity of the mixture is determined from an empirical relationship developed for proppant slurry flows in channels for the oil and gas industry (Barree and Conway 1995).  Viscosity is computed by an approximate slurry formula based on the volume fraction of waste solids:

                                                                                                       (PA.158)

where h0 is a base mixture viscosity (9.17 ´ 10 -3 Pa s), w = Vw/ V is the current volume fraction of waste solids, wmax is an empirically determined maximal volume fraction above which flow is choked (0.615), and s is an empirically determined constant (-1.5) (Hansen, Pfeifle, and Lord 2003).

Initial conditions in the wellbore approximate mixture flow conditions just prior to waste penetration.  The wellbore is assumed to contain only mud and salt.  Initial conditions for the pressure, fluid density, volume fractions of mud and salt, and the mixture velocity are set by the following algorithm:

Step 1.    Set pressure in the wellbore to hydrostatic: P(z) = Patm – rm,0 gz.

Step 2.    Set mud density using Equation (PA.150).

Step 3.    Set mixture velocity: u(z) = Rm/ A(z), where Rm is the volume flow rate of the pump (0.0202 m3/s), and A(z) is the cross-sectional area of the wellbore.

Step 4.    Set volume of salt in each cell: Vs,i= Rdrill AbitDzi/ ui, where Rdrill is the rate of drilling (0.004445 m/s),  is the area of the bottom of the wellbore, and dbit is the diameter of the bit (0.31115 m).

Step 5.    Set volume fraction of mud in each cell: Vm,i  = Vi – Vs,i.

Step 6.    Recalculate mixture density using Equation (PA.151), assuming no waste or gas in the wellbore.

The initial conditions set by this algorithm approximate a solution to the wellbore flow (Equation (PA.147) and Equation (PA.148)) for constant flow of mud and salt in the well.  The approximation rapidly converges to a solution for wellbore flow if steady-state conditions are maintained (WIPP Performance Assessment 2003d).

For simplicity, DRSPALL does not model flow of mud down the pipe to the bit.  Mass can enter the wellbore below the drill bit and exit at the wellbore outlet.  Below the bit, mud, salt, gas, and waste can enter the wellbore.  PA assumes a constant volume of mud flow down the drilling pipe; therefore, the source term for mud, Sm,in, is set by the volumetric flow rate of the pump Rm (0.0202 m3/s) and the density of the mud at the bottom of the wellbore:

                                                                                                               (PA.159)

Until the drill bit penetrates the repository, salt enters the wellbore at a constant rate:

                                                                                                       (PA.160)

Additional mass enters the wellbore by gas flow from the repository (Sgas,in) and spalling of waste material (Sw,in); these mass sources are discussed in Section PA-4.6.2.3.  The outlet of the wellbore is set to atmospheric pressure.  Mass exiting the wellbore is determined from the mixture velocity, the area of the outlet Aout (0.066 m2), and the density and volume fraction of each phase at the outlet of the wellbore:

                                                                                                  (PA.161)

Finally, the net change in mass and momentum for phase q is

                                                                                                        (PA.162)

                                                                                            (PA.163)

The outlet of the wellbore is set to atmospheric pressure.  Momentum exiting the wellbore is determined from the fluid velocity and the area of the outlet Aout (0.066 m2):

                                                                                                (PA.164)

No momentum is added by mass flow into the wellbore from the repository; thus

                                                                                          (PA.165)

The repository is modeled as a radially symmetric domain.  A spherical coordinate system is used for this presentation and for most DRSPALL calculations.  In a few circumstances, cylindrical coordinates are used in PA calculations, where spall volumes are large enough that spherical coordinates are not representative of the physical process (Lord, Rudeen, and Hansen 2003).  Cylindrical coordinates are also available; the design document for DRSPALL (WIPP Performance Assessment 2003e) provides details on implementing the repository flow model in cylindrical coordinates.

Flow in the repository is transient, compressible, viscous, and single-phase (gas) flow in a porous medium.  Gas is treated as isothermal and ideal.  The equations governing flow in the repository are the equation of state for ideal gases (written in the form of Boyle’s law for an ideal gas at constant temperature), conservation of mass, and Darcy’s law with the Forchheimer correction (Aronson 1986, Whitaker 1996):

                                                                                                                 (PA.166)

                                                                                                   (PA.167)

                                                                                                      (PA.168)

where

          P  = pressure in pore space (Pa)

         rg   =  density of gas (kg/m3)

          u  =  velocity of gas in pore space (m/s)

          f =  porosity of the solid (unitless)

        hg   = gas viscosity (8.934 ´ 10-6 Pa s)

          k = permeability of waste solid (m2)

          F = Forchheimer correction (unitless)

The Forchheimer correction is included in Equation (PA.168) to account for inertia in the flowing gas, which becomes important at high gas velocities (Ruth and Ma 1992).  When the Forchheimer coefficient is zero, Equation (PA.168) reduces to Darcy’s law.  A derivation of Equation (PA.168) from the Navier-Stokes equations is given by Whitaker (1996); the derivation suggests that F is a linear function of gas velocity for a wide range of Reynolds numbers.

In PA, the Forchheimer correction takes the form

                                                                      F = bnd ru                                                 (PA.169)

where bnd is the non-Darcy coefficient, which depends on material properties such as the tortuosity and area of internal flow channels, and is empirically determined (Belhaj et al. 2003).  DRSPALL uses a value from Li et al. (2001) that measured high-velocity nitrogen flow through porous sandstone wafers, giving the result

                                                                                                         (PA.170)

Equation (PA.166) through Equation (PA.168) combines into a single equation for pressure in the porous solid:

                                                                            (PA.171)

where

                                                                                                 (PA.172)

and the Laplacian operator in a radially symmetric coordinate system is given by

                                                                                               (PA.173)

where n = 2 and n = 3 for polar and spherical coordinates, respectively.

In DRSPALL, the permeability of the waste solid is a subjectively uncertain parameter that is constant for waste material that has not failed and fluidized.  In a region of waste that has failed, the permeability increases as the waste fluidizes by a factor of 1 + Ff, where Ff is the fraction of failed material that has fluidized and is based on the fluidization relaxation time. This approximately accounts for the material bulking as it fluidizes.

Initial pressure in the repository is set to a constant value Pff.  A no-flow boundary condition is imposed at the outer boundary (r = R):

                                                                      ÑP(R) = 0                                                (PA.174)

At the inner boundary (r = rcav), the pressure is specified as P(rcav, t) = Pcav( t), where Pcav( t) is defined in the next section.  The cavity radius rcav increases as drilling progresses and waste material fails and moves into the wellbore; calculation of rcav is described in Section PA-4.6.2.3.3.

Prior to penetration, a cylinder of altered-permeability salt material with diameter equal to the drill bit is assumed to connect the bottom of the wellbore to the repository.  At the junction of the repository and this cylinder of salt, a small, artificial cavity is used to determine the boundary pressure for repository flow.  After penetration, the cavity merges with the bottom of the wellbore to connect the wellbore to the repository.

The cylinder of salt connecting the wellbore to the repository is referred to as the DDZ in Figure PA-23.  The permeability of the DDZ, kDDZ, is 1 ´ 10-14 m2.  The spallings model starts with the bit 0.15 m above the repository; the bit advances at a rate of Rdrill = 0.004445 m/s.

To couple the repository to the DDZ, the model uses an artificial pseudo-cavity in the small hemispherical region of the repository below the wellbore with the same surface area as the bottom of the wellbore (Figure PA-26).  The pseudo-cavity is a numerical device that smoothes the discontinuities in pressure and flow that would otherwise occur upon bit penetration of the repository.  The pseudo-cavity contains only gas, and is initially at repository pressure.  The mass of gas in the cavity mcav is given by

                                                                                                    (PA.175)

where

          Srep  = gas flow from repository into pseudo-cavity (kg/s); see Equation (PA.176)

         Sg, in  = gas flow from pseudo-cavity through DDZ into wellbore (kg/s); see Equation (PA.177)

Flow from the repository into the pseudo-cavity is given by

                                                                                                (PA.176)

where

        rg,rep  = gas density in repository at cavity surface (kg/m3) =

          urep  = gas velocity (m/s) in repository at cavity surface =

             f = porosity of waste (unitless)

          Acav  = surface area of hemispherical part of the cavity (m2)

                 = , where dbit is the diameter of the bit (m)

Flow out of the pseudo-cavity through the DDZ and into the wellbore is modeled as steady-state using Darcy’s Law:

                                                                            (PA.177)

where

           hg  = viscosity of H2 gas (8.934 ´ 10 -6 Pa s)

          Mw  = molecular weight of H2 gas (0.00202 kg / mol)

             R = ideal gas constant (8.314 J/mol K)

             T = repository temperature (constant at 300 K (27 ºC; 80 ºF))

             L = length (m) of DDZ (from bottom of borehole to top of repository)

         Pcav  =  pressure in pseudo-cavity (Pa)

          PBH  =  pressure at bottom of wellbore (Pa)

A justification for using this steady-state equation is provided in the design document for DRSPALL (WIPP Performance Assessment 2003e).  The pseudo-cavity is initially filled with gas at a pressure of Pff.  The boundary pressure on the well side (PBH) is the pressure immediately below the bit, determined by Equation (PA.147) and Equation (PA.148).  The pressure in the pseudo-cavity (Pcav) is determined by the ideal gas law:

                                                                                                            (PA.178)

where mcav is the number of moles of gas in the cavity and the cavity volume Vcav is given by

                                                                                                                 (PA.179)

In PA, the drilling rate into the ground is assumed constant at 0.004445 m/s; thus L = Li – 0.004445t until L = 0, at which time the bit penetrates the waste. The term Li is the distance from the bit to the waste at the start of calculation (0.15 m).

After waste penetration, the bottom of the wellbore is modeled as a hemispherical cavity in the repository, the radius of which grows as drilling progresses and as material fails and moves into the cavity.  Gas, drilling mud, and waste are assumed to thoroughly mix in this cavity; the resulting mixture flows around the drill collars and then up the annulus between the wellbore and the drill string.  Gas flow from the repository into the cavity is given by Equation (PA.176); however, Acav is now dependent on the increasing radius of the cavity (see Section PA-4.6.2.3.3).  Mudflow into the cavity from the wellbore is given by Equation (PA.159).  Waste flow into the cavity is possible if the waste fails and fluidizes; these mechanisms are discussed in Section PA-4.6.2.3.4 and Section PA-4.6.2.3.5.  Pressure in the cavity is equal to that at the bottom of the wellbore, and is computed by Equation (PA.178).

The cylindrical cavity of increasing depth created by drilling is mapped to a hemispherical volume at the bottom of the wellbore to form the cavity.  This mapping maintains equal surface areas in order to preserve the gas flux from the repository to the wellbore.  The cavity radius from drilling is thus

                                                                                              (PA.180)

where DH is the depth of the drilled cylinder.  In PA, the drilling rate into the ground is assumed constant at 0.004445 m/s; thus DH = 0.004445t until DH = H, the height of compacted waste (m).  Since the initial height of the repository is 3.96 m, H is computed from the porosity f by , where f0 is the initial porosity of a waste-filled room.

The cavity radius rcav is increased by the radius of failed and fluidized material rfluid, which is the depth to which fluidization has occurred beyond the drilled radius.  That is,

                                                                                                        (PA.181)

Gas flow from the waste creates a pressure gradient within the waste, which induces elastic stresses in addition to the far-field confining stress.  These stresses may lead to tensile failure of the waste material, an assumed prerequisite to spallings releases.  While the fluid calculations using Equation (PA.166), Equation (PA.167), and Equation (PA.168) are fully transient, the elastic stress calculations are assumed to be quasi-static (i.e., sound-speed phenomena in the solid are ignored).  Elastic effective stresses are (Timoshenko and Goodier 1970)

                                 (PA.182)

                             (PA.183)

where b is Biot’s constant (assumed here to be 1.0) and sff is the confining far-field stress (assumed constant at 14.8 MPa).

The flow-related radial and tangential stresses (ssr and ss q , respectively) are computed by equations analogous to differential thermal expansion (Timoshenko and Goodier 1970):

                                                               (PA.184)

                               (PA.185)

where Pff is the initial repository pressure and u is Poisson’s ratio (0.38).

Since stresses are calculated as quasi-static, an initial stress reduction caused by an instantaneous pressure drop at the cavity face propagates instantaneously through the waste.  The result of calculating Equation (PA.182) can be an instantaneous early-time tensile failure of the entire repository if the boundary pressure is allowed to change suddenly.  This is nonphysical and merely a result of the quasi-static stress assumption, combined with the true transient pore pressure and flow-related stress equations.  To prevent this nonphysical behavior, tensile failure propagation is limited by a tensile failure velocity (1000 m/s; see Hansen et al. 1997).  This limit has no quantitative effect on results, other than to prevent nonphysical tensile failure.

At the cavity face, Equation (PA.182) and Equation (PA.184) evaluate to zero, consistent with the quasi-static stress assumption.  This implies that the waste immediately at the cavity face cannot experience tensile failure; however, tensile failure may occur at some distance into the waste material.  Consequently, the radial effective stress sr is averaged from the cavity boundary into the waste over a characteristic length Lt (0.02 m).  If this average radial stress  is tensile and its magnitude exceeds the material tensile strength (|| > TENSLSTR), the waste is no longer capable of supporting radial stress and fails, permitting fluidization.  The waste tensile strength is an uncertain parameter in the analysis (see TENSLSTR in Table PA-15).

Equation (PA.183) and Equation (PA.185) evaluate shear stresses in the waste.  DRSPALL does not use the waste shear stresses to calculate waste failure for spall releases.  These stresses are included in this discussion for completeness.

Failed waste material is assumed to be disaggregated, but not in motion; it remains as a porous, bedded material lining the cavity face, and is treated as a continuous part of the repository from the perspective of the porous flow calculations.  The bedded material may be mobilized and enter the wellbore if the gas velocity in the failed material (see Equation (PA.168)) exceeds a minimum fluidization velocity, Uf.  The minimum fluidization velocity is determined by solving the following quadratic equation (Cherimisinoff and Cherimisinoff 1984, Ergun 1952)

                            (PA.186)

where

       a = particle shape factor (unitless)

      dp  = particle diameter (m)

Fluidization occurs in the failed material to the depth at which gas velocity does not exceed the fluidization velocity; this depth is denoted by rfluid and is used to determine cavity radius (Section PA-4.6.2.3.3).  If fluidization occurs, the gas and waste particles mix into the cavity at the bottom of the wellbore.  Because this mixing cannot be instantaneous, which would be nonphysical (much as allowing instantaneous tensile failure propagation would be nonphysical), a small artificial relaxation time, equal to the cavity radius rcav divided by the superficial gas velocity u(rcav), is imposed upon the mixing phenomenon.  The fluidized material is released into the cavity uniformly over the relaxation time.

The numerical model implements the conceptual and mathematical models described above (Section PA-4.6.2).  Both the wellbore and the repository domain calculations use time-marching finite differences.  These are part of a single computational loop and therefore use the same time step.  The differencing schemes for the wellbore and repository calculations are similar, but not identical.

The wellbore is zoned for finite differencing, as illustrated in Figure PA-27, which shows zones, zone indices, grid boundaries, volumes, and interface areas.  The method is Eulerian:  zone boundaries are fixed, and fluid flows across the interfaces by advection.  Quantities are zone-centered and integration is explicit in time.

To reduce computation time, an iterative scheme is employed to update the wellbore flow solution.  The finite-difference scheme first solves Equation (PA.147) and Equation (PA.148) for the mass of each phase in each grid cell and the momentum in each grid cell.

The updated solution to Equation (PA.147) and Equation (PA.148) is then used to compute the volume of each phase, the pressure, and the mixture velocity in each grid cell.

All of the materials (mud, salt, gas, and waste) are assumed to move together as a mixture.  Because fluid moves through the cell boundaries, the calculation requires a value for the flow through each cell boundary during a time step.  These values are obtained by averaging the fluid velocities at the zone centers, given by

Figure PA-27.  Finite-Difference Zoning for Wellbore

                                                                                               (PA.187)

The mass transport equation, prior to any volume change, becomes

                      (PA.188)

Here, the source terms Sm,i correspond to material entering or exiting at the pump, cavity, and surface. The “upwind” zone-centered densities are used for the interfaces values,  and .

Finally, any changed volumes are incorporated and numerical mass diffusion is added for stability:

                                                                          (PA.189)

where

and zq is the diffusion coefficient for phase q.  The density rfq for phase q being diffused is calculated from the mixture density, r, and the mass fraction, fq , of phase q in the referenced cell (fq = rVq,i/ rVi). The numerical diffusion coefficient zq is chosen empirically for stability.  Separate diffusion coefficients could be used for the different materials (mud, gas, etc.); however, sufficient stability is obtained by diffusing only mud and salt using the same coefficient (zm = zs = 0.0001 and zw = zg = 0).

Momentum is differenced as