Title 40 CFR Part 191
Subparts B and C
Compliance Recertification Application 2014
for the
Waste Isolation Pilot Plant
Appendix TFIELD2014
Transmissivity Fields
United States Department of Energy
Waste Isolation Pilot Plant
Carlsbad Field Office
Carlsbad, New Mexico
Compliance Recertification Application 2014
Appendix TFIELD2014
TFIELD1.0 Overview of the Tfield Development, Calibration, and Mining Modification Process
TFIELD2.1 Culebra Hydrogeologic Setting
TFIELD2.2 Refinement of Geologic Boundaries
TFIELD2.2.1 Rustler Halite Margins
TFIELD2.2.2 Salado Dissolution
TFIELD2.3 Confinement and Recharge in the Culebra
TFIELD2.3.1 Surface Drainage Basins
TFIELD2.3.2 Culebra Confinement
TFIELD2.3.3 Gypsum Cements in the Culebra
TFIELD3.0 TField Conceptual Model Refinement
TFIELD3.2 Overburden Thickness
TFIELD3.3 Fracture Interconnection
TFIELD3.5 Rustler Halite Margins
TFIELD3.6 Transmissivity Regression Model
TFIELD3.7 Culebra Conceptual Model Peer Review
TFIELD4.0 Base TField Construction
TFIELD4.1 Step 1  Linear Regression Analysis
TFIELD4.2 Step 2  Creation of Soft Data
TFIELD4.2.3 Diffusivity and Hydraulic Connections
TFIELD4.2.4 Combined Soft Data
TFIELD4.3 Step 3  Indicator Variography
TFIELD4.4 Step 4  Conditional Indicator Simulation
TFIELD4.5 Step 5  Construction of Transmissivity Fields
TFIELD5.0 TField Calibration
TFIELD5.1 Model Calibration Targets
TFIELD5.2 Step 1  Calibration Setup and Configuration
TFIELD5.2.1 Creation of Parameter Zones
TFIELD5.2.1.1 Transmissivity Zones
TFIELD5.2.1.2 Horizontal Anisotropy Zones
TFIELD5.2.1.3 Storativity Zones
TFIELD5.2.1.5 Flow, NoFlow, and FixedHead Zones
TFIELD5.2.2 Selection of Pilot Point Locations
TFIELD5.2.3 TransmissivitySpecific Pilot Point Settings
TFIELD5.2.4 AnisotropySpecific Pilot Point Settings
TFIELD5.2.5 StorativitySpecific Pilot Point Settings
TFIELD5.2.6 RechargeSpecific Pilot Point Settings
TFIELD5.2.7 Selection of Initial Values
TFIELD5.2.7.1 Parameter Initial Values
TFIELD5.2.7.2 Initial Head Field
TFIELD5.2.8 Creation of Transmissivity Fields
TFIELD5.2.9 Observations and Residuals
TFIELD5.2.10 MODFLOW Numerical Model
TFIELD5.3 Step 2 Calibration Process
TFIELD5.3.1 PEST Calibration Process
TFIELD5.3.2 Calibrated Correction of SteadyState Head Values
TFIELD5.3.2.1 Localized Recalibration in the Vicinity of SNL8
TFIELD5.3.2.2 Continued Recalibration Activity
TFIELD5.3.3 Evaluation of Impact of Multiple Calibration Processes
TFIELD5.3.4 Selection of BestCalibrated Fields
TFIELD5.4 Step 3  PostCalibration Analysis
TFIELD5.4.1 Statistical Analysis of Resulting T, A, and S Fields
TFIELD5.4.1.1 Final Transmissivity and Anisotropy Fields
TFIELD5.4.1.2 Final Storativity Values
TFIELD5.4.1.3 Final Recharge Values
TFIELD5.4.2 Forward Model Results Using the Calibrated Fields
TFIELD6.0 Culebra TField Mining Modifications
TFIELD6.2 Model Domain, Boundary, and Initial Conditions
TFIELD6.2.1 Boundary and Initial Conditions
TFIELD6.2.2 Determination of Potential Mining Areas
TFIELD6.2.3 Use of Mining Zones in Forward Simulations
TFIELD6.2.4 ParticleTracking Simulations
TFIELD6.3 ParticleTracking Results
TFIELD6.3.1 Particle Travel Times
TFIELD6.3.4 ParticleTracking Discussion
Figure TFIELD 21. Generalized Stratigraphy Near the WIPP
Figure TFIELD 23. Rustler Formation Stratigraphic Nomenclature
Figure TFIELD 24. M1/H1 Halite Margin In the Lower Los Medaños Member
Figure TFIELD 25. M2/H2 Halite Margin In the Upper Los Medaños Member
Figure TFIELD 26. M3/H3 Halite Margin In the Tamarisk Member
Figure TFIELD 27. M4/H4 Halite Margin In the Fortyniner Member
Figure TFIELD 47. Standard Deviation of Indicator Values Across All 1000 Base Realizations
Figure TFIELD 48. Sample Log10 T (m2/s) Base Field Realization r123
Figure TFIELD 49. Mean Log10 T (m2/s) Values Across All 1000 Base Realizations
Figure TFIELD 51. Complete Calibration Process for a Single Realization
Figure TFIELD 512. Flowchart Illustrating the PEST Calibration Process
Figure TFIELD 518. HighT Zone Membership Calculated for the Calibrated T Values
Figure TFIELD 519. Number of T Fields Where Low T Became High T Through PEST Calibration
Figure TFIELD 520. Number of T Fields Where High T Became Low T Through PEST Calibration
Table TFIELD 41. βvalues for Regression Equation TFIELD 4.1
Table TFIELD 51. Freshwater Head Observations Used as Steadystate Calibration Targets
Table TFIELD 52. Summary of Transient Observations Used as Calibration Targets
Table TFIELD 53. Initial Values of Parameters at Pilot Points
Table TFIELD 56. Cutoff Values for Final Field Selection
Table TFIELD 57. Final Selected Field Identifiers
A transmissivity anisotropy [dimensionless]
AMSL above mean sea level
AP analysis plan
BLM Bureau of Land Management
CB Cabin Baby
CCA Compliance Certification Application
CDF cumulative distribution function
CFR Code of Federal Regulations
CRA Compliance Recertification Application
D diffusivity [m^{2}/s]
DOE U.S. Department of Energy
EPA U.S. Environmental Protection Agency
km kilometer
LWB Land Withdrawal Boundary
m meter
m^{2} square meters
M/H mudstone/halite margin
m^{2}/s square meters per second
m^{3}/s cubic meters per second
m/yr meters per year
PA performance assessment
PABC performance assessment baseline calculation
R recharge [m/s]
RMSE root mean squared error
S storativity [dimensionless]
SNL Sandia National Laboratories
SSE sum of squared errors
SVD singular value decomposition
T transmissivity [m^{2}/s]
USGS U.S. Geological Survey
UTM Universal Transverse Mercator map projection
WIPP Waste Isolation Pilot Plant
Modeling the transport of radionuclides through the Culebra Dolomite Member of the Rustler Formation (hereafter referred to as the Culebra) is one component of the performance assessment (PA) performed for the U.S. Department of Energy (DOE) Waste Isolation Pilot Plant (WIPP) 2014 Compliance Recertification Application (CRA). Transport modeling in PA requires flow velocity results from the Culebra groundwater flow model. This Appendix describes the process used to develop and calibrate the input parameter fields for the Culebra flow model. Calibrated model parameters are referred to broadly as "Tfields" (transmissivity fields), although more parameters than just transmissivity (T) were calibrated as part of the CRA2009 Performance Assessment Baseline Calculation (PABC) model (Clayton et al. 2010). This appendix describes the process followed for the CRA2009 (PABC), which was a major change from the process followed for CRA2004 PABC (Leigh et al. 2004), and involved a hydrology conceptual model peer review. The Tfields developed for CRA2009 PABC were used unchanged in CRA2014. Figures illustrating each calibrated Tfield are given in Attachment A to this Appendix.
The work described in this Appendix was performed under two Sandia National Laboratories (SNL) analysis plans (APs): AP114 (Beauheim 2008) and AP144 (Kuhlman 2009). AP114 (evaluation and recalibration of Culebra transmissivity fields) dealt with the development and calibration of the Tfields (including T, storativity (S), horizontal anisotropy (A), and vertical recharge (R)), in addition to development of Tfield acceptance criteria. AP144 (calculation of Culebra flow and transport) dealt with the modification of Tfields for the potential future effects of potash mining for use in the PA Culebra radionuclide transport calculations. The PA Culebra radionuclide transport calculations are not described in this Appendix, which focuses on the development and modification of the Tfields.
West of the WIPP, Culebra T is high where the Culebra overlies areas where the Salado Formation has been removed by dissolution (mostly in Nash Draw). East of the WIPP, Culebra T is low when the Culebra is bounded either above or below by halite in adjoining Rustler units. Further to the east, Culebra T is very low when the Culebra is bounded both above and below by halite in the Rustler. At the WIPP, between the high T in the west and low T in the east, Culebra T is observed to change significantly over short distances and is simulated in the WIPP Culebra flow model using a random mixture (i.e., stochastic patches) of high and low T zones, consistent with geologic and hydrologic observations. The geologic data discussed in Section TFIELD2.0 are used to specify the boundaries of these Culebra conceptual model zones (Section TFIELD3.0), which are then carried forward into the numerical implementation of the Culebra groundwater model (Section TFIELD4.0).
The starting point in the Tfield development process was to assemble and update information on geologic factors potentially affecting Culebra T (Section TFIELD2.0). These factors include dissolution of the upper Salado Formation located below the Culebra, presence of gypsum cements, the thickness of overburden above the Culebra, and the spatial distribution of halite in the Rustler Formation both above and below the Culebra. Geologic information is available from hundreds of oil and gas wells and potash exploration holes in the vicinity of the WIPP site, while estimates of Culebra T are available from only 64 well locations. Details of the geologic data compilation are given in Powers (Powers 2002a, Powers 2002b and Powers 2003), updated in Powers (Powers 2007a and Powers 2007b), and summarized in Section TFIELD2.0.
A twopart geologically based approach was used to generate base Culebra Tfields. In the first part (Section TFIELD3.0), a conceptual model for geologic controls (i.e., soft data) on Culebra T was formalized and the hypothesized geologic controls were regressed against Culebra T estimates at monitoring wells to determine linear regression coefficients. The regression includes one continuously varying function, Culebra overburden thickness, and three indicator functions that assume values of 0 or 1 depending on the occurrence of open, interconnected fractures; Salado dissolution; and the presence or absence of halite in Rustler units bounding the Culebra.
In the second part (Section TFIELD4.0), a method was developed for applying the linear regression model to predict Culebra T across the WIPP model area between the sparse observations at wells. The regression model was combined with the maps of geologic factors to create 1,000 stochastically varying base Culebra Tfields. Details about the development of the regression model and the creation of the base Tfields are given in Hart et al. (Hart et al. 2008). The conceptual model embodied in these 1,000 base Culebra Tfields was subject to peer review before model calibration proceeded (Section TFIELD3.7). The peer review panel concluded the justification and scientific rigor of the methodology for preparing base Tfields were adequate (Burgess et al. 2008). A sample of 200 out of the 1,000 created base Tfields were calibrated following the process outlined in Section TFIELD5.0, with the 100 best calibrated Tfields eventually chosen for use in PA radionuclide transport calculations.
Section TFIELD5.0 presents details on the modeling approach used to calibrate the Tfields to both steadystate heads across the model domain and transient drawdown measurements from multiwell pumping tests. Heads measured in 42 Culebra observation wells around May 2007 were used to represent steadystate conditions in the Culebra, and drawdown responses in 67 total observation wells (62 unique locations) across nine pumping tests were used to provide transient calibration data. See Appendix HYDRO2014 for more information on the Culebra monitoring well network and recent trends observed in Culebra water levels. Details on the steadystate heads are described in Johnson (Johnson 2009a and Johnson 2009b), and the transient drawdown data are summarized in Hart et al. (Hart et al. 2009). Assumptions made in modeling, the definition of an initial head distribution, assignment of boundary conditions, discretization of the spatial and temporal domain, weighting of the observations, and the use of PEST (Doherty 2000) with MODFLOW2000 to calibrate the Tfields using a pilotpoint method are described in detail in Hart et al. (Hart et al. 2009) and summarized in Section TFIELD5.0. Section TFIELD5.3.4 addresses the development and application of acceptance criteria to select the 100 best Tfields from the 200 calibrated Tfields. Acceptance was based on a combination of objective fit to both the steadystate and transient pumping test calibration data. Section TFIELD5.4 provides summary statistics and other information for the 100 Tfields that were judged to be acceptably calibrated. Attachment A presents T, S, diffusivity (D), and modelpredicted flow speed for each of the chosen 100 realizations.
The data used in the construction (Sections TFIELD3.0 and TFIELD4.0) and calibration (Section TFIELD5.0) of the Tfields were divided into three groups:
This first group included geologic data (i.e., Salado dissolution, presence of gypsum cements in the Culebra, Culebra overburden thickness, and the locations of Rustler halite margins) and hydrologic data (i.e., pairs of pumping and observation wells interpreted to have high diffusivity from multiwell pumping tests). The soft data were included in Tfield generation because they were used to define boundaries also used in the calibration phase. The second group only included estimated T from singlewell pumping tests, which were used directly in indicator kriging, indirectly in the overburden regression analysis, and directly as fixed pilot points in the calibration phase. The third group only included head and drawdown observations used as calibration targets. Some of these third group data were used in other analyses to estimate diffusivity values used as soft data (first group), but the drawdown data from multiwell pumping tests only appeared directly in the calibration phase.
Section TFIELD6.0 discusses modifications of the Tfields performed to account for the effects of potash mining both within and outside the WIPP land withdrawal boundary. Potentially miningaffected areas were delineated, random transmissivity multipliers were applied to the transmissivity field in those areas, and particle tracks and travel times were computed (Kuhlman 2010). The flow fields produced by these miningaffected Tfields were input to the radionuclide transport model SECOTP2D used to compute both CRA2009 PABC and CRA2014 longterm PA releases (Appendix PA2014). Section TFIELD7.0 provides an executive summary of the development and modification of the Culebra Tfields.
The work outlined in Section TFIELD2.0 was performed as Task 1 under AP114, Analysis Plan for Evaluation and Recalibration of Culebra Transmissivity Fields (Beauheim 2008). There were no changes to the model between CRA2009 PABC and CRA2014. Geologic data were updated to improve definition of geologic boundaries used to define zones as part of the process of creating new Tfields. Geologic boundaries were refined for CRA2009 PABC using data from field investigations and newly obtained oil and gas well log data (Section TFIELD2.2). The Salado dissolution margin bounds the highT Culebra zone to the west in Nash Draw, and was only modified slightly for CRA2009 PABC. The Rustler halite margins bound the very lowT Culebra zone to the east, and were modified significantly for CRA2009 PABC. The confinement of the Culebra in the southeastern portion of Nash Draw was also investigated in AP114 Task 1, to constrain the Culebra flow model inputs (Section TFIELD2.3). Previous to CRA2009 PABC, the Culebra groundwater flow model was only steadystate and did not include input parameters related to confinement or recharge. A model was developed regarding distribution of gypsum cements in the Culebra from available core data (Section TFIELD2.3.3).
The Culebra Member of the Rustler Formation is considered as a potential longterm release pathway in WIPP PA because it is the most permeable laterally continuous geologic unit above the WIPP repository level (see Figure TFIELD 21 for general stratigraphy). Potential future human intrusion into the repository might connect the repository with the Culebra, which would then transport radionuclides to the accessible environment under natural flow conditions. The accessible environment is defined to be where the WIPP Land Withdrawal Boundary (LWB) intersects the Culebra in the subsurface.
The ability of the Culebra to advect groundwater and radionuclides to the accessible environment is affected by both depositional and postdepositional effects. The Culebra is believed to have been deposited quite uniformly over a wide area (the vertical thickness of the Culebra is quite uniform over lateral distances of many miles (e.g., Holt and Powers 1988)), but depositional effects include the presence of mudstone or halite layers in the Rustler Formation immediately above and below the Culebra. Postdepositional processes include dissolution of halite from the underlying Salado Formation and precipitation of vug and porefilling evaporates within the Culebra (see Figure TFIELD 22).
Understanding of the spatial distribution and thickness of halite in the Rustler Formation was improved for CRA2009 PABC (compared to CRA2009 (U.S. DOE 2009)) due to data obtained from analysis of geophysical logs from oil and gas wells. The lateral extent of Salado dissolution was modified slightly for CRA2009 PABC, but remained largely similar to CRA2009, with minor adjustments due to additional information for a few new WIPP wells (Powers 2007a and Powers 2007b).
Groundwater flow in the Culebra is generally from north to south at the WIPP site. Water levels in Culebra wells in Nash Draw (several miles to the west of the WIPP site) respond more rapidly to precipitation and behave differently than Culebra wells in the immediate vicinity of the WIPP site (Appendix HYDRO2014, Section 7.1 ). Based on the northsouth gradient currently observed at the WIPP, particletracking predictions from the WIPP waste panels through the Culebra result in flow towards the southern edge of the WIPP LWB.
Figure TFIELD 21. Generalized Stratigraphy Near the WIPP
The locations of the Rustler halite margins (Section TFIELD2.2.1) and the Salado dissolution margin (Section TFIELD2.2.2) both affect the conceptual model of Culebra T (see Figure TFIELD 21 for general relationships between the Rustler, the Salado and the WIPP repository). These were updated as part of the CRA2009 PABC geologic study.
The presence of halite in the nondolomite members of the Rustler Formation correlates strongly with estimates of Culebra T. Halite and anhydrite are found as porefilling cements in the Culebra (reducing open fractures) when halite exists in layers above the Culebra, below the Culebra, or both (Figure TFIELD 22). When halite is found either above (H3) or below (H2) the Culebra, observed Culebra T is low. When halite exists both above and below the Culebra, observed Culebra T is extremely low.
North, south and west of the WIPP site, Cenozoic dissolution has affected the upper Salado Formation. Where this dissolution has occurred, the rocks overlying the Salado, including the Culebra, are strained (leading to larger apertures in existing fractures), fractured, collapsed, or brecciated (e.g., Beauheim and Holt 1990). All WIPP Culebra wells within the Salado dissolution zone have been interpreted to have high T. It is hypothesized that all regions affected by Salado dissolution have wellinterconnected fractures and therefore high T.
Figure TFIELD 22. WIPP Culebra Dolomite Conceptual Model. Culebra T decreases to the east (increasing overburden and halite) and increases to the west (fracturing due to underlying Salado dissolution). Halite appears both above (H3) and below (H2) the Culebra in the east. Primary groundwater flow direction through the Culebra is south.
The Rustler Formation stratigraphic column given in Figure TFIELD 23 shows two types of geologic variability. Vertical stratigraphy places older formations below younger formations at the same location in space (e.g., the Los Medaños Member is older than the Culebra Member), while facies change place two units of similar age at different spatial locations, due to changes in depositional environments (e.g., Mudstone 4 (M4) and Halite 4 (H4) of the Fortyniner Member are of the same age, but occur in different locations).
Figure TFIELD 23. Rustler Formation Stratigraphic Nomenclature
Powers (Powers 2002a and Powers 2002b) provided geologic data across the CRA2004 PABC Culebra modeling domain that included maps of halite margins within the Rustler Formation. Those margins were largely based on work in Powers and Holt (Powers and Holt 1995), modified by some data collected from potash drillholes, especially in the northern area of the Culebra modeling domain. The observed distribution and thickness of halite in the Rustler is interpreted to be the result of sedimentary structures and facies relationships controlled by deposition, rather than the result of dissolution alone (Holt and Powers 1988; Powers and Holt 1999; Powers and Holt 2000). Before Holt and Powers (Holt and Powers 1988), many researchers incorrectly believed a uniform thickness of Rustler halite was deposited and later removed by dissolution in the areas near Nash Draw, leaving the observed mudstone layers as dissolution residue. Definitive data collected during WIPP airintake shaft geologic mapping provided the basis for the current faciesbased conceptual model used at the WIPP (Holt and Powers 1988). Some minor zones adjacent to the depositional margins have been interpreted as having undergone some postdepositional dissolution of halite, specifically the halite in the Tamarisk Member, but the extent of this Rustler halite dissolution is relatively minor (Beauheim and Holt 1990).
Significant changes to the locations of the M3/H3 and M2/H2 margins have been made for CRA2009 in some areas since CRA2004 (U.S. DOE 2004) as part of Task 1A of AP114. The Rustler halite margins used since CRA2009 PABC are shown over a wide area in Figure TFIELD 24 through Figure TFIELD 27, as defined in Powers (Powers 2007a and Powers 2007b). Changes in the location of the halite margins were based mostly on newly obtained geophysical log data obtained from oil and gas exploration (both new and old wells), and a few hydrologic wells drilled by the WIPP program since 2003 (Powers 2007a).
Wells H17 and H12 (see Figure TFIELD 28), located where halite occurs in the Tamarisk Member (H3 interval; Figure TFIELD 26) but not in the Los Medaños Member (M2 interval; Figure TFIELD 25) of the Rustler Formation show low transmissivity. We assume hightransmissivity zones do not occur in the Culebra where H2 or H3 also occur. Margins near the WIPP remain nearly unchanged, and all modifications to the margins do not change the basic interpretation that the margins are the result of deposition and local syndepositional dissolution of halite, not regional halite dissolution from the Rustler (Holt and Powers 1988; Powers and Holt 2000; Powers et al. 2006). Core evidence from well SNL8 shows limited brecciation of anhydrite 3 in the Tamarisk (Figure TFIELD 23) that is interpreted as an extension of a narrow margin along the H3 margin where a limited amount of halite was dissolved after deposition (Powers 2009).
After refining the Rustler halite margin locations, all mudstone/halite margins now show similar gross trends (compare Figure TFIELD 24 through Figure TFIELD 27 and Figure TFIELD 28). Southeast of the WIPP, the margins are elongate roughly northwest to southeast. The gross trends of these margins are similar to the trend in the elevation of the top of Culebra (Figure TFIELD 29). As previously described (e.g., Holt and Powers 1988); Powers et al. (Powers et al. 2003), this northwesttosoutheast trending anticlinal feature is called the Divide anticline. Mudstone dominates along this trend in three of the mudstonehalite units of the Rustler (i.e., all except M1/H1).
Figure TFIELD 24. M1/H1 Halite Margin In the Lower Los Medaños Member 
Figure TFIELD 25. M2/H2 Halite Margin In the Upper Los Medaños Member 
Figure TFIELD 26. M3/H3 Halite Margin In the Tamarisk Member 
Figure TFIELD 27. M4/H4 Halite Margin In the Fortyniner Member 
Figure TFIELD 28. Salado Dissolution Margin and Rustler Mudstone/Halite (M/H) Margins. WIPP Culebra wells with high or low transmissivity (T) are indicated. WIPP Culebra model extents indicated with large black rectangle. Wells mentioned in text are labeled using larger font.
Figure TFIELD 29. Top Elevation (m Above Mean Sea Level (AMSL)) of the Culebra. WIPP LWB indicated with blue dashed line. Township (T) and Range (R) corners indicated with crosses.
A margin marking the lateral extent of significant dissolution of upper Salado Formation halite for CRA2004 PABC was inferred from significant local changes in thickness of the interval between the Culebra Dolomite and the Vaca Triste Sandstone Member of the Salado (Powers 2002a and Powers 2002b). For CRA2009 PABC, the margin was modified to reflect information indicating embayments of the dissolution margin. Additional data were added south of the WIPP, with log cross sections, to delineate the margin more accurately (Powers 2003). Some of these data are reflected in the simplified maps included in Powers et al. (Powers et al. 2003) and Holt et al. (Holt et al. 2005).
The Salado dissolution margin was updated for CRA2009 (see Appendix G from the Analysis Report for Task 5 of AP114 (Hart et al. 2008)) based on reinterpretation of geophysical borehole logs from oil and gas wells in the vicinity of H9, which were not available for CRA2004 PABC. This analysis placed H9 east of the dissolution line, where previously it was considered to be within the area affected by Salado dissolution. The Salado dissolution margin is shown with the Rustler halite margins in Figure TFIELD 28, reflecting the change near H9.
Field and map studies were performed to identify potential recharge locations south and west of the WIPP in the southeastern arm of Nash Draw (Powers 2006). This work also identified Culebra unconfined regions in the same geographic area. The boundaries to the west and south correspond to the model domain; the northern and eastern boundaries included the southeastern arm of Nash Draw and an area beyond the apparent eastern extent of the draw.
Five elements were identified as contributing to understanding recharge, which might be useful for modeling the possible effects of recharge to the Culebra in the study area:
1. extent of and relationship between surface drainage basins,
2. areas with differing Culebra confinement,
3. location and character of drainage channels within drainage basins,
4. location of specific recharge points (e.g., sinkholes), and
5. soil characteristics and rainfall infiltration across the study area.
Of these, the estimate of Culebra confinement is the most interpretive element. Drainage basins, channels and specific points of recharge are identified using surface topography features identifiable from maps, aerial photos, or field reconnaissance. Existing maps of soils, combined with surface reconnaissance and aerial photographs, permitted relatively direct assignment of soil properties controlling runoff. The degree of confinement of the Culebra in the study area, however, was not directly determinable from the surface data. As a result, a variety of surface features and well data were combined to estimate areas where the Culebra is less confined compared to conditions at the WIPP site, where there are more welltest and drillhole data.
Drainage basin size and characteristics are important elements to determine how rainfall, infiltration, and runoff may contribute to recharge of nearsurface Rustler hydrologic units in Nash Draw. Topographic maps, aerial photographs, and some field checking were used to define separate surface drainage basins.
The drainage basins are mainly separated by topographic divides and local lows or concentration points that can be distinguished on 7.5minute quadrangle topographic maps supplemented by study of aerial photographs. Because Nash Draw is an area of significant evaporite karst (e.g., Powers and Owsley (Powers and Owsley 2003)), collapse features, caves, or sinkholes may capture local drainage in smaller basins or subbasins wholly enclosed by another basin (Figure TFIELD 210). An example is drainage basin 7, which is wholly enclosed in drainage basin 6 (Figure TFIELD 210). These quite localized closed drainage basins in Nash Draw represent potential recharge locations for the Rustler Formation. Mapping the basins is the first step in understanding the complex geology and hydrology inside Nash Draw, which expresses itself as waterlevel fluctuations in some Culebra wells in and near Nash Draw (see hydrographs and references in Appendix HYDRO2014).
Figure TFIELD 210. Closed Drainage Subbasins Identified in Southeastern Nash Draw. White areas are either outside Nash Draw or the study area.
Across the WIPP site, the Culebra can be considered confined, with little potential for direct vertical recharge for the relatively short time period covered in the WIPP Culebra model calibration (i.e., the length of multiwell pumping tests). Within portions of Nash Draw, the Culebra is very shallow (i.e., only covered by portions of the highly fractured Tamarisk) and observed water levels show the Culebra responds to precipitation events in a very short time (Appendix HYDRO2014). Due to the interpretative nature of the confinement estimate, no numerical values of storativity were predicted from this geologic analysis, only zones of relatively higher or lower confinment, with an intermediate transition zone. The confined area has a relatively unambiguous definition, whereas the boundary between transition and unconfined is much more subjective.
The area of the Culebra considered confined (red in Figure TFIELD 211) is defined approximately by the interpreted margin of upper Salado halite dissolution (Powers et al. 2003). There is a significant increase in Culebra T values west and south of this margin, and this change is attributed to changes in fracture aperture associated with strain induced by dissolution. The transition zone (green in Figure TFIELD 211) includes areas where some data from wells indicate there is some vertical isolation of the Culebra, but information is less conclusive.
Most of the Culebra unconfined zone (blue in Figure TFIELD 211) is in central Nash Draw and out of the AP114 Task 1 study area. The strategy for estimating relative Culebra confinement was to select areas where the Culebra is known or believed to be very shallow (≤30 meters (m) below ground surface) and where observed recharge points (caves, sinkholes, alluvial dolines) are believed to access units below the Magenta. Some large caves and sinkholes are developed in the Tamarisk gypsum beds and have a greater likelihood of providing hydraulic connection to the Culebra than similar openings in the Fortyniner gypsum beds. Many potash exploration holes within Nash Draw encountered lostcirculation zones, but the stratigraphic relationships of these zones to the Culebra are not well constrained. Thus, apart from the location from Livingston Ridge (the escarpment marking the eastern edge of the surface expression of Nash Draw) and the upper Salado dissolution margin, the factors determining confinement of the Culebra are generally qualitative.
Figure TFIELD 211. Culebra Confinement Map for Southern Nash Draw Study Area. White areas are outside the Nash Draw geologic study area. Zones are shown over the entire model area in Figure TFIELD 53.
The amount of gypsum cements in fractures and vuggy porosity within the Culebra is believed to be inversely related to Culebra T (Beauheim and Holt 1990). They postulated gypsum fracture fillings limited Culebra T by closing fracture apertures, filling critical fracture junctions. The postulated relationship remained qualitative because too few well locations had both measured T values and describable core. Since 1990, the Culebra has been cored and hydraulically tested at 24 additional locations, providing sufficient data to construct a quantitative model linking Culebra T with the presence of gypsum cements. No soft data on gypsum cements was used in Tfield construction or calibration before CRA2009 PABC.
In Appendix F of Hart et al. (Hart et al. 2008), a simple quantitative model was constructed relating Culebra gypsum content to T. Using units defined by Holt (Holt 1997), maps were developed to illustrate spatial occurrence of gypsum in the Culebra. The maps used a gypsum index accounting for the relative Culebra gypsum content (Figure TFIELD 212 and Figure TFIELD 213). Using a critical value of the gypsum index, the highT/lowT status of Culebra well locations were predicted with an accuracy >97% for WIPP well locations where both sufficient core and T estimates exist. These maps revealed that regions of no gypsum occur predominantly where Salado dissolution has affected the Culebra. The lowgypsum region in the southern WIPP LWB (Figure TFIELD 213) is similar to the highdiffusivity region defined by Beauheim (Beauheim 2007) (Figure TFIELD 42). Soft data were used to incorporate information about the influence of gypsum content on predicted Culebra T.
Figure TFIELD 212. Areas Where No Gypsum Has Been Found in Core Samples, Corresponding to a Greater Likelihood of Having Higher Culebra T Values
Figure TFIELD 213. Areas Where Wells Have Either No or Low Gypsum Content. The areas not shaded are likely to have high gypsum content and lower T.
The work outlined in Section TFIELD3.0 was performed for CRA2009 PABC under AP114, Analysis Plan for Evaluation and Recalibration of Culebra Transmissivity Fields (Beauheim 2008), and still applies to CRA2014. The conceptual model for base field creation was originally explained in Holt and Yarbrough (Holt and Yarbrough 2002), as Task 2, Subtask 1 of AP088 for CRA2004 PABC. Since then, the data supporting the conceptual model were updated and improved, but the model itself has changed very little. Any deviations of the CRA2009 PABC model from the CRA2009 model due to updates in data or process are discussed in this section. No updates have occurred for CRA2014 since CRA2009 PABC.
Figure TFIELD 22 illustrates the current geologic and hydrologic conceptual model of the Culebra dolomite in the vicinity of the WIPP site. Geologic controls on Culebra T were identified and a linear regression model relating these controls to T was constructed. The geology and geologic history of the Culebra has been described in detail elsewhere in the literature (Holt and Powers 1988; Beauheim and Holt 1990; Holt 1997). The following conceptual model was developed from this published work. Specifically, the model follows Holt (Holt 1997) in assuming variability in Culebra T is due strictly to postdepositional processes. Throughout the following discussion, the informal stratigraphic subdivisions of Holt and Powers (Powers 1988) are used to identify geologic units within the Rustler Formation, as listed in Figure TFIELD 23 and shown in map view for the Culebra model area in Figure TFIELD 28. The Culebra conceptual model given in this section passed a peer review (Burgess et al. 2008) before the calibration process in Section TFIELD5.0 was begun.
It is hypothesized that Culebra T spatial distribution is a function of several geologic factors, some of which can be determined at a location using mapped geologic data, including:
1. Culebra overburden thickness,
2. fracture interconnection,
3. presence of gypsum cements in fractures and vuggy porosity,
4. dissolution of the upper Salado Formation below the Culebra, and
5. occurrence of halite in Rustler units above or below the Culebra.
HighT regions near the WIPP cannot be predicted using geologic data, as they represent areally persistent zones of wellinterconnected fractures, and fracture interconnection cannot be observed or inferred from core or geophysical log data. Fracture interconnection is therefore treated as a stochastic process. Presence of gypsum cements in the Culebra, occurrence of Rustler halite, and Culebra overburden thickness instead varies slowly in space. These properties can be meaningfully mapped at the scale of the groundwater flow model.
The CRA2009 PABC model domain was expanded to the east relative to the domain used for the CRA2004 (U.S. DOE 2004) to reach an area where halite is present in all of the nondolomite members of the Rustler Formation. This change was made to simplify the specification of the eastern boundary condition of the model. The current extent of the model domain is 601,700 to 630,000 m UTM X NAD27 and 3,566,500 to 3,597,100 m UTM Y NAD27. The domain was discretized into 100m square cells, yielding a model 284 cells wide by 307 cells high. The Culebra was modeled as a single layer of uniform 7.75m thickness (U.S. DOE 1996). The area covered by Figure TFIELD 28 corresponds to the model domain, showing the WIPP site boundary, the relevant geologic margins, and various Culebra monitoring wells. The model domain for CRA2014 has not changed since CRA2009 PABC.
An inverse relationship was hypothesized between Culebra overburden thickness and Culebra T. Overburden thickness is a metric for two different controls on Culebra T. First, fracture apertures can be related to overburden thickness (e.g., Currie and Nwachukwu 1974), as lower T are found where Culebra depths are greater (Beauheim and Holt 1990; Holt 1997). Second, erosion of overburden leads to stressrelief fractures, and the amount of Culebra fracturing increases as the overburden thickness decreases (Holt 1997). The structure contour map of Culebra elevation (Figure TFIELD 29) has been constructed using geophysical logs from hundreds of oil and gas wells, and geologic information from more than 100 WIPPrelated boreholes. The difference between the land surface elevation (as obtained from U.S. Geological Survey (USGS) topographic maps) and Culebra elevation is the overburden thickness (Figure TFIELD 31). Culebra overburden thickness ranges from near zero in the southern end of Nash Draw, to over 550 m in the northeastern corner of Figure TFIELD 31. The depth to the Culebra from the land surface defined the value of d(x,y) for each cell in the Culebra model domain.
Figure TFIELD 31. Culebra Overburden Thickness Contours (m). Square is the WIPP LWB; irregular black outline west of WIPP is Nash Draw.
HighT zones within the Culebra are associated with interconnected fractures and occur randomly between areas bounded on the west by the Salado dissolution margin and on the east by H2 and/or H3 (the central area Zone 4 in Figure TFIELD 32). In these zones, fractures are wellinterconnected, and fracture interconnectivity is controlled by a complicated history of fracturing with several episodes of cement precipitation and dissolution (Beauheim and Holt 1990; Holt 1997). Unfortunately, no geologic metric for fracture interconnectivity was identifiable in cores or from subsurface geophysical logs, and fracture interconnectivity has only been identified from in situ hydraulic test data.
Because of this lack of a corresponding easytomap geologic metric for fracture interconnectivity, the spatial location of highT zones was considered to be a stochastic process that could not be predicted deterministically. The spatial layout of these zones was simulated for CRA2009 PABC using geostatistical indicator kriging with conditioning data (this was not changed for CRA2014 since CRA2009 PABC). This stochastic development of zones was a change from CRA2004 PABC (Holt and Yarbrough 2002), where the only conditioning information was based on the T at wells. Information was added to the geostatistical model to increase the likelihood of high T being placed between two wells that hydraulic testing has revealed to be associated with larger diffusivity values. North of the WIPP site (i.e., south of WIPP30) evidence exists for both high levels of gypsum in the Culebra and relatively high D between pumping/observation well pairs. In this unique region, the geologic conceptual model indicates there is slightly lower probability of being in a highT zone than in other areas where a high D or high T estimate exists. Section TFIELD4.2 discussed the process of merging hydraulic hard and soft data (singlewell T estimates and multiwell D estimates, respectively) with geologic soft data on gypsum.
Figure TFIELD 32. Conceptual Model Zones With Indicator Values and Zone Numbers (Equation TFIELD 3.2). Zones 3 and 4 are distributed randomly between the Salado dissolution margin and westernmost M2/H3 or M3/H3 Rustler halite margins.
The Culebra T estimates at WIPP wells used in the CRA2009 PABC modeling were the same as those used by Holt and Yarbrough (Holt and Yarbrough 2002), supplemented by more recent data reported from subsequent pumping tests (Roberts 2006 and 2007; Bowman and Roberts 2008). The log_{10} T data show a bimodal distribution in Figure TFIELD 33. Closely spaced wells sometimes show very different values; higherT values are hypothesized to reflect the presence of wellinterconnected fractures absent at lowerT locations. For example, wells WQSP2 and WIPP12 are only 454 m apart, but have T values differing by over two orders of magnitude (see blue star labeled W12 and red circle labeled WQSP2 in the north portion of the WIPP LWB in Figure TFIELD 28). Thus, the fractures present at WQSP2 apparently do not extend to WIPP12 or are not intersected by the WIPP12 borehole. Wellinterconnected fractures occur in regions affected by Salado dissolution (e.g., Nash Draw) and in areas with complicated cement dissolution and precipitation histories (e.g., highT zones near the WIPP site). The natural break between the measured log_{10} T square meters per second (m^{2}/s) values at −5.4 (Holt and Yarbrough 2002) is illustrated with a vertical black line in Figure TFIELD 33. The fractureinterconnection indicator (I_{f} ) is defined in terms of this break (Equation TFIELD 3.1).
Figure TFIELD 33. Histogram of Log_{10} Culebra Transmissivity (T) Estimates at WIPP Wells from Singlewell Tests
Slight modification was made to the Salado dissolution margin used in CRA2009 PABC, compared to CRA2009, as outlined in Section TFIELD2.2.2. No modifications were made for CRA2014 since CRA2009 PABC. The indicator variable for Salado dissolution is I_{D} , and was defined to be 1 in areas of the model domain where dissolution has occurred, and 0 elsewhere. The Salado dissolution margin is plotted with the Rustler halite margins in Figure TFIELD 28.
The M2/H2 and M3/H3 Rustler halite margins were modified for CRA2009 PABC compared to CRA2009, as outlined in Section TFIELD2.2.1. No modifications were made for CRA2014 since CRA2009 PABC. The margins are shown individually in Figure TFIELD 25 and Figure TFIELD 26, and together with the M1/H1 and M4/H4 Rustler halite margins and Salado dissolution margin in Figure TFIELD 28.
Wells SNL6 and SNL15 were drilled since Holt and Yarbrough (Holt and Yarbrough 2002). They are located east of the M2/H2 and M3/H3 halite margins, where halite is present in both intervals (see Figure TFIELD 28). As predicted by Holt (Holt 1997), the Culebra itself was partially cemented with halite at these locations, and estimated T were extremely low (Roberts 2007; Bowman and Roberts 2008). Based on these observations, Culebra T is assumed lower in the region where halite occurs both above (in the M3/H3 interval) and below (in the M2/H2 interval), than the Culebra T where halite occurs in only one of these intervals. The indicator term I_{H} was defined to be 1 at any point where halite is present in both the M2/H2 and M3/H3 margins, and to be 0 elsewhere.
The following linear model for Y(x,y) = log_{10} T(x,y) was constructed
where β _{1} through β _{5} are regression coefficients, the twodimensional location vector (x,y) consists of NAD27 UTM Zone 13 x and y coordinates, d(x,y) is the Culebra overburden thickness (Figure TFIELD 31), I_{f} is an indicator of whether interconnected fractures are present in the Culebra, I_{D} is the Salado dissolution indicator, and I_{H} is the halite bounding indicator. In this model,  means logical or, while & means logical and. Regression coefficient β _{1} is the intercept value for the linear model. Coefficient β _{2} is the slope of Y(x,y)/d(x,y). The coefficients β _{3}, β _{4}, and β _{5} represent adjustments to the intercept for the occurrence of interconnected fractures, Salado dissolution, and halite bounding, respectively. Although other types of linear models could have been developed, Equation TFIELD 3.2 is consistent with the conceptual model relating Culebra T to geologic controls, can be tested using published WIPP geologic and T estimates, and can be potentially verified with new Culebra wells.
Because there are only two data points for T in the zone where Culebra is bounded by halite, and both are significantly lower than any other T values in the model, the β _{5} I_{H} term in Equation TFIELD 3.2 was included to take into account the very low T zone. This was done to keep the conceptual model consistent for all zones, recognizing the base fields are primarily a starting point for subsequent calibration.
The combined results of the regression and the indicator kriging (Section TFIELD4.3) were 1,000 base Tfields that shared certain geologic features, but were different from one another. This difference was provided by the stochastic placement of highT areas in the central zone. These areas were placed using the GSLIB Sequential Indicator Simulation (SISIM) routine (qualified for use in WIPP PA according to NP 191 (Chavez 2006)). This routine used geostatistical methods to create stochastic indicator (Boolean value) fields.
The Culebra conceptual model given in this section passed peer review before proceeding with the CRA2009 PABC calibration of the Culebra T fields (Burgess et al. 2008). The peer review panel found the methodology presented here to be adequate, accurate, and valid enough to justify proceeding with the numerical implementation and calibration of the Culebra Tfields. The panel found the CRA2009 PABC conceptual model to be greatly improved, compared to the Culebra conceptual model used in the Compliance Certification Application (CCA) (U.S. DOE 1996). The panel found the understanding of the physical processes connecting the Culebra groundwater geochemistry with the Culebra hydraulic properties to be insufficient. The peer review panel did not feel this particular lack of understanding would be a problem in Tfield development and calibration, due to the relatively high density of Culebra hydrologic data available at the WIPP site.
The work outlined in Section TFIELD4.0 was performed under AP114, Analysis Plan for Evaluation and Recalibration of Culebra T Fields (Beauheim 2008). This section discusses details associated with the incorporation of soft and hard data into the base Tfield construction process. The base Culebra Tfields were the starting point for the calibration processes outlined in Section TFIELD5.0. Aside from the definitions of some fixed parameter zones, all the parameters specified in the base Tfield construction were allowed to be modified during the calibration process to produce model output that better matched observed steadystate and transient pressure observations. Inside the WIPP LWB there is a large amount of hard data to constrain the parameters of the groundwater model during calibration, while distant to the WIPP LWB the hard data are not sufficient to uniquely constrain the calibration. To help alleviate this problem, base Tfield construction used soft data to provide additional constraints that could not be incorporated directly into the calibration process. Specification of soft data was used to create physically realistic starting points for the calibration. The starting point for the calibration has the most impact at locations distant to the WIPP LWB.
Kriging is a linear estimation process in the field of geostatistics that predicts an average value at locations without observations, using available observations and a model describing the variability of the function (i.e., the variogram, which is itself estimated from data). Indicator kriging is a specific form of the kriging where cutoffs are estimated (i.e., is the value above or below 1.0?), rather than a continuous value. Conditional stochastic simulation is a geostatistical approach for generating realizations that will have a common specified statistical structure specified through a variogram and data, but are otherwise random. Kriging predicts the mean and variance of a field, resulting in smooth mean fields. Kriging would be conceptually similar to generating many stochastic simulations and averaging the results. Conditional stochastic simulation with indicator kriging was used to predict location of high and lowT areas
(is log_{10} T > −5.4 or < −5.4?), taking the model indicator variogram and various hard and soft data into account.
The constraints used to construct the base Tfields included a classbased linear regression relationship between log_{10} T and Culebra overburden within each type of well (see Section TFIELD4.1) and geologic soft data such as the presence of halite in nearby units or gypsum cements in the Culebra (see Section TFIELD4.2). The indicator variograms were constructed from these data (see Section TFIELD4.3) and used to stochastically simulate the cutoff between high and low Culebra T (see Section TFIELD4.4). The indicator kriging simulation result was a component of the base Tfield construction (see Section TFIELD4.5).
The best fit to estimate T from single well tests was based on a multiline regression analysis. The wells were separated into three groups: wells in the Salado dissolution zone, wells with lowT pumping test results, and wells with highT pumping test results. Figure TFIELD 41 shows the log_{10} T values from pumping test results, the Culebra overburden thickness, and the regression lines fit to each group's data individually. The cutoff between low and high log_{10} T is −5.4. Wells located where the Culebra is bounded above and below by halite (SNL6 and SNL15) were considered outliers and were not included in the regression analysis. Instead, the β _{5} I_{H} term was chosen to yield values close to those interpreted from tests at SNL6 and SNL15 (presented in Appendix F of Hart et al. (Hart et al. 2008), Table F1); this value was directly modified during the calibration stage in AP114 Task 7 (Hart et al. 2009). The final regression equation for Y = log_{10} T (Equation TFIELD 4.1) and a table of the β values (Table TFIELD 41) resulted in a fit characterized by R ^{2} = 0.92 and F = 216.
The remainder (ε) represents the misfit between the regression model and observed data.
Table TFIELD 41. βvalues for Regression Equation TFIELD 4.1
β _{1} 
β _{2} 
β _{3} 
β _{4} 
β _{5} 
−5.69805 
−3.48357×10^{−3} 
2.06581 
0.68589 
−4.75095 
The data and calculations were provided in Appendix A of Hart et al. (Hart et al. 2008).
Figure TFIELD 41. Regression Lines for LowT Wells (Blue), HighT and Nondissolution Wells (Green), and Wells Within the Salado Dissolution Zone (Red). Open diamonds are wells new to the CRA2009 PABC regression analysis (i.e., not included in CRA2004 PABC).
Geologic and hydraulic information are included as soft data to maintain the geologic conceptual model through the stochastic indicator kriging simulations in Section TFIELD4.4. Soft data define probabilities (P _{low}) a new well at a given point would have a low T value. For model cells that include wells where log_{10} T (m^{2}/s) has been estimated from singlewell hydraulic tests, the observation is referred to as hard data to distinguish it from more indirect contributions to T values associated with soft data. Model cells where hard data (singlewell testderived log_{10} T) is greater than −5.4 are assigned P _{low} = 0, while P _{low} = 1 for all cells containing lowT pumping test results. Estimated T used as hard data are presented in Table TFIELD 42, including coordinates, depth, and log_{10} T values used in regression model (from Listing A.1 of Appendix A in Hart et al. (Hart et al. 2008)).
Table TFIELD 42. Listing of Coordinates, Culebra Depth, and Log_{10} T Estimates from Singlewell Tests (Hard Data) Used in Regression Model (Equation TFIELD 4.1)
Well 
UTM X NAD27, 
UTM Y NAD27, 
depth to 
log_{10}
T

H10b 
622975 
3572473 
419.25 
−7.4 
P15 
610624 
3578747 
129.24 
−7.0 
WIPP12 
613710 
3583524 
250.7 
−7.0 
AEC7 
621126 
3589381 
269.14 
−6.8 
H15 
615315 
3581859 
265.79 
−6.8 
WQSP3 
614686 
3583518 
260.38 
−6.8 
H12 
617023 
3575452 
254.97 
−6.7 
H5c 
616903 
3584802 
277.82 
−6.7 
WIPP30 
613721 
3589701 
195.69 
−6.7 
H17 
615718 
3577513 
219.03 
−6.6 
SNL8 
618523 
3583783 
291.5 
−6.6 
WIPP21 
613743 
3582319 
225.85 
−6.6 
WQSP6 
612605 
3580736 
180.31 
−6.6 
CB1 
613191 
3578049 
157.27 
−6.5 
H14 
612341 
3580354 
170.23 
−6.5 
SNL10 
611217 
3581777 
182.58 
−6.5 
WIPP18 
613735 
3583179 
243.08 
−6.5 
SNL13 
610394 
3577600 
118.26 
−6.4 
WIPP22 
613739 
3582653 
229.51 
−6.4 
ERDA9 
613696 
3581958 
218.08 
−6.3 
C2737 
613597 
3581401 
205.74 
−6.2 
H2c 
612666 
3581668 
192.94 
−6.2 
WIPP19 
613739 
3582782 
233.93 
−6.2 
H16 
613369 
3582212 
217.46 
−6.1 
H4c 
612406 
3578499 
153.31 
−6.1 
H1 
613423 
3581684 
209.55 
−6.0 
P17 
613926 
3577466 
173.89 
−6.0 
WQSP5 
613668 
3580353 
200.67 
−5.9 
D268 
608702 
3578877 
115.98 
−5.7 
H18 
612264 
3583166 
213.57 
−5.7 
SNL5 
611970 
3587285 
194.16 
−5.3 
H19b0 
614514 
3580716 
229.2 
−5.2 
DOE1 
615203 
3580333 
253.44 
−4.9 
WQSP4 
614728 
3580766 
236.42 
−4.9 
H3b1 
613729 
3580895 
207.87 
−4.7 
WQSP2 
613776 
3583973 
249.72 
−4.7 
WQSP1 
612561 
3583427 
215.79 
−4.5 
H6c 
610610 
3584983 
187.61 
−4.4 
SNL9 
608705 
3582238 
167.64 
−4.4 
Engle 
614953 
3567454 
204.22 
−4.3 
H11b4 
615301 
3579131 
223.93 
−4.3 
SNL14 
614973 
3577643 
198.12 
−4.3 
WIPP13 
612644 
3584247 
217.17 
−4.1 
DOE2 
613683 
3585294 
254.51 
−4.0 
H9c 
613974 
3568234 
201.78 
−4.0 
SNL18 
613606 
3591536 
163.98 
−3.9 
SNL2 
609113 
3586529 
138.99 
−3.8 
WIPP25 
606385 
3584028 
140.06 
−3.6 
WIPP28 
611266 
3594680 
131.98 
−3.6 
P14 
609084 
3581976 
178 
−3.5 
SNL17 
609863 
3576016 
101.19 
−3.5 
SNL19 
607816 
3588931 
103.94 
−3.4 
WIPP11 
613791 
3586475 
256.95 
−3.4 
SNL12 
613210 
3572728 
166.73 
−3.3 
USGS1 
606462 
3569459 
162.44 
−3.3 
WIPP27 
604426 
3593079 
92.97 
−3.3 
SNL1 
613781 
3594299 
181.66 
−3.2 
SNL3 
616103 
3589047 
229.51 
−3.0 
WIPP29 
596981 
3578694 
8.23 
−3.0 
SNL16 
605265 
3579037 
58.83 
−2.9 
WIPP26 
604014 
3581162 
60.2 
−2.9 
H7c 
608095 
3574640 
77.88 
−2.8 
Two geologic margins, M2/H2 and M3/H3, were updated by Powers (2007a and 2007b), as summarized in Section TFIELD2.2.1. Wells penetrating the Culebra in areas that are bounded both above and below by halite (e.g., SNL6 and SNL15) have been found to have very low T estimates, less than 10^{−11} m^{2}/s (Roberts 2007). Wells bounded by only one margin (e.g., H12 and H17) have lower than average T estimates.
Because highT fractures are not predicted where halite is present in the Rustler, model cells located on the combined M2/H2 and M3/H3 margin were assigned P _{low} = 1. This ensured that no highT areas were placed on the boundary itself, largely a cosmetic consistency fix. Additionally, regression results for all model cells in the halite zone were replaced with values directly from the regression equation; indicator values were only used in Zones 3 and 4, and were not used east of the Rustler halite margins in Zone 1 (Figure TFIELD 32).
In all cases where sufficient core and T estimates exist, wells with no gypsum (Figure TFIELD 212) have high T, due to wellinterconnected fractures. To account for this relationship, cells were assigned P _{low} = 0.05 where no gypsum is present. As seen in Figure TFIELD 212, this is a fairly large area. Rather than give all the cells in the area such a low P _{low} value, cells were selected from a regular grid at 1300m spacing to receive soft data assignments (Figure TFIELD 43). A grid of 1300m spacing was chosen to provide sufficient definition of the boundaries without overwhelming the SISIM geostatistical simulation with too many samples. The size of the matrix decomposed by SISM during estimation is proportional to the number of samples considered at each estimation location.
It was observed that in all cases where sufficient core and T estimates exist, wells outside of the lowgypsum region (Figure TFIELD 43) have low T because fracture interconnectivity is limited by gypsum cements. In the indicator kriging, areas outside of the lowgypsum region were assigned P _{low} = 0.95 to increase the likelihood of predicting low T in the simulation.
By definition, the areas of nogypsum and highgypsum content cannot overlap, therefore the highgypsum data were sampled on the same grid used by the nogypsum data. By using fractional likelihoods and sparse sampling, these soft data did not overwhelm the random sampling algorithm of SISIM and allowed for greater variation between base field realizations. The high/lowgypsum content map is shown in Figure TFIELD 213. The lowgypsum region was not sampled, since it overlapped the nogypsum region. Instead, the highgypsum region was used. The area of high gypsum directly north of the WIPP LWB was sampled at 300m spacing, to compensate for the diffusivity soft data described in the next section and produce model results more similar to observed Culebra behavior during pumping tests.
Figure TFIELD 42. Diffusivity Values Calculated Between Wells From Pumping Test Data. Connections where log_{10} D > 0.2 are included as conditioning data with P _{low} = 0.25.
Figure TFIELD 43. Soft Data Points (Open Symbols) Generated During Step 2. Hard data points (filled symbols) are located at wells with singlewell estimates of T. The black square is the WIPP LWB.
Soft data were used to incorporate the degree of hydraulic connection observed between pairs of wells (pumping and observation wells) into the construction of the base fields. The diffusivity D (m^{2}/s) associated with the pumping and observation well pair was calculated from the results of many hydraulic tests conducted at the WIPP site (Beauheim 2007). A map of these values is shown in Figure TFIELD 42, showing colored lines connecting pumping and observation wells involved in the nine pumping tests used in Tfield calibration (pumping/observation wells listed in Table TFIELD 52). The model cells falling on a straight line connecting two wells with a calculated log_{10} D > 0.2 (i.e., all red connecting lines and some orange dashed lines) were assigned P _{low} = 0.25 to account for the increased likelihood a cell on the connecting line would be high T (inverted maroon triangles in Figure TFIELD 43). Using P _{low} = 0 would have forced SISIM to create a direct path connecting two wells where a strong response to pumping was observed, and there is no geologic reason that these connections must be straight.
In addition to the highT connection lines, a set of lowT points was placed roughly parallel to the SNL14/SNL12/H9/Engle connection path to keep the highT connection relatively narrow (blue circles in south central portion of Figure TFIELD 43). These points were assigned a P _{low} = 1, to ensure they would impact the indicator kriging simulation. Pumping SNL14 in 2005 produced a strong response at H9c nearly ten kilometers (km) to the south. During model development it was found the only way to recreate the observed response with the MODFLOW model was to incorporate a relatively narrow connecting zone of high T. Without adding some lowT points along the flanks of this path, SISIM tended to create a wide highT area, which did not allow the drawdown response to propagate significantly from SNL14 to H9, as was observed. The observed response would be consistent with a narrow linear geological feature, which is difficult to simulate using the current MODFLOW model with 100m grid spacing. These lowT points did not force the simulation to create a narrow highT pathway in each realization, as many base fields still had large areas of high T that extend past these points. Fields generated with and without the narrow highT pathway were modified through the calibration process, which included the SNL14 pumping test data as calibration targets. This exercise was performed in an attempt to improve the ability of the base T fields to match observed data, since this might lead to fewer PEST iterations and quicker calibration times.
The final combined soft data field is shown in Figure TFIELD 43. The soft data input files and calculation scripts are provided in Appendix B of Hart et al. (Hart et al. 2008).
Singlewell estimates of T are hard data shown on the figure with filled diamonds (data listed in Table TFIELD 42). Hard data are combined here with soft data in base Tfield creation, but appear again (without soft data) as fixed pilot points in the Tfield calibration process (Section TFIELD5.0). Filled green diamonds are wells with log_{10} T estimates ≤ −5.4 (m^{2}/s), and black filled diamonds are wells with log_{10} T estimates > −5.4 (m^{2}/s).
The grid of inverted open blue triangles in the east indicate areas with "high gypsum" (white area in Figure TFIELD 213), while the grid of open red triangles in the west indicates areas with "no gypsum" (peach area in Figure TFIELD 212).
Lines of closely spaced inverted brown triangles represent the connections between pumping and observation wells interpreted with a log_{10} diffusivity (D) > 0.2 m^{2}/s, including all solid red and some dashed orange lines in Figure TFIELD 42.
The open blue circles ("Halite" in the figure legend) are used in two ways to enforce high probability for low T in two different locations. The first use is the line of closely spaced open blue circles corresponding to M2/H2 or M3/H3; locations east of this line have either halite above or below the Culebra. This boundary marks the eastern edge of the random placement of highT and lowT zones (Zones 3 and 4 in Figure TFIELD 32). The second group, the open blue circles straddling the line connecting Engle, H9, SNL12, and SNL14 (running south to north from the bottom middle of the figure), represents a lowT zone added to increase the contrast of the highT zone along this line of southtonorth connectors, to better represent results observed in the SNL14 multiwell pumping test.
The geostatistical indicator simulations done as part of the base Tfield development are only utilized in the central section of the model domain, between the Salado dissolution area to the west and the lowT halitesandwiched region to the east. Therefore, only wells in this middle section are used for construction of the indicator variogram. A total of 46 wells that provide information regarding log_{10}
T were used in the calculation of the indicator variograms. The indicator value is determined by comparing each log_{10}
T value to a threshold log_{10}
T value,
T_{t}
= −5.4,
Where I(x,y) denotes the unitless indicator value at well location (x,y). The experimental indicator variogram was fit with a spherical variogram model, whose model parameters are given in Table TFIELD 43. Figure TFIELD 44 illustrates the experimental and model indicator variograms. The proportion of lowT values in the data set is 0.652. The variance of an indicator value is (1 − p)p, where p is the proportion of high or low values. The variance for these indicator data is 0.227 and is used directly as the sill in the variogram modeling (dashed horizontal line in Figure TFIELD 44). The parameters in Table TFIELD 43 are used as input to the SISIM program for creation of the stochastic component of the base Tfields. In an attempt to identify anisotropy, model variograms were calculated in both the NESW and NWSE directions (see Appendix C of Hart et al. (Hart et al. 2008)). These directional variograms analyses were inconclusive, only omnidirectional (i.e., isotropic) variograms were used in the final analysis.
Table TFIELD 43. Variogram Parameters for Isotropic Fit to Indicator Data Variogram. Omnidirectional variogram calculated with a lag spacing of 500 m.
Parameter 
Value 
Model Type 
Spherical 
Nugget 
0.0 
Sill 
0.227 
Range 
2195 m 
Figure TFIELD 44. Experimental Variogram (Dots) and Spherical Model (Line) for Indicator Values. xaxis is lag distance [meters], yaxis is the unitless indicator; numbers by dots indicate the number of pairs represented at each lag.
With previous sections describing the indicator variogram model, hard T data values, and soft geologic and hydrologic data, stochastic realizations of highT zones were constructed using the GSLIB program SISIM (Deutsch and Journel 1998). An example indicator field is given in Figure TFIELD 45. Maps summarizing statistics for the 1,000 resulting base Tfields are presented in Figure TFIELD 46 and Figure TFIELD 47. These figures show the impact the conditioning information had on the overall fields. The combined M2/H2 and M3/H3 margins have a standard deviation of 0 and are constant at the proper value as desired. Areas designated as higher likelihood of high T do show an average value that trends towards the highT value (in this case, 0), but they still have a standard deviation that is nonzero, indicating that there is some variability in those areas. The same is true in areas outside the lowgypsum region. Additionally, areas with no conditioning information have higher standard deviations, indicating that highT zone placement in those locations was allowed to be fully variable. Though there are some visible artifacts from the grids used in the average and standard deviation fields (locations of soft data points in Figure TFIELD 43 are discernible in Figure TFIELD 46 and Figure TFIELD 47), the individual realizations, such as Figure TFIELD 45, do not show these artifacts. Additionally, the majority of the artifacts occur outside the central zone, which is the only place the indicator fields are used. The indicator fields created by this process are the best possible combination of hydraulic and geologic conditioning given current data.
There is a relatively high density of hard data available inside the WIPP LWB, which significantly constrains the estimation process there. Geostatistical simulation is most useful to fill in large gaps near the edges of the model domain where a small number of observations exist. It must also be considered that these baseT fields are just a first guess for the model calibration process, which utilizes the singlewell T and both the steadystate and multiwell pumping test drawdown data as calibration targets. Ultimately, these data drive the calibration, adjusting input parameters to better match observed data. See Section TFIELD5.4 for figures and discussion illustrating the extent to which the input fields were adjusted to match the data (e.g., see Figure TFIELD 519).
Figure TFIELD 45. Sample Indicator Field for Realization r123. 1 indicates low T and 0 indicates high T.
Figure TFIELD 46. Average Indicator Values Across All 1000 Base Realizations. 1 indicates low T and 0 indicates high T.
Figure TFIELD 47. Standard Deviation of Indicator Values Across All 1000 Base Realizations
Once the indicator fields were created, the T values were assigned using Equation TFIELD 4.1 using a Perl script. Equation TFIELD 4.3 was used to calculate Y = log_{10} T at each cell,
Y(x,y) = b [Z(x,y)] + a [Z(x,y)] d(x,y) (TFIELD 4.3)
where b and a represent combinations of the βcoefficient based on the zone (Z(x,y)) of the cell. Table TFIELD 44 shows how the variables in the original linear regression equation (Equation TFIELD 4.1) were related to Equation TFIELD 4.3. Figure TFIELD 32 shows the indicator zone distribution.
Table TFIELD 44. Correlation of β and I Values from Equation TFIELD 4.1 to a and b Values in Equation TFIELD 4.3

Zone 0 Salado 
Zone 1 Halite 2 
Zone 2 Halite 
Zone 3 Central low T 
Zone 4 Central high T 
I_{f} 
1 
0 
0 
0 
1 
I_{D} 
1 
0 
0 
0 
0 
I_{H} 
0 
1 
0 
0 
0 
I_{h} 
0 
1 
1 
0 
0 
b 
β _{1} + β _{3} + β _{4} 
β _{1} + β _{5} 
β _{1} 
β _{1} 
β _{1} + β _{3} 
a 
β _{2} 
β _{2} 
β _{2} 
β _{2} 
β _{2} 
The Perl script was executed on all 1,000 realizations. A sample final base T field is presented in Figure TFIELD 48 for realization r123 (a random representative selection from the possible fields). The mean log_{10} Tfield across all 1,000 realizations is presented in Figure TFIELD 49. The standard deviation of log_{10} T is presented in Figure TFIELD 410. Very low standard deviation occurs across the 1000 base realizations outside the stochastically sampled areas, including the higher T areas to the west and the lower T areas to the east (Figure TFIELD 410). These stochastically sampled areas are the main source of variability between the 1000 base realizations. The variability of the indicator variable across the realizations (Figure TFIELD 47) is one component of the variability observed in the final T values in the base T fields (plotted normalized to the range {0,1} in Figure TFIELD 410). The regression analysis produced variability in the predicted T values in a zone related to the variability in the overburden thickness.
Figure TFIELD 48. Sample Log_{10} T (m^{2}/s) Base Field Realization r123
Figure TFIELD 49. Mean Log_{10} T (m^{2}/s) Values Across All 1000 Base Realizations
Figure TFIELD 410. Normalized Standard Deviation of Log_{10} T (m^{2}/s) Values Across All 1000 Base Realizations
The work outlined in Section TFIELD5.0 was performed under AP114, Analysis Plan for Evaluation and Recalibration of Culebra T Fields (Beauheim 2008). The calibration of the Tfields used 200 of the 1,000 base fields from the results of AP114 Task 5 (Hart et al. 2008, summarized in Section TFIELD4.0) as starting points for the calibration process. More than 200 fields could not be calibrated, due to time constraints. Calibration is the process of systematically adjusting the input parameters to the MODFLOW model (fields of T, A, S, and R) to reduce the sum of the squared differences between field observations and MODFLOW model output (steadystate and transient head).
The pilot point calibration method was implemented using the parameter estimation software PEST. Automatic model calibration was utilized to make the process more easily documentable and reproducible, compared to manual calibration (i.e., trial and error). The MODFLOW model used to simulate groundwater flow through the Culebra contains a large number of active model cells for T, A, and S fields. Estimating each model element independently would require estimating hundreds of thousands of unknown parameters. The pilot point approach makes the calibration process more tractable by lowering the number of parameters to estimate. Instead of estimating each parameter in each model cell independently, parameter values are estimated at strategically placed pilot point locations. Parameter values at each model cell are then interpolated from the pilot points using kriging as a preprocessing step between parameter assignment by PEST and MODFLOW model execution. The pilot point approach allows mixing estimated and fixed pilot points (e.g., T pilot points at wells with singlewell hydraulic test estimates of T). The pilot point approach was also used in CRA2004 PABC and a variant of it was used (without PEST) in the CCA. Both steadystate and transient head calibration targets are discussed in Section TFIELD5.1. The parameter zones, pilot point locations for each parameter, initial conditions, and boundary conditions are specified in Section TFIELD5.2. The components of the MODFLOW model used with PEST, and the utilities required to preprocess and postprocess the inputs and outputs from the model during the calibration are discussed in Section TFIELD5.3. Finally, some postcalibration analysis of the results is presented in Section TFIELD5.4.
The initial T values at pilot points were taken from the base fields. In addition to T, the horizontal T anisotropy (A), the storativity (S), and a linear section of recharge (R) were also calibrated. The same zone definitions used for developing the base T fields (see TFIELD4.0) were used for T pilot points (although zone numbers changed), and similar zone definitions were used for anisotropy. Zones for storativity and recharge were based on other analyses completed in the area surrounding the WIPP site (see Section TFIELD5.2.1). Pilot points were selected for each parameter and initial values were selected that were consistent with the conceptual model used to create the base fields (see Section TFIELD5.2.2 and Section TFIELD5.2.7).
A model variogram for T was created using the estimated T from singlewell hydraulic tests. This variogram was also used for all parameters, as it was the only one that could be created from field data (see Section TFIELD5.2.8). This variogram was used to create kriging factors that were then used to create continuous fields from the pilot point values. The T, A, S, and R fields were then used as inputs to the MODFLOW numerical model to produce simulated head and drawdown results (see Section TFIELD5.2.10).
Once the MODFLOW models produced simulated drawdown and head results, the modeled results were compared to the field data for the tests that were modeled. The residual of an observation is calculated in PEST as the weighted difference between measured and modeled data. Observation weights were selected to make the sum of the weighted steadystate head errors approximately equal to the sum of errors of four observation wells in a transient pumping test to approximately balance the steadystate and transient modeltodata misfits. The PEST optimization uses a single objective function, which is the sum of the steadystate and transient residuals. Because the improvement of model fit for steadystate heads might come at the expense of fit to transient pumping tests, a decision was made to balance the importance of the two groups in the calibration effort. The residuals were used by PEST to construct a finitedifference approximation of the Jacobian matrix. The Jacobian matrix is a measure of the sensitivity each model prediction has to each adjustable parameter, and is used to optimize the pilot point parameter values. The goal of the optimization is to minimize the objective function value, a measure of the misfit between model predictions and observed data (see Section TFIELD5.2.9).
Because traditional construction of the Jacobian matrix requires at least N_{p} + 1 forward model calls (N_{p} parameters estimated in the calibration), using 1,100plus parameters would be impossible without additional efficiency in the optimization. The PEST singular value decomposition (SVD) assist approach reduces the size of the Jacobian matrix by only using the most significant "superparameters" that correspond to the eigenvectors with the largest singular values, estimated using the SVD of the Jacobian matrix. The SVD process required initial calculation of a full Jacobian matrix, but then reduced the subsequent number of required forward calls by a factor of four to six. The result was that three calls to PEST were required to calibrate the fields (see Figure TFIELD 51):
1. a single full Jacobian calculation, which required 1,100+ forward model calls;
2. an SVD calibration using the reduced parameter set that ran up to 50 iterations, requiring between 100 and 400 forward model calls per iteration; and
3. a final PEST run with the best parameter results to create the final fields corresponding to the best parameter values calculated during the SVDassisted calibration.
Total calibration time for a single base field was approximately seven days using six processors (one master node and five slave nodes).
Figure TFIELD 51. Complete Calibration Process for a Single Realization
After approximately 140 fields had been calibrated, a few steadystate calibration targets were found to be incorrect by several meters. A total of 150 fields were calibrated using the incorrect targets, and an additional 50 fields were started using the corrected heads (Beauheim 2009; Johnson 2009a and Johnson 2009b). To deal with the incorrect values, a limited recalibration was performed on the results of the first 150 calibrations (see Section TFIELD5.3.2). The same process that has been described was followed for the limited secondary calibration, but since the initial parameter values were taken from the calibrated results, only the necessary pilot point locations near the updated steadystate head values were allowed to be changed, and the SVD portion of the PEST recalibration was limited to 10 iterations. The end result was 200 fields, with 150 of these fields having undergone a secondary calibration to incorporate corrected field observation data. The impact of this change is discussed in Section TFIELD5.3.3.
The end result of the calibration was the first 200 of the original 1,000 field sets (T, A, S, and R) calibrated to one set of steadystate heads and nine transient, multiobservation well pumping tests. More than 200 of the 1,000 base Tfields could not be calibrated due to time constraints associated with the calibration effort. Since the 1,000 base Tfields were generated randomly, the first 200 were in effect a random sample. The 100 bestcalibrated fields (those with the smallest residuals) were selected as the final results from this task (see Section TFIELD5.3.4). Several statistical analyses were performed on the fields themselves to generate average values and variances in the field results (Section TFIELD5.4.1).
The complete calibrations were performed under qualityassurance run control with inputs and outputs stored in a version control repository on a central server. The calibrations required approximately six months of continuous runtime on a total of 80 processors. See Hart et al. (Hart et al. 2009) for details.
Two sets of head values were used for calibration of the Culebra flow model. The first dataset consists of 42 freshwater head values measured in May 2007, which were used as steadystate calibration targets (see Table TFIELD 51 for values used and see Johnson (Johnson 2009a and Johnson 2009b) for details). Appendix HYDRO2014 shows the behavior of water levels in wells since 2007. The data in Appendix HYDRO2014 support the continued use of 2007 water levels to represent steadystate conditions at the WIPP up to the CRA2014 (see additional discussion in Appendix HYDRO2014). The second dataset consists of drawdown data collected during 9 independent multiwell pumping tests conducted over a span of more than 20 years (see Table TFIELD 52 for lists of observation wells and relevant references).
The steadystate MODFLOW simulation is calibrated against data in Table TFIELD 51. Each of the nine transient pumping tests in Table TFIELD 52 are simulated as a separate MODFLOW model simulation. A "single" forward model run actually involved 10 individual MODFLOW simulations (1 steadystate and 9 transient simulations), each using the same input parameter fields, but different pumping configurations. Some wells appeared in multiple pumping tests (e.g., H15, H18, or ERDA9). These wells had calibration targets associated with multiple transient forward simulations. The MODFLOW model is discussed more in Section TFIELD5.2.10.
SNL6 and SNL15 are listed as steadystate targets in Table TFIELD 51, but they are located in the constanthead portion of the model domain and therefore their corresponding modelpredicted values are not adjustable through changes in the T field.
Table TFIELD 51. Freshwater Head Observations Used as Steadystate Calibration Targets
Well Name 
May 2007 
C2737 
921.23 
ERDA9 
924.88 
H2b2 
929.62 
H3b2 
918.68 
H4b 
916.34 
H5b 
939.12 
H6b 
936.44 
H7b1 
914.58 
H9c 
912.8 
H10c 
922.02 
H11b4 
917.09 
H12 
916.53 
H15 
920.32 
H17 
916.24 
H19b0 
918.84 
IMC461 
928.95 
SNL1 
941.86 
SNL2 
937.65 
SNL3 
939.81 
SNL5 
938.59 
SNL6 
1110 
SNL8 
929.94 
SNL9 
932.05 
SNL10 
931.54 
SNL12 
915.24 
SNL13 
918.19 
SNL14 
916.33 
SNL15 
1060 
SNL16 
918.68 
SNL17A 
916.78 
SNL18 
939.87 
SNL19 
937.58 
USGS4 
911.11 
WIPP11 
940.65 
WIPP13 
939.78 
WIPP19 
933.66 
WIPP25 
937.57 
WIPP30 
939.37 
WQSP1 
938.28 
WQSP2 
939.87 
WQSP3 
936.43 
WQSP4 
919.5 
WQSP5 
918.18 
WQSP6 
921.96 
Table TFIELD 52. Summary of Transient Observations Used as Calibration Targets
Pumping Well(s) 
Observation Wells^{1} 
Total # obs. 
Reference 
WQSP2 
H18, DOE2, WQSP1, WIPP13 
77 
Beauheim and Ruskauff 1998 
H19 and H11 
WQSP5, H1, H15, DOE1, ERDA9, WIPP21, H3b2 
143 
Beauheim and Ruskauff 1998 
WQSP1 
H18, WIPP13 
36 
Beauheim and Ruskauff 1998 
WIPP11 
WQSP2, WQSP3, WIPP12, SNL9, SNL5, H6b, SNL3, SNL2, SNL1, WIPP30, WQSP1, WIPP13 
250 
Toll and Johnson 2006b; Roberts 2006 
H11 
H4b, H12, H14, H15, H17, DOE1, CB1, P15, P17, H3b2 
130 
Beauheim 1989 
WIPP13 
WIPP19, WIPP18, H2b2, H6b, WIPP12, WIPP25, DOE2, WIPP30, P14 
167 
Beauheim 1987b 
SNL14 
H9c, H4b, SNL13, SNL12, H15, 
252 
Toll and Johnson 2006a; Roberts 2006 
P14 
D268, H18, WIPP26, H6b, WIPP25 
82 
Beauheim and Ruskauff 1998 
H3 
DOE1, H2b2, H1, H11b2 
69 
Beauheim 1987a 
1. ^{ 1 } See Figure TFIELD 42 for locations of pumping and observation wells and diffusivity values estimated from pumping tests. 
This step comprised the setup and configuration processes that were the same for every calibrated base field. Step 1 included the creation and definition of zones for each of the parameters and the selection of pilot point locations and initial values. Because of the stochastic nature of the transmissivity fields, unique zones are associated with each field, as defined in Hart et al. (Hart et al. 2008). The process to set up each field was the same, but certain elements, such as the exact number of pilot point locations used, were unique to each field.
The model domain and extent are identical to the domain defined in Section TFIELD3.1. The base Tfields were taken from the results of Section TFIELD4.5. Some model input files were created statically, and were used for every calibration. Some model input files were created dynamically for each calibration, but used the same variables and parameters in their creation.
Parameter zones were defined for each of the calibrated parameters. These zones or regions were defined to be consistent with the conceptual model for Culebra flow defined in Hart et al. (Hart et al. 2008) and summarized in Section TFIELD3.0. The parameter zones were chosen to organize common pilot points groups together in the PEST calibration. Numerical values of parameter zones (i.e., zone numbers) are often different from those used in the conceptual model. Figure TFIELD 28 shows the different margins that define geologic zones and the locations of the Culebra wells that have been drilled in the vicinity of the WIPP.
The T zones were defined to be consistent with the zones used in the geologic model and soft data analysis (Section TFIELD4.2). As shown in Figure TFIELD 52, a highT zone exists to the west (zone 2), corresponding to the area of Salado dissolution, and a mixture of higher and lower T values corresponding to the stochastic zones (zone 0 and 1) provided in the base fields defined the center of the model domain. Unlike AP114 Task 5 (Section TFIELD4.5), where it was a separate zone, the area between the H2/M2 and H3/M3 margins (green in Figure TFIELD 52) was included in the same zone (zone 1) as the lower T stochastic areas provided by the base fields. The area east of both the H2/M2 and H3/M3 margins  where the Culebra is bounded both above and below by halitecemented elements  was defined to be its own zone (zone 3), as was done in AP114 Task 5. A noflow boundary roughly follows the axis of Nash Draw defining the final region (zone 4), which is inactive and applies to all parameters.
Figure TFIELD 52. Transmissivity Zone Map for a Single Realization. Zones 0 and 1 are the stochastic zones created in Hart et al. (2008); Zone 2 is the highT Salado dissolution area; Zone 3 is the very lowT halitebounded area; Zone 4 is the northwestern inactive area.
The T anisotropy field used the same zones defined for T. The T values in the northsouth (y) direction were calculated by multiplying the transmissivity value for the eastwest (x) direction (given in the Tfield) by the anisotropy value (A) at a given cell
T _{NS} = T _{EW} A. (TFIELD 5.1)
The map shown in Figure TFIELD 52 represents the zonation used for both A and T.
Besides the noflow area (zone 4), four zones were defined for storativity estimation. The westernmost region (zone 2) is unconfined, as described in Powers et al. (Powers et al. 2006) (summarized in Section TFIELD2.3.2). The largest area (zone 0) is fully confined, with its western boundary roughly following the Salado dissolution margin. The area between these two regions (zone 1) is the transition zone, where the Culebra is uncertainly confined. As with the T and A zones, the area east of both the H2/M2 and H3/M3 margins is a separate region (zone 3), but in this case storativity is simply held constant at the initial confinedzone value. A map of the storativity zones is shown in Figure TFIELD 53.
The unconfined zone was implemented as a zone of high confined storativity to simplify the numerical model by approximating the nonlinear unconfined problem with a linear storativity one. By defining a much higher storativity value in the unconfined part of the domain, unconfined behavior can be approximately modeled using a confined numerical model, which is linear and executes quicker than the unconfined nonlinear model. Since the unconfined or transition zone does not exist inside the WIPP LWB, this choice has little impact on interpretation of transient hydraulic tests there, and this choice has no impact on steadystate flow simulations (S is only used in transient simulations) used to predict radionuclide transport in PA.
Figure TFIELD 53. Storativity Zones. Zone 0 is confined; Zone 1 is a transition between confined and unconfined; Zone 2 is unconfined; Zone 3 is confined and uncalibrated from the initial confined value; Zone 4 is inactive (no flow).
The conceptual model presented in Holt (Holt 1997) and Hart et al. (Hart et al. 2008) indicates that a groundwater divide exists somewhere southwest of the WIPP site. Previously in CRA2004 PABC and CRA2009, this groundwater divide was represented by extending the noflow zone all the way to the southern boundary (McKenna and Hart 2003). Because the model used in this current task included the unconfined zone, it was decided to model the groundwater divide using recharge instead of a noflow boundary. The areas of possible recharge were defined in AP114 Task 1B (Powers 2006), summarized in Section TFIELD2.3.1. Recharge values had to be extremely small (on the order of 10^{−11} m/s) to ensure convergence of steadystate MODFLOW simulations. Rather than try to determine which of the configurations presented in Task 1B was the "best" approximation, a similar simple approximation to the older noflow approach was used. A recharge zone consisting of a line of cells extending NW to SE along the axis of the largest topographic feature (and roughly following the old noflow boundary from McKenna and Hart (McKenna and Hart 2003) was used as the recharge zone. See Figure TFIELD 54 for a map of the recharge zone.
Figure TFIELD 54. Zone 1, the Zone of Nonzero Recharge. Zone 1 is the linear feature directed southeast from the center of the west edge of the model domain. The remaining model area has no recharge.
While the boundary conditions were not variable parameters in the calibration, the definition of the specifiedhead boundary conditions was an important part of the setup and configuration step. The noflow zone in the northwest was defined to be the same as was used in AP088 for CRA2004 PABC. Though the T in this area is extremely low (10^{−13} to 10^{−11} m^{2}/s), there should be some flow exiting along the zone margin, however minute. Testing at SNL6 and SNL15 indicates that the hydraulic heads in this area may be recovering ultimately to levels above the land surface (Roberts 2007; Bowman and Roberts 2008; Appendix HYDRO2014, Section 7.1 ). The "halitesandwiched" zone east of either M2/H2 or M3/H3 was simultaneously made extremely low T and set to fixedhead values at the ground surface elevation. While this meant that the head values at SNL6 and SNL15 were no longer estimable, it was considered the simplest way to model the nearly stationary nature of the water in this zone using MODFLOW2000. The flow zones are shown in Figure TFIELD 55, and the selection of the initial values is discussed in Section TFIELD5.2.7. The northern, western, and southern flow boundaries were all fixedhead boundaries.
Figure TFIELD 55. Flow Zones. The fixedhead zone is green; the noflow zone is salmon; the white area is normal flow. The fixedhead zone includes one cell along the northern, southern, and western boundaries.
Once the zones were defined, pilot point locations were selected. There were two types of pilot points, fixed and variable, and two placement approaches, gridded and linear. Selection of the points for each parameter required a combination of both types and approaches. The exact algorithm used to calculate placement is detailed in Hart et al. (Hart et al.2009), and resulted in the pilot point locations used and shown in Figure TFIELD 56 through Figure TFIELD 59.
In addition to pumping test wells, extra pilot points were placed in the transmissivity fields. These were included along the northern and southern boundaries to try to limit the effects that the fixedhead boundaries would have on transient pumping and the steadystate model results. The estimated T values from singlewell pumping and slug tests were used as fixed T points that corresponded to the tested wells (see Table TFIELD 42 for testderived transmissivity values). See Table TFIELD 53 for the ranges of pilot point values used, and see Figure TFIELD 56 for pilot point locations.
Figure TFIELD 56. Transmissivity Pilot Point Locations. Blue points are fixed values and red points are variable parameters. Zones correspond to a single realization.
Table TFIELD 53. Initial Values of Parameters at Pilot Points
Parameter 
Zone 
log_{10} Value^{1} 
Pilot Point log_{10} Calibration Limits 
Transmissivity 
Zone 0 
−0.003484 d(x,y) − 3.6322 
[−19.0,−1.0] 
Zone 1 
−0.003484 d(x,y) − 5.6981 
[−19.0,−1.0] 

Zone 2 
−0.003484 d(x,y) − 2.9463 
[−19.0,−1.0] 

Zone 3 
−0.003484 d(x,y) − 10.4490 
[−19.0,−1.0] 

Anisotropy 
All Zones 
0.0 
[ −0.5, 0.5] 
Storativity 
Zone 0 
−5.0 
[ −5.5,−4.5] 
Zone 1 
−4.0 
[ −6.0,−0.5] 

Zone 2 
−1.5 
[ −2.5,−0.5] 

Zone 3 
−5.0 
Fixed 

Recharge 
Zone 1 
−11.0 
[−19.0,−1.0] 
^{1} d(x,y) is Culebra overburden thickness. 
Anisotropy was unique because no fixed values and therefore no fixed pilot points were used. This result is due to the singlewell tests not providing any estimate of anisotropy, and the multiwell tests providing too localized an estimate of anisotropy (only valid for a single cell or two in the model). See Figure TFIELD 57 for the locations of anisotropy pilot points.
The only variable storativity pilot points in the confined zone were along straight lines connecting pumping and observation wells in transient pumping tests. The gridded points were set as fixed values, since it was assumed there was no information allowing for effective calculation of storativity outside the transient tests. All pilot points located within the unconfined and transition zones were defined as variable. See Figure TFIELD 58 for the location of storativity pilot points.
Figure TFIELD 57. Anisotropy Pilot Point Locations. All anisotropy pilot points were variable parameters. Zones correspond to a single realization.
Figure TFIELD 58. Storativity Pilot Point Locations. Only pilot points along lines between wells in transient pumping tests and points in the unconfined zones (zones 1 and 2) were variable (red dots); the remaining pilot points were fixed (blue dots).
Because the recharge zone was a line, only four pilot points were needed in the entire zone. In this case, the pilot point nearest the western domain boundary was set as a fixed value of
10^{−30} m/s, which was interpreted as zero by the preprocessors to MODFLOW, and the other three were variable. See Figure TFIELD 59 for the location of the recharge pilot points.
Figure TFIELD 59. Recharge Pilot Point Locations. The pilot point along the model domain boundary was fixed, while the other three points were variable.
Section TFIELD5.2.7.1 discusses the initial values assigned to parameters before calibration, while TFIELD5.2.7.2 discusses the assignment of a head field to the initial condition and certain specifiedhead boundary conditions in the groundwater flow model.
The initial values for each of the pilot points were defined according to the conceptual model and the values presented in Hart et al. (Hart et al. 2008), and summarized in Section TFIELD3.0. For T, this meant that the same equation used to create the base T fields was used to define the initial values for the pilot points, based on their zone. Anisotropy was set to isotropic conditions (A = 1) for all points. Storativity was defined to start at 10^{−5} for the confined zone (the same value that was used for S in CRA2004 PABC model AP088), at 10^{−4} in the transition zone, and at 10^{−1.5} in the unconfined zone. Recharge was initialized as 10^{−11} m/s, a value found to be sufficiently small to allow MODFLOW to perform an initial run prior to PEST calibration. The zonebyzone initial values for each parameter, and the limits placed on the range the values could take in calibration, are presented in Table TFIELD 53. See Table TFIELD 42 for fixed T values.
Initial heads (H _{0}) were created using a multivariate equation based on normalized x and y coordinates (1 ≤ {x,y} ≤ +1). The equation was designed to keep head along the northern boundary just above the measured head at SNL1 and head along the southern boundary below the level measured at H9c. These constraints were the defining factors on the constants in the equation that follows. This process was done only once, and the result was used as a static input file for all calibrations. The field was defined by the following equations
where sign() is the sign of its argument (either +1 or 1) and y is absolute value. For values east of both the H2/M2 and H3/M3 boundaries, the groundsurface elevation was used as the initial head value; see Appendix A of Hart et al. (Hart et al. 2009) for details. The resulting initial head field is shown in Figure TFIELD 510.
Figure TFIELD 510. Initial Head Values for Use in MODFLOW Steadystate Solution. Brick red head values were fixed at the ground surface elevation (>1,000 m AMSL).
Transmissivity fields are created from pilot points using kriging. Some pilot points are adjusted using PEST, while other pilot points are held fixed, because they correspond to estimated T values at wells with pumping tests. A variogram is needed to interpolate and extrapolate from the pilot points onto every element of the MODFLOW grid.
The transmissivity variogram (different from the indicator kriging variogram discussed in Section TFIELD4.3) was created using transmissivity values estimated from well tests at 62 of the wells around the WIPP site. Wells outside the model domain and values at SNL6 and SNL15 were excluded from the calculation. The values at SNL6 and SNL15 are both several orders of magnitude lower than at the other wells, and are in a geologically distinct zone. While initial calculations showed that there was some statistical anisotropy, there were not sufficient measurements to create an anisotropic variogram. The complete steps for creating the variogram are presented in Hart et al. (Hart et al. 2009), Appendix B. The final parameters used are shown in Table TFIELD 54.
Table TFIELD 54. Parameters for T Model Variogram, Fitted to Transmissivity Data Using an Omnidirectional Variogram with Lag Spacing of 1,500 m
Parameter 
Value 
Model Type 
Exponential 
Nugget 
0.02 (log_{10} T)^{2} 
Sill 
1.95 (log_{10} T)^{2} 
Range 
9,500 meters 
The observations (steadystate freshwater heads and pumping test drawdowns) used as calibration targets for PEST are summarized in Section TFIELD5.1. Residuals are calculated as the difference between measured and modelgenerated freshwater heads or drawdowns. The PEST utility program MOD2OBS is used to extract the observations from model output at times and locations associated with each steady and transient observation.
Inverse modeling (i.e., automatic calibration) requires a numerical model which generates results to compare against observed information. In this task, a MODFLOW 2000 (Harbaugh et al. 2000) flow model was developed for the Culebra that could use the base fields generated in Hart et al. (Hart et al. 2008) as inputs. As was done in CRA2004 PABC (McKenna and Hart 2003), the link algebraic multigrid (Mehl 2001) solver was used to increase speed and performance compared to other available solvers. In addition to T, it was decided to calibrate the local horizontal T anisotropy, storativity, and a strip zone of recharge as parameters in the calibration. Having these four parameters  T, A, S, and R  required a slightly more complex MODFLOW model implementation than was used in CRA2004 PABC AP088 (McKenna and Hart 2003). Specifically, both storativity and anisotropy were single values previously, and changing these to cellbycell values required the use of the layer property flow package instead of the block centered flow package used previously. Using recharge also required the addition of the recharge package; both packages are part of the standard MODFLOW distribution (Harbaugh et al. 2000). For the known information, steadystate heads from 2007 and drawdown results from nine different pumping tests performed between 1985 and 2008 were used as the measured data. A conceptual diagram of the MODFLOW model with its inputs and outputs is shown in Figure TFIELD 511.
Figure TFIELD 511. Flow Chart Showing the Forward Model Used In the Model Calibration. T, A, S, and R are parameter fields. H represents the steadystate flow solution. DD _{1}DD _{9} represent transient drawdown computed for the 9 individual pumping tests from 9 separate forward simulations.
The calibration process used multiple forward model calls to evaluate the impact that perturbing an input parameter has on model predictions at locations and times corresponding to observations. This process was computationally intensive, and involved 80 processors on 2 different computing clusters running for 6 months to calibrate 200 of the 1,000 base T fields.
The calibration process was done using the PEST inverse modeling software suite and its groundwater utilities. The steps involved in each forward model run during the PEST calibration are illustrated in Figure TFIELD 512; the complete calibration process is shown in Figure TFIELD 51.
The completed PEST simulation included the creation of the fields from the kriging factors and pilot points (PPK2FAC, FAC2REAL, REAL2MOD), the MODFLOW calls, and finally the observation extraction utilities (MOD2OBS and OBS2REAL), which extract modeled cell head or drawdown values from a binary MODFLOW output file. For SVD iterations, another preprocessor, PARCALC, is used to create the pilot point values from the super parameters (i.e., eigenvectors related to the largest eigenvalues  see description in TFIELD5.0). The model script (model.sh), the REAL2MOD script, and the OBS2REAL script were written for this task, and are included in Appendix G of Hart et al. (Hart et al. 2009). PPK2FAC, FAC2REAL, and MOD2OBS are PEST utilities.
The first call to the PEST program was a single outer iteration to estimate the Jacobian matrix. This required over 1,100 forward model calls, one for each variable parameter value. Once the Jacobian matrix was calculated, the SVDAPREP program decomposed the Jacobian matrix into eigenvectors and kept the super parameters corresponding to the largest singular values. The result was a set of 100 to 300 super parameters that were then used with a 50iteration PEST calibration. The termination criteria were:
Once termination criteria had been reached, the PEST program would output the best parameters to a file. This file was then used to create one final PEST control file, which issued a single model run with the best parameters as input. The results of this final call were then used to calculate the measures of fit and the final fields.
The run control details regarding the calibration process are presented in Appendix G of Hart et al. (Hart et al. 2009). Using the ReadScript.py run control system allowed automatic checkout of input files, execution, and checkin of the results to version control following calibration.
Figure TFIELD 512. Flowchart Illustrating the PEST Calibration Process
Because some of the original steadystate calibration targets were incorrect, the fields that had already been calibrated to the incorrect data needed to be recalibrated to the new data. The two wells with the most significant changes, ERDA9 and SNL8, had more than one meter change from the old to new values. At ERDA9 the calibrations had consistently been unable to match the incorrect head value, which was too low compared to the higher corrected value. Without any recalibration, correcting the value for ERDA9 produced better model fits.
The new calibration target for SNL8 was based on a recalculation of the freshwater head (Johnson 2009a and Johnson 2009b). Because SNL8 was not an observation well in any of the transient pumping tests, and because it was to the east and upgradient from the WIPP LWB, only a section of the fields were recalibrated to correct for the change in the calibration target at SNL8. It was hoped that this would allow the T, S, A and R fields to change to match the SNL8 head without requiring the weeklong recalibration for each of the affected fields if the entire domain was recalibrated.
The recalibration process involved fixing all the parameters that had previously been calibrated, except for those parameters in a rectangular area around and upgradient from SNL8. The complete area definition was 14 km eastwest by 9 km northsouth with the southwest corner at 616000 m X, 3580000 m Y UTM NAD27, and is shown in red on Figure TFIELD 513. All other aspects of the automatic calibration, including the forward model and the SVD assist process, were left the same. The resulting fields had significantly better fits to the steadystate heads, and little impact was seen on the transient test results (Table TFIELD 55).
Figure TFIELD 513. Recalibration Boundary Shown in Red To the Northeast of the WIPP Site. Recalibration boundary limits are UTM X and Y NAD27 Zone 13 (m).
After examination of the acceptance criteria (discussed in the next section), some fields were recalibrated again, using the same recalibration process but holding none of the parameter values fixed at previously calculated values. This process essentially added some additional calibration iterations to these fields. This was only done on 15 fields that were now completely within the acceptance criteria for steadystate heads, and just outside the acceptance criteria for transient tests (Table TFIELD 55). The intent of this additional calibration was to increase the quality of the transient fits to get a total of 100 fields that met both the steadystate and transient calibration requirements. This secondary recalibration process was continued until 100 fields were obtained that met the requirements, and it did not always improve fits (i.e., in some cases the fields were already as fully calibrated, given the number of pilot points and observations and initial conditions).
Table TFIELD 55. Summary of Statistics Regarding Average Steadystate and Transient Errors Across Three Calibration Groups

Steadystate Average Error (cm) 
Pumping Response Average Error (cm) 



Minimum 
Average 
Maximum 
Minimum 
Average 
Maximum 
n 
A 
45.6 
68.9 
195.0 
12.9 
16.4 
34.2 
135 
B 
50.5 
73.8 
115.3 
12.5 
16.5 
23.0 
50 
C 
57.8 
66.1 
72.8 
13.9 
15.4 
17.4 
15 
A: Realizations with targeted recalibration near SNL8. B: Realizations with correct SNL8 data. C: Realizations recalibrated twice. 
Because some fields were calibrated only once (set B: 50 fields following correction of the steadystate values), some fields were calibrated once and then underwent a localized recalibration (set A: the 135 first fields calibrated), and some fields even underwent a second round of calibration (set C: 15 fields), the impact this may have had on the final selection of fields was evaluated. Summaries of statistics for these calibration groups are given in Table TFIELD 55, while Appendix E of Hart et al. (Hart et al. 2009) contains the complete list of fields.
Because the final selection process did not look at which set of fields the results were taken from, the mix of fields should be similar to a random selection if the calibration processes were producing equivalent results. The random selection of fields from set B can be modeled as a binomial distribution with the pvalue of 0.25 and n = 100. If the results are within the 95% confidence interval for a random selection of fields, then there should be between 17 and 33 fields selected from the set B. The final results used 83 fields from sets A and C, and the remaining 17 were selected from set B. This is within the confidence interval, so it is concluded that the different processes had no impact on the selection of the final fields. The selection of fields from set C versus those from set A can be modeled the same way, with a pvalue of 0.10 and n = 83. The final selection included 10 from set C, which is within the confidence interval of 3 to 13 fields, and again the calibration process did not impact the field selection.
Because this mix of final fields is acceptable and came strictly from the cutoff values, and not from any deliberate attempt to select from one group or another, all 100 fields meeting the acceptance criteria are equally good and equally probable representations of the Culebra.
The selection criteria for the "best" calibrated fields consisted of comparing the absolute error of the modeled steadystate heads (sum of absolute values of residuals between model and data) to a cutoff value, and comparing the absolute average error of the modeled transient responses (sum of absolute values of errors at individual observation wells averaged through time) to a cutoff value. The steadystate and transient criteria were evaluated separately, and only fields that were less than the cutoff value for both sets of tests were selected as the final fields. The final cutoff values used were the mean value of the errors taken across all 200 fields, and are presented in Table TFIELD 56. The cutoffs were selected to choose approximately half of the fields. Using the mean values resulted in a set of 102 fields, so the two fields with the largest sum of the two metrics were discarded. In Figure TFIELD 514, the sum of the steadystate average errors was graphed against the sum of the transient pumping tests' average errors, and the selected and unselected fields are shown. The trend line shows graphically how PEST allows tradeoffs while keeping the improvement in errors as balanced as possible. The final field IDs are presented in Table TFIELD 57.
Table TFIELD 56. Cutoff Values for Final Field Selection
Test Type 
Average Error Selection Cutoff 
Steady State 
0.699 m 
Transient Pumping Response 
0.164 m 
Table TFIELD 57. Final Selected Field Identifiers
r001 
r055 
r207 
r652 
r002 
r058 
r256 
r655 
r004 
r059 
r260 
r657 
r006 
r060 
r273 
r664 
r007 
r061 
r276 
r669 
r009 
r064 
r279 
r694 
r010 
r070 
r298 
r707 
r012 
r073 
r327 
r727 
r013 
r074 
r328 
r752 
r017 
r076 
r361 
r791 
r024 
r078 
r431 
r806 
r027 
r082 
r440 
r808 
r028 
r083 
r465 
r809 
r029 
r084 
r486 
r814 
r032 
r090 
r489 
r823 
r034 
r092 
r506 
r861 
r037 
r095 
r508 
r883 
r038 
r097 
r511 
r902 
r040 
r098 
r515 
r910 
r041 
r102 
r522 
r921 
r045 
r104 
r568 
r922 
r051 
r137 
r571 
r940 
r052 
r142 
r631 
r981 
r053 
r191 
r634 
r982 
r054 
r203 
r640 
r984 
Figure TFIELD 514. Selection of Best Fields From All Fields by Weighted Sum of Steadystate Errors and Sum of Average Pumping Test Average Errors
The postcalibration analysis consisted of performing statistical analyses on the selected fields, and examining the calibrated forward model outputs. The full results of the steadystate forward model outputs are presented in Appendix C of AP114 Task 7 (Hart et al. 2009), pumping test results are presented in Appendix D, and tabular results are presented in Appendix E of the same report. Calibrated model inputs and outputs for the 100 selected fields are presented in Attachment A to Appendix TFIELD.
Plots of mean and standard deviation of the final 100 fields are given in Section TFIELD5.4.1.1. The bulk T of the final calibrated fields are also compared to the bulk T of the base fields, using membership in the high or low T categories. Similarly, Sections TFIELD5.4.1.2 and TFIELD5.4.1.3 show summary statistics regarding the calibrated S and R fields, respectively.
The T values presented in this section are the effective T values (T _{e}), which include A. Effective transmissivity was calculated as
log_{10} T _{e} = log_{10} T _{EW} + ½ log_{10} A (TFIELD 5.3)
which is the average of log_{10}T in the northsouth and eastwest directions (see Equation TFIELD 5.1 in Section TFIELD5.2.1.2). The bulk T _{e}, which is the average log_{10} T _{e} value of all cells in a given zone or zones, was calculated for the central and Salado dissolution region (zones 02) and compared to the bulk T _{e} of the same zones from the base fields. The eastern, very lowT region (zone 3) was compared separately. The bulk T _{e} values are shown in Table TFIELD 58. The mean effective T _{e} and the standard deviation of T _{e} are presented in map form in Figure TFIELD 515 and Figure TFIELD 516. The mean effective transmissivity map does not show the very low T zone east of the halite margins to improve the colormap contrast across the area around the WIPP site.
Because pilot point parameter values were essentially unconstrained for T (they were allowed to change across 18 orders of magnitude as shown in Table TFIELD 53), some areas in zones 0 and 1 could change from a lowT zone into the range generally considered highT and vice versa. The defining value for highT was set in AP114 Task 5 (Hart et al. 2008) to be the bulk transmissivity value of the base fields: −5.41 log_{10} m^{2}/s. At each cell, the number of fields whose initial and final T values were in the highT zone was calculated, and the maps of those numbers for the base and calibrated fields are presented in Figure TFIELD 517 and Figure TFIELD 518, respectively. The total number of fields where transmissivity effectively changed zones is represented graphically in Figure TFIELD 519 and Figure TFIELD 520. In these figures, the white regions define areas where no fields changed groups. The two measures shown in these sets of maps provide an indication of how the geologically based conceptual model used to create the base fields was altered by the steadystate and transient hydraulic information.
Table TFIELD 58. Bulk Log_{10} T _{e} Values Comparison
Base field bulk log_{10} T _{e} (Zones 02) 
−5.41 log_{10} (m^{2}/s) 
Calibrated field bulk log_{10} T _{e} (Zones 02) 
−5.02 log_{10} (m^{2}/s) 
Base field bulk log_{10} T _{e} (Zone 3) 
−11.74 log_{10} (m^{2}/s) 
Calibrated field bulk log_{10} T _{e} (Zone 3) 
−10.47 log_{10} (m^{2}/s) 
Figure TFIELD 515. Mean Effective Transmissivity (T _{e}) for Zones 02 Across the 100 Final Selected Fields. All 100 calibrated T _{e} fields are plotted in Appendix TFIELD Attachment A.
Figure TFIELD 516. Standard Deviation of Effective Transmissivity (T _{e}) for All Zones Across the 100 Final Selected Fields. All 100 calibrated T _{e} fields are plotted in Appendix TFIELD Attachment A.
Figure TFIELD 517. HighT Zone Membership Calculated for the Base 100 T Fields Corresponding to the 100 Selected Calibrated Fields
Figure TFIELD 518. HighT Zone Membership Calculated for the Calibrated T Values
Figure TFIELD 519. Number of T Fields Where Low T Became High T Through PEST Calibration
Figure TFIELD 520. Number of T Fields Where High T Became Low T Through PEST Calibration
The mean and standard deviation of the final S fields are presented in Figure TFIELD 521 and Figure TFIELD 522. The mean S fields indicate that the overall S values in the confined and transitional zones did not change much from their initial values. Figure TFIELD 522 highlights the area northwest of P14 with a red dashed oval. This area has high variability in estimated S across the 100 selected realizations. This may have some relation to the relatively poorer fits of the transient data for the WIPP25 response to the P14 pumping test in the model.
Figure TFIELD 521. Mean Storativity Values Across the 100 Final Calibrated Fields. Individual S fields are plotted in Appendix TFIELD Attachment A.
Figure TFIELD 522. Standard Deviation of Storativity Values Across the 100 Final Calibrated Fields. Red oval shows P14 to WIPP25 area of influence. Individual S fields are plotted in Appendix TFIELD Attachment A.
The final recharge values were all less than the initial values of 10^{−11} m/s (3.2 × 10^{−4} meters per year (m/yr)). Compared to the other parameters, there was very little change in recharge. Because the recharge zone was linear, in addition to the cellbycell mapping, a view of the average, minimum and maximum recharge values is shown as a cross section in the crossdirection (across a row) as if looking from the south to the north through the domain in Figure TFIELD 523.
Figure TFIELD 523. Recharge as Viewed Through Columns From the South. The initial value was set at 10^{−3.5} m/year. The sharp dropoff to the west is the transition to the single fixedrecharge point of 10^{−11.5} m/year (interpreted as zero by REAL2MOD).
The two main divisions of the results are the steadystate head results and the pumping test results. The results presented here only represent the 100 final selected fields, and therefore the maximum error is limited by the selection criteria described in Section TFIELD5.3.4: an average steadystate error of less than 0.699 m and an average pumping test observation error of less than 0.164 m.
Figure TFIELD 524 shows the modeled steadystate head values plotted against the measured head values. The onetoone correspondence line shows the ideal match, and the modeled results are presented as boxandwhisker plots at each observation well. Figure TFIELD 525 shows all 4,200 head errors across all 100 fields as a histogram of error values for steadystate head. Additional figures and tables summarizing steadystate calibration are presented in Appendix C of Hart et al. (Hart et al. 2009). The modelmeasurement misfit can be modeled as a zeromean Gaussian distribution with a standard deviation of 0.10 m (McKenna and Wahi 2006). The measurement error distribution curve is included in Figure TFIELD 525.
Graphs for each of the transient pumping test results are presented in Appendix D of Hart et al. (2009). The average value of the error in the final fields ranged from 0.12 m to 0.164 m across all tests, with an average error of 0.15 m. The maximum error for a single observation well ranged from 0.005 m to 2.5 m, with an average of 0.36 m as the maximum error at a given observation well.
Figure TFIELD 524. Results for 42 Total Steadystate Head Measurements for the 100 Selected Fields (No SNL6 or SNL15). Observed heads are red ×'s along the diagonal line.
Figure TFIELD 525. Histogram of Steadystate Head Errors for the 100 Selected Fields (No SNL6 or SNL15). Red dashed line is the ±3σ section of the measurement error PDF. The slight skew to right is an artifact of the binning.
The work described in Section TFIELD6.0 was completed under AP144, Analysis Plan for the Calculation of Culebra Flow and Transport for CRA2009 PABC (Kuhlman 2009). The modifications used for CRA2009 PABC still apply to CRA 2014, and are therefore discussed here.
PA models two categories of miningimpacted transmissivity fields: partial mining with only mining outside the LWB, and full mining with regions both inside and outside the LWB mined.
The CRA2009 PABC Culebra Tfield mining modifications basically follow the procedure used in CRA2004 PABC (Lowry and Kanney 2005), with two exceptions: 1) a new definition of the region containing minable potash is used, and 2) the new Tfields in Sections TFIELD3.0 through TFIELD5.0 are used as inputs. The procedure for the mining modification portion of the analysis is summarized below:
1. Obtain the sampled values for the random mining modification factor (100 vectors × 3 replicates);
2. Map potential areas of future potash mining onto the groundwater modeling domain for both full and partialmining scenarios;
3. Apply the mining modification factor to the 100 stochastically calibrated Tfields from AP114 Task 7 (Hart et al. 2009), producing 600 miningmodified Tfields (100 vectors × 2 mining scenarios × 3 replicates);
4. Perform steadystate flow simulations for each of the 600 miningmodified Tfields using MODFLOW2000; and
5. Perform particle tracking using the new miningaffected flowfields to determine advective travel times to the WIPP LWB.
Potash mining in the region surrounding the WIPP involves underground excavation in the McNutt Potash zone of the upper Salado Formation, which is located stratigraphically above the WIPP repository horizon but below the Culebra Member of the Rustler Formation (see Figure TFIELD 21). It is hypothesized that subsidence due to collapse of the underground voids created in the McNutt potash zone during mining will lead to increased permeability in the Rustler Formation, due to increased fracturing, similar to the Salado dissolution zone effects in Figure TFIELD 22. The purpose of the mining scenario calculations is to determine the impact of potash mining on groundwater flow directions and transport velocities in the Culebra. This analysis largely represents a reapplication of the methods used in CRA2004 PABC (Lowry and Kanney 2005), with a few minor exceptions:
1. The definition of the regions where minable potash is believed to exist has been updated, as obtained from the Bureau of Land Management (BLM) (Cranston 2009).
2. The configuration of the MODFLOW model that mining modifications are being applied to has changed (see Sections TFIELD3.0 through TFIELD5.0).
3. The way the miningmodified areas interact with specified head areas of the flow model has changed, due to the change in the boundary conditions (there were no specified head areas inside the CRA2004 PABC MODFLOW model).
This section describes the CRA2009 PABC effort in characterizing mining effects in the Culebra and highlights the differences and additions relative to past calculations (Ramsey and Wallace 1996; Lowry 2003a, Lowry 2003b, and Lowry 2004). The reader is encouraged to review the past documents for further background information. There have been no changes to the mining modifications procedure or data for CRA2014 since CRA2009 PABC.
Starting with the 100 calibrated Tfields from Section TFIELD5.0, Tfields are modified to reflect the effects of mining by multiplying the transmissivity value in cells that lie within designated mining zones by a random factor uniformly sampled between 1 and 1000. The range of this factor is set by the U.S. Environmental Protection Agency (EPA) in regulation 40 CFR 194.32(b) (U.S. EPA 1996). The scaling factor for each Tfield is provided by Latin hypercube sampling (Kirchner 2010).
A forward steadystate flow simulation is run for each new Tfield under each mining scenario (full and partial) across three replicates of mining factors, resulting in 600 simulations. Particle tracking is performed on both the 100 original and 600 modified flowfields to compare the flow path and groundwater travel time from a point above the center of the WIPP disposal panels to the LWB. Cumulative distribution functions (CDFs) are produced for each mining scenario and compared to the undisturbed scenario. The CDFs indicate the probability a conservative tracer (i.e., a marked water particle) would reach the WIPP LWB during a given interval of time, flowing in the Culebra under a natural gradient. In addition to comparing travel times, particletracking directions are also examined to determine the effect on the regional flow direction in the WIPP area due to mining.
The eastern limit of the MODFLOW model domain used in the CRA2009 PABC analysis (Hart et al. 2008) was extended eastward, compared to the MODFLOW domain used in the CRA2004 PABC analysis. This change was made to locate the boundary in an area where halite is present in all of the nondolomite members of the Rustler Formation, simplifying the specification of the eastern model boundary condition. See Section TFIELD3.1 for CRA2009 PABC MODFLOW model construction details.
Like the model domain and discretization, the boundary and initial conditions used in CRA2009 PABC are described fully in AP114 Task 7 (Hart et al. 2009). Regional flow rates within the flow model are controlled by the boundary conditions and the hydraulic conductivity distribution. The regional gradient across the domain was approximately
(943.9 m  911.6 m)/30.7 km = 0.00105 (TFIELD 6.1)
The gradient across the model domain was computed by averaging the constant heads along the northernmost portion of the northern boundary (row 1, columns 1140, 943.9 m), subtracting the average heads along the entire southern boundary (911.6 m), and dividing by the northsouth model domain extent (30.7 km). It was assumed that mining impacts would not significantly change this regional gradient and thus the specified initial conditions for the mining scenarios are identical to those in AP114 Task 7 (Hart et al. 2009). In addition, the CCA, CRA2004, and CRA2004 PABC all used this same conceptualization (keeping the outer boundary conditions fixed between the mining and nonmining scenarios). The same conceptualization was maintained in the CRA2009 PABC model to facilitate comparisons between the different models.
Figure TFIELD 61. Comparison of Minable Potash Area to the Flow and Transport Modeling Domains. Green hatches represent minable potash area (Cranston 2009).
The 2009 version of the BLM map delineating distribution of minable potash ore was obtained from BLM as an ESRI shapefile (Cranston 2009), plotted in Figure TFIELD 61. The process to convert this shapefile to a grid of integers corresponding to Culebra MODFLOW model cells (indicating whether a model cell was affected by mining or not) is explained in Appendix 1 of Kuhlman (Kuhlman 2010).
Since the potashmining horizon is located in the Salado Formation below the Culebra, the areas disturbed by mining activities in the Culebra are larger than mined areas in the Salado due to angleofdraw effects. Subsidence effects will not propagate up vertically, but instead will propagate up and out at 45° angles between horizontal and vertical. Based on an average distance from the McNutt potash zone to the Culebra, a 253mwide collar was added to the miningimpacted areas (consistent with that done previously; see Ramsey and Wallace (Ramsey and Wallace 1996) and Bertram (Bertram 1995). This was considered a conservative estimate of the angleofdraw effects. To accommodate the angle of draw, the mining zone boundaries, as overlaid on the current model grid, were extended outward three cells (300 m) in the x and ydirections, and two cells (283 m) in the diagonal directions (see Figure TFIELD 62 for an illustration of the miningexpansion stencil).



A 




A 
A 
A 
A 
A 


A 
A 
A 
A 
A 

A 
A 
A 
M 
A 
A 
A 

A 
A 
A 
A 
A 


A 
A 
A 
A 
A 




A 



Figure TFIELD 62. Stencil Used to Model Cells Affected by Miningrelated Subsidence (Blue Cells with A) Due to Mining in Red "M" Cell, Using 45° Angle of Draw
Figure TFIELD 63 shows the CRA2009 PABC modeling domain and mining zones for the fullmining case in comparison to the 1996 CCA and the CRA2004 delineations. The comparison of the current and previous partialmining cases is shown in Figure TFIELD 64. A closeup of the WIPP site and the distribution of minable potash is shown in Figure TFIELD 65, which illustrates how the definition inside the WIPP LWB has changed significantly since CRA2004 PABC. For CRA2004 PABC, the closest minable potash was approximately 1,230 m from the center of the WIPP panels in the southeast direction; for CRA2009 PABC, this distance was reduced to approximately 670 m (in a more easterly direction).
The output of this miningarea delineation was a mining effects indicator field. A value of 1 means the model cell lies within a potential miningaffected zone, and a value of 0 means that it is outside a potential miningaffected zone.
Figure TFIELD 63. Definitions of Miningaffected Areas in Fullmining Scenario Between Current and Previous Models. Base image is Figure 3.2 from Lowry and Kanney (2005). CRA2009 PABC mining area (red stippled area) and model domain (green line) definitions are current definitions used in CRA2014.
Figure TFIELD 64. Definitions of Partialminingaffected Areas Between Current and Previous Applications. Base image is Figure 3.3 from Lowry and Kanney (2005). CRA2009 PABC mining area (red stippled area) and model domain (green line) definitions are current definitions used in CRA2014.
Figure TFIELD 65. Comparison of Minable Potash Distribution Inside the WIPP LWB for CRA2004 PABC (Dark Gray) and CRA2009 PABC (Translucent Green). The WIPP repository plan is shown for comparison, from Figure 3.6 of Lowry and Kanney (2005).
The calibration process in Section TFIELD5.0 produced 100 sets of T, A, S, and R fields that each minimize the error between observed and modelcalculated head distributions. To simulate the effects of mining, the field of T values from each realization was multiplied by its own unique mining scaling factor in areas of potential mining, and MODFLOW was rerun with these miningmodified Tfields to produce the miningaffected head and flow distributions. The other parameter fields (S, A, and R) were not modified in the process. The cellbycell flow budget files were used in particle tracking and radionuclide transport calculations. Three different sets of mining factors were used, each set forming a replicate (given here as R1, R2, and R3). Thus, for each mining scenario (full and partial), three sets of 100 miningaltered Tfields and related cellbycell flow budgets were produced.
In each realization, a single conservative particle was tracked from the UTM NAD27 coordinate
x = 6,135,975 m, y = 35,813,852 m (i.e., the Culebra location of monitoring well C2737, directly above the center of the WIPP waste panels) to the LWB for each combination of Tfield, replicate, and mining scenario using the SNLdeveloped particletracking code DTRKMF. Two main outputs were generated from the suite of particle tracks. First, plots were constructed showing the individual tracks for all 100 Tfields in each scenario for each replicate (six plots total). This allows visual comparison of the prevailing flow directions for the full and partialmining scenarios and the qualitative comparison of the variability of the tracking direction. Secondly, CDFs were constructed for each replicate and scenario, which describe the probability that a water particle will cross the LWB in a given amount of time. The six plots and the CDFs are presented in Section TFIELD6.3.
Particle tracks were computed using DTRKMF (Rudeen 2003), which uses the binary cellbycell flow budget files produced by MODFLOW2000. In flow calculations, the full 7.75 m thickness of the Culebra is used, while for transport and particletracking purposes the thickness is reduced to 4.0 m to focus all flow through the lower, more permeable portion of the Culebra (Holt 1997). A value of 16% porosity was used for the particletracking calculations (parameter DPOROS). Porosity directly affects transport, but is not needed for the calibration of the flow model.
Particle tracking was performed to allow plotting and quantitative comparison between the two mining scenarios and the nonmining scenario, which was not used in the PA SECOTP2D radionuclide transport calculations. The particletracking results illustrate the advective pathway taken by a marked water particle. They do not take into consideration retardation, dispersion, or molecular diffusion, which may be accounted for in PA by radionuclide transport. The particle tracks also allow easier comparison of the 600 results (each a 1D trace) in a single plot, in contrast to showing 600 plots of concentrations (each a 2D field) produced from SECOTP2D.
Compared to the nonmining scenario (results already given in AP114 Task 7 (Hart et al. 2009)), the travel times for the partialmining scenarios are longer, while travel times for the fullmining scenarios are shorter. The median travel time across all three replicates for the fullmining scenario is approximately 0.689 times the median travel time of the nonmined scenario (see Table TFIELD 61, Figure TFIELD 66, and Figure TFIELD 67 for summary statistics and comparison to CRA2004 PABC results). All advective particle travel times are plotted, but it should be noted that the regulatory limit for radionuclide transport modeling is 10,000 years, taking into consideration retardation, diffusion, and dispersion (which do not apply to particle track modeling). The median travel time across all three replicates for the partialmining scenario is 3.034 times greater than for the nonmining scenario. For CRA2004 PABC, travel times in both the full and partialmining scenarios were slower (longer) than for the nonmining scenario. The CDFs for the full, partial, and nonmining scenarios are shown in Figure TFIELD 66. Travel times from CRA2009 PABC for particles in the nonmining scenario are closer to CCA travel times than those from CRA2004 PABC (Figure TFIELD 67).
Table TFIELD 61. Particletracking Traveltime Statistics from Center of the WIPP Panels to the WIPP LWB (Years). Global statistics for full and partial mining include 300 realizations, while no mining only includes 100 realizations.

CRA2009 PABC 
CRA2004 PABC 

Replicate 
Statistic 
Full 
Partial 
No Mining 
Full 
Partial 
No Mining 
R1 
Median 
5,138 
22,581 
N/A 
64,026 
117,815 
N/A 
Max 
200,260 
91,119 
2,175,165 
2,727,191 

Min 
1,591 
5,042 
2,130 
5,185 

R2 
Median 
4,956 
21,999 
80,801 
148,489 

Max 
94,852 
84,929 
2,059,263 
1,667,084 

Min 
1,421 
5,037 
2,463 
4,855 

R3 
Median 
5,560 
22,537 
74,315 
118,919 

Max 
93,172 
86,758 
1,779,512 
3,128,693 

Min 
1,421 
4,505 
2,507 
3,314 

Global 
Median 
5,084 
22,376 
7,374 
70,170 
131,705 
18,289 
Max 
200,260 
91,119 
73,912 
2,175,165 
3,128,693 
101,205 

Min 
1,421 
4,505 
2,618 
2,130 
3,314 
3,111 
Figure TFIELD 66. CDF of Advective Particle Travel Times From the Center of the WIPP Waste Panels To the WIPP LWB for Full, Partial, and Nonmining Scenarios
Figure TFIELD 67. Comparison of Advective Particle Travel Time CDFs for Nonmining Scenarios of CRA2009 PABC, CRA2004 PABC, and CCA. Travel times are from the center of the WIPP waste panels to the WIPP LWB.
The particle track directions for the non, full, and partialmining scenarios for CRA2009 PABC are illustrated in Figure TFIELD 68 to Figure TFIELD 611. Figure TFIELD 613 shows the nonmining case particle tracks all the way to the edge of the MODFLOW model domain, rather than only to the WIPP LWB. Like past mining scenario calculations (i.e., CRA2004 PABC), there is a strong similarity between the three replicates (R1, R2, and R3) for each scenario (full or partial mining), although the travel directions for the CRA2009 PABC are different than for the CRA2004 PABC (Lowry and Kanney 2005). A larger amount of minable ore exists inside the WIPP LWB, especially the ore immediately to the east of the particle release point. This leads to different effects of full mining on travel times compared to CRA2004 PABC.
The highT pathway in the southeastern portion of the WIPP site (Figure TFIELD 515) was more accurately represented in CRA2009 PABC results, compared to CRA2004 PABC results. This difference in the underlying flow field calibration is another source of difference between CRA2004 PABC and CRA2009 PABC particletracking results.
In CRA2009 PABC, nearly all particles immediately go east to this boundary and then move south along it towards to the edge of the LWB at approximately x = 612.75 km (see Figure TFIELD 69 and Figure TFIELD 612). This is in contrast to the partialmining scenario where the tracking directions are more similar to the nonmining scenario, but more evenly distributed spatially along the southern boundary. In the nonmining scenario, most of the particles exit near the hightransmissivity zone at approximately x = 615 km.
Figure TFIELD 68. 100 Particle Tracks for Nonmining Scenario. Black box is the WIPP LWB, green circles are Culebra monitoring well locations.
Figure TFIELD 69. 100 Particle Tracks Each for R1 Full and Partial Mining Scenarios. Small magenta squares and blue crosses indicate centers of MODFLOW cells located within potash and miningaffected areas, respectively; thin black line is minable potash.
Figure TFIELD 610. 100 Particle Tracks Each for R2 Full and Partialmining Scenarios. Small magenta squares and blue crosses indicate centers of MODFLOW cells located within potash and miningaffected areas, respectively; thin black line is minable potash.
Figure TFIELD 611. 100 Particle Tracks Each for R3 Full and Partialmining Scenarios. Small magenta squares and blue crosses indicate centers of MODFLOW cells located within potash and miningaffected areas, respectively; thin black line is minable potash.
HighT areas corresponding to the miningaffected zones create preferential pathways through the system (e.g., see oranges and yellows in Figure TFIELD 614). These preferential pathways result in higher velocities and flow rates through the mining zone and therefore relatively slower velocities in the nonmined areas. In the partialmining scenario, where there is no mining inside the WIPP LWB, the preferential pathway goes "around" the LWB, rather than through it (similar to behavior seen in both mining scenarios for CRA2004 PABC). In the fullmining scenario, the potentially mined regions are closer to the release point than in CRA2004 PABC (see Figure TFIELD 65 for comparison), giving the particles a highT pathway from the release point to the LWB, resulting in shorter travel times than the nonmined scenario. This behavior is different from that predicted with the CRA2004 PABC model. A comparison of the median, maximum, and minimum travel times for each scenario is presented in Table TFIELD 61.
Figure TFIELD 612. Histograms of Particle xcoordinates at Exit Point From LWB. Full and partialmining include all three replicates (note different vertical scales between plots; no mining contains 100 particles while mining scenarios each include 300 particles).
Figure TFIELD 613. Particle Counts in Each Cell Across All 100 Selected Realizations for Nonmining Scenario
Figure TFIELD 614. Magnitude of Darcy Flux for a Single Realization (r440) for No, Partial, and Fullmining Scenarios Using Cellbased Coordinates. LWB (black) and SECOTP2D transport model domains (red) shown for reference.
Instantaneous speeds (the magnitude of particle velocities) were calculated from the DTRKMF particle locations and times using backwards finite differences,
where a subscript i indicates the previous time step (a record or line in the DTRKMF output file) and a subscript i+1 is the current time step. This approach assumes a straight line connects the locations at the beginning and ends of the step, so it is potentially underestimating speeds, but step sizes are small and error should be minimal. These values should be used for qualitative comparisons between realizations and scenarios, rather than quantitative estimates of true particle velocities.
In Figure TFIELD 615 through Figure TFIELD 618, the color of the diamond indicates the particle velocity; the dots are located at the midpoint of the step, e.g.,
x
_{midpt} = ½[x(t_{i}
) + x(t_{i}
_{+1})], y
_{midpt} = ½[y(t_{i}
) + y(t_{i}
_{+1})]. In the nomining case (Figure TFIELD 615), the highest particle velocities are in the southeastern portion of the particle swarm, corresponding to the hightransmissivity pathway (Hart et al. 2009) exiting the LWB at approximately
x = 614,750 m (Figure TFIELD 612). This highT pathway has been observed in multiwell pumping tests, and was intentionally included in the soft data used to create the base Tfields (Section TFIELD4.2.3). This highT pathway was not as prevalent in the CRA2004 PABC model calibration results.
The effects of the fullmining scenario are clearly evident in the left halves of Figure TFIELD 616 through Figure TFIELD 618. High particle velocities (yellows and oranges) are found along the margin of the miningaffected areas, where particles enter the increasedT region. For comparison, in the partialmining scenario the particles are relatively slowed down and more evenly distributed in the region between the release point and the southern WIPP LWB, with the only high velocities found in the highT pathway in the southeast, similar to the nomining scenario (Figure TFIELD 615).
Figure TFIELD 615. Particle Speeds for Nonmining Scenario Computed from DTRKMF Results. Open symbols are Culebra well locations.
Figure TFIELD 616. Particle Speeds for R1, Computed From DTRKMF Results. Open symbols are Culebra well locations.
Figure TFIELD 617. Particle Speeds for R2, Computed From DTRKMF Results. Open symbols are Culebra well locations.
Figure TFIELD 618. Particle Speeds for R3, Computed from DTRKMF Results. Open symbols are Culebra well locations.
Correlation analysis for the CRA2009 PABC particletracking calculations shows that travel time and the random mining factor have weak positive correlation with the fullmining (Figure TFIELD 619) or partialmining (Figure TFIELD 620) scenarios. Larger mining factors are weakly correlated with longer travel times. This is similar to what was observed for CRA2004 PABC (Lowry and Kanney 2005). The high scatter in Figure TFIELD 619 and Figure TFIELD 620 indicates that the transmissivity spatial distribution plays the more significant role in determining the travel time than the mining factor does. See Appendix 1 of Kuhlman (Kuhlman 2010) for a crosssectional comparison of transmissivity for each mining type, showing that the variability in the transmissivity due to calibration is on the same order as that due to mining for a single realization. The mining factor plays a weak but slightly larger role in explaining the observed variance for the partialmining realizations (Figure TFIELD 620) than the fullmining realizations (Figure TFIELD 619), based on the larger (but still relatively small) R^{2} value for the straightline fit of log_{10} travel times to mining factors.
Figure TFIELD 619. Correlation of Mining Factor and Travel Time to the WIPP LWB for Fullmining Scenario (All Replicates)
Figure TFIELD 620. Correlation of Mining Factor and Travel Time to the WIPP LWB for Partialmining Scenario (All Replicates)
The 100 transmissivity fields resulting from calibration to both steadystate and transient observed freshwater heads in the Culebra (Section TFIELD5.1) were modified to account for potential effects due to mining potash from the Salado Formation above the repository. A definition of the areal extent of minable potash was obtained from the BLM (Cranston 2009) and used to define areas where Culebra transmissivity was increased by a randomly sampled mining factor (1 ≤ MINP_FACT ≤ 1000). Two mining scenarios were developed: a fullmining scenario with all minable potash removed and a partialmining scenario with only potash outside the WIPP LWB removed.
The miningmodified transmissivities were inputs to a MODFLOW flow model, which produced budget files used by DTRKMF to compute advective particle tracks from a release point at the center of the WIPP waste panels (C2737) to the edge of the WIPP LWB. Results show that for the partialmining scenario, the median particle travel time of 22,376 years is 3.03 times greater than for the nonmining scenario (7,374 years); the median particle travel time for the partialmining scenario in CRA2004 PABC was 7.06 times greater than for the nonmining scenario. In contrast to the CRA2004 PABC, the fullmining scenario decreased median travel time to 5,084 years, a factor of 1.45 faster than for the nonmining scenario; the median particle travel time for the fullmining scenario in CRA2004 PABC was 3.84 times slower than for the nonmining scenario. For the partialmining scenario, the increase in transmissivity due to mining increases the relative flow rate through the mining zones, with a corresponding decrease in flow through the nonmining zones. This decrease in flow through the nonmining zones produces longer travel times for the partialmining scenario. For the fullmining scenario, the potash definition from BLM (Cranston 2009) locates minable potash ore much closer to the C2737 release point than in CRA2004 PABC (see Figure TFIELD 65). This new shortened distance from the release point to the nearest minable potash (in the fullmining scenario) reverses the slowingdown trend observed in the CRA2004 PABC analysis.
As in the CRA2004 PABC calculations, a very weak positive correlation was found between the travel times and the random mining factor (the higher the random mining factor, the longer the particle travel time  see Figure TFIELD 619 and Figure TFIELD 620). As the mining factor is increased, the flow through the nonmining areas (including the C2737 particle release location) is decreased, producing longer travel times and the positive correlation. Most of the advective particle travel time variability is due to differences in the base Tfields and their subsequent calibration, and not the random mining factor.
Observed Culebra transmissivities (T) have been related to four deterministic factors: the thickness of overburden above the Culebra, the presence or absence of dissolution of the upper Salado, the presence of gypsum cements, and the presence or absence of halite in units above and below the Culebra. Culebra T is also related to the occurrence of open, interconnected fractures that cannot be mapped as easily as the other three factors and therefore must be treated stochastically. A linearregression model for Culebra T has been developed based on these factors that provided an excellent match to the observed data, and can be tested through the collection of additional data. This model was used to create 1000 stochastic realizations of the distribution of Culebra T (base Tfields) in the vicinity of the WIPP site.
A MODFLOW2000 modeling domain was defined extending 30.7 km northsouth and 28.4 km eastwest, roughly centered on the WIPP site. This domain was discretized into 87,188 uniform 100m square twodimensional finitedifference cells. An inactive portion of the northwest corner of the domain is used to represent a noflow boundary along the axis of Nash Draw. A lowpermeability constanthead portion of the eastern section of the domain is used to represent the lithostatic pressure portion of the Culebra sandwiched above and below by Rustler halite units. Freshwater head observations in 42 monitoring wells from May 2007 were used as steadystate calibration targets. Drawdown observations in 62 observation wells, in response to 9 unique pumping tests, were used as transient calibration targets. A subset consisting of 100 of the 200 calibrated Culebra model realizations were selected based on their ability to simulate these observed heads.
The EPA requires that the potential effects of future potash mining be taken into account when evaluating the performance of the WIPP disposal system. Accordingly, T in the areas within the model domain where current or future mining might affect the Culebra were scaled by a random multiplier between 1 and 1,000 obtained from Latin hypercube sampling. A single multiplier was used for each Tfield, applied first to the areas outside the WIPP LWB that might be mined to create a partialmining Tfield, and then to all areas (both inside and outside the WIPP LWB) to create a fullmining Tfield. Three statistically similar replicates of mining multipliers were generated, leading to a total of 600 unique Tfields (100 calibrated realizations, 2 mining scenarios, and 3 replicates). The MODFLOW2000 flow budgets were used from each realization as input for both advective particle tracking (DTRKMF) summarized here and radionuclide solute transport (SECOTP2D) used in WIPP PA.
The nonmined travel times from the center of the WIPP waste panels to the WIPP LWB are similar to those computed for the CCA and therefore faster than those computed for CRA2004 PABC. The decrease in travel time to the LWB can be attributed to the presence of a consistent hightransmissivity pathway leaving the southeast portion of the LWB. The presence of this pathway is supported by observed drawdown data from the SNL14 pumping test.
In the partialmining case, particle tracks show increased travel times from the center of the WIPP waste panels to the WIPP LWB, compared to the nonmining scenario. In the fullmining case, particle tracks showed decreased travel times to the WIPP LWB, due to the close proximity of minable potash to the center of the WIPP waste panels.
(*Indicates a reference that has not been previously submitted.)
Beauheim, R.L. 1987a. Analysis of Pumping Tests of the Culebra Dolomite Conducted at the H3 Hydropad at the Waste Isolation Pilot Plant (WIPP) Site. Albuquerque, NM: Sandia National Laboratories. SAND862311. [PDF / Author]
Beauheim, R.L. 1987b. Interpretation of the WIPP13 Multipad Pumping Test of the Culebra Dolomite at the Waste Isolation Pilot Plant (WIPP) Site. Albuquerque, NM: Sandia National Laboratories. SAND872456. [PDF / Author]
Beauheim, R.L. 1989. Interpretation of H11b4 Hydraulic Tests and the H11 Multipad Pumping Test of the Culebra Dolomite at the Waste Isolation Pilot Plant (WIPP) Site. Albuquerque, NM: Sandia National Laboratories. SAND890536. [PDF / Author]
Beauheim, R.L. 2007. Diffusivity mapping of fracture interconnections. In Proceedings of the 2007 US EPA/NGWA Fractured Rock Conference. Westerville, OH, 2007. National Ground Water Association. [Author]
Beauheim, R.L. 2008. Analysis Plan for Evaluation and Recalibration of Culebra Transmissivity Fields AP114, Revision 1. Carlsbad, NM: Sandia National Laboratories. [PDF / Author]
Beauheim, R.L. 2009. Changes to Culebra TField Calibration Procedures under AP114 Task 7. Carlsbad, NM: Sandia National Laboratories. ERMS 551437. [PDF / Author]
Beauheim, R.L. and R.M. Holt. 1990. Geological and Hydrological Studies of Evaporites in the Northern Delaware Bains for the Waste Isolation Pilot Plant (WIPP). In D.W. Powers, R.M. Holt, R.L. Beauheim and N. Rempe, eds. Guidebook 14. Geological Society of America (Dallas Geological Society). pp.4578. [Author]
Beauheim, R.L. and G.J. Ruskauff. 1998. Analysis of Hydraulic Tests of the Culebra and Magenta Dolomites and Dewey Lake Redbeds Conducted at the Waste Isolation Pilot Plant Site. Albuquerque, NM: Sandia National Laboratories. SAND980049. [PDF / Author]
Bertram, S.G. 1995. Record of FEP Screening Work, FEP ID#NS11: Subsidence Associated with Mining Inside or Outside the Controlled Area. Carlsbad, NM: Sandia National Laboratories. ERMS 230761. [PDF / Author]
Bowman, D.O. and R.M. Roberts. 2008. Analysis Report for AP070, Analysis of Hydraulic Tests Performed in IMC461, SNL6, H11b2, H15, and C2737. Carlsbad, NM: Sandia National Laboratories. ERMS 539221. [Author]
Burgess, A., T. Doe, T. Lowenstein, and J.A. Thies. 2008. Waste Isolation Pilot Plant Culebra Hydrogeologicy Conceptual Model Peer Review Final Report. Carlsbad, NM: U.S. Department of Energy. ERMS 551513. [PDF / Author]
Chavez, M. 2006. NP 191 Software Requirements, Revision 12. Carlsbad, NM: Sandia National Laboratories. [PDF / Author]
Clayton, D.J., R.C. Camphouse, J.W. Garner, A.E. Ismail, T.B. Kirchner, K.L. Kuhlman, and M.B. Nemer. 2010. Summary Report of the CRA2009 Performance Assessment Baseline Calculation. Carlsbad, NM: Sandia National Laboratories. ERMS 553039. [PDF / Author]
Cranston, C.C. 2009. Minable Potash Ore. Carlsbad, NM: Bureau of Land Management. ERMS 551120. [PDF / Author]
Currie, J.B. and S.O. Nwachukwu. 1974. Evidence on incipient fracture porosity in reservoir rocks at depth. Bulletin of Canadian Petroleum Geology. 22, pp.4258. [Author]
Deutsch, C.V. and A.G. Journel. 1998. GSLIB: Geostatistical Software Library and User's Guide. 2nd ed. New York: Oxford University Press. [Author]
Doherty, J. 2000. PEST Manual. Brisbane, Australia: Watermark Numerical Computing. [Author]
Harbaugh, A.W., E.R. Banta, M.C. Hill, and M.G. McDonald 2000. MODFLOW 2000: The U.S. Geological Survey Modular GroundWater Model User Guide. Reston, VA: U.S. Geological Survey. OFR 0092. [PDF / Author]
Hart, D.B., R.L. Beauheim, and S.A. McKenna. 2009. Analysis Report for Task 7 of AP114: Calibration of Culebra Transmissivity Fields. Carlsbad, NM: Sandia National Laboratories. ERMS 552391. [PDF / Author]
Hart, D.B., R.M. Holt, and S.A. McKenna. 2008. Analysis Report for Task 5 of AP114: Generation of Revised Base Transmissivity Fields. Carlsbad, NM: Sandia National Laboratories. ERMS 549597. [PDF / Author]
Holt, R.M. 1997. Conceptual Model for Transport Processes in the Culebra Dolomite Member of the Rustler Formation. Albuquerque, NM: Sandia National Laboratories. SAND970194. [PDF / Author]
Holt, R.M., R.L. Beauheim, and D.W. Powers. 2005. Predicting fractured zones in the Culebra Dolomite. In B. Faybishenko, P.A. Witherspoon and J. Gale, eds. Dynamics of Fluids and Transport in Fractured Rock. American Geophysical Union. pp.10316. [Author]
Holt, R.M. and D.W. Powers. 1988. Facies variability and postdepositional alteration within the Rustler Formation in the vicinity of the Waste Isolation Pilot Plant, southeastern New Mexico. Carlsbad, NM: Department of Energy. WIPPDOE88004. [PDF / Author]
Holt, R.M. and L. Yarbrough. 2002. Analysis Report, Task 2 of AP088, Estimating Base Transmissivity Fields. Carlsbad, NM: Sandia National Laboratories. ERMS 523889. [PDF / Author]
Johnson, P.B. 2009a. Potentiometric Surface, Adjusted to Equivalent Freshwater Heads, of the Culebra Dolomite Member of the Rustler Formation near the WIPP Site, Revision 1. Carlsbad, NM: Sandia National Laboratories. ERMS 548746. [PDF / Author]
Johnson, P.B. 2009b. Routine Calculations Report In Support of Task 6 of AP114: Potentiometric Surface, Adjusted to Equivalent Freshwater Heads, of the Culebra Dolomite member of the Rustler Formation near the WIPP Site, Revision 2. Carlsbad, NM: Sandia National Laboratories. ERMS 551116. [PDF / Author]
Kirchner, T.B. 2010. Generation of the LHS Samples for the AP145 (PABC09) PA Calculations. Carlsbad, NM: Sandia National Laboratories. ERMS 552905. [PDF / Author]
Kuhlman, K.L. 2009. AP144 Analysis Plan for the Calculation of Culebra Flow and Transport for CRA2009PABC. Carlsbad, NM: Sandia National Laboratories. ERMS 551676. [PDF / Author]
Kuhlman, K.L. 2010. Analysis Report for the CRA2009 PABC Culebra Flow and Transport Calculations (AP144). Carlsbad, NM: Sandia National Laboratories. ERMS 552951. [PDF / Author]
Leigh, C.D., J.F. Kanney, L.H. Brush, J.W. Garner, G.R. Kirkes, T. Lowry, M.B. Nemer, J.S. Stein, E.D. Vugrin, S. Wagner, and T.B. Kirchner. 2004. 2004 Compliance Recertification Application Performance Assessment Baseline Calculation, Revision 0. Carlsbad, NM: Sandia National Laboratories. ERMS 541521. [PDF / Author]
Lowry, T.L. 2003a. Analysis Report, Task 5 of AP088: Evaluation of Mining Scenarios. Carlsbad, NM: Sandia National Laboratories. ERMS 531138. [PDF / Author]
Lowry, T.L. 2003b. Analysis Report, Tasks 2 & 3 of AP100: Grid Size Conversion and Generation of SECOTP2D Input. Carlsbad, NM: Sandia National Laboratories. ERMS 531137. [PDF / Author]
Lowry, T.S. 2004. Analysis Report for Inclusion of Omitted Areas in Mining Transmissivity Calculations in Response to EPA Comment G11. Carlsbad, NM: Sandia National Laboratories. ERMS 538218. [PDF / Author]
Lowry, T.S. and J.F. Kanney. 2005. Analysis Report for the CRA2004 PABC Culebra Flow and Transport Calculations. Carlsbad, NM: Sandia National Laboratories. ERMS 541508. [PDF / Author]
McKenna, S.A. and D.B. Hart. 2003. Analysis report for Task 4 of AP088: Conditioning of Base T Fields to Transient Heads. Carlsbad, NM: Sandia National Laboratories. ERMS 531124. [PDF / Author]
McKenna, S.A. and A. Wahi. 2006. Local Hydraulic Gradient Estimator Analysis of LongTerm Monitoring Networks. Ground Water. 44(5), pp.72331. [Author]
Mehl, S.W. 2001. MODFLOW2000, The US Geological Survey Modular GroundWater Model User Guide to the LinkAMG (LMG) Package for Solving Matrix Equations using an Algebraic Multigrid Solver. USGS (US Geological Survey). OFR 01177. [PDF / Author]
Powers, D.W. 2002a. Analysis report Task 1 of AP088, Construction of geologic contour maps. Carlsbad, NM: Sandia National Laboratories. ERMS 522085. [PDF / Author]
Powers, D.W. 2002b. Addendum to Analysis Report, Task 1 of AP088, Construction of Geologic Contour Maps. Carlsbad, NM: Sandia National Laboratories. ERMS 523886. [PDF / Author]
Powers, D.W. 2003. Addendum 2 to Analysis report Task 1 of AP088, Construction of geologic contour maps. Carlsbad, NM: Sandia National Laboratories. ERMS 522085. [PDF / Author]
Powers, D.W. 2006. Analysis Report Task 1B of AP114: Identify possible area of recharge to the Culebra west and south of WIPP. Carlsbad, NM: Sandia National Laboratories. ERMS 547094. [PDF / Author]
Powers, D.W. 2007a. Analysis Report for Task 1A of AP114: Refinement of Rustler Halite Margins Within the Culebra Modeling Domain. Carlsbad, NM: Sandia National Laboratories. ERMS 547559. [PDF / Author]
Powers, D.W. 2007b. Analysis Report for Task 1A of AP114: refinement of Rustler Halite margins within the Culebra modeling domain. Carlsbad, NM: Sandia National Laboratories. ERMS 547559. [Author]
Powers, D.W. 2009. Basic Data Report for Drillhole SNL8 (C3150) (Waste Isolation Pilot Plant). DOE/WIPP 053324. Carlsbad, NM: U.S. Department of Energy. [PDF / Author]
Powers, D.W. and R.M. Holt. 1995. Regional geological processes affecting Ruster hydrogeology. IT Corporation for Westinghouse Electric Corporation. [PDF / Author]
Powers, D.W. and R.M. Holt. 1999. The Los Medanos Member of the Permian Rustler Formation. New Mexico Geology. 21(4), pp.97103. [Author]
Powers, D.W. and R.M. Holt. 2000. The salt that wasn't there: mudflat facies equivalents to halite of the Permian Rustler Formation, southeastern New Mexico. Journal of Sedimentary Research, 70(1), pp.2936. [Author]
Powers, D.W., R.M. Holt, R.L. Beauheim, and S.A. McKenna. 2003. Geologic factors related to the transmissivity of the Culebra Dolomite Member, Permian Rustler Formation, Delaware Basin, southeastern New Mexico. In K.S. Johnson and J.T. Neal, eds. Evaporite Karst and engineering/environmental problems in the United States. 109th ed. Oklahoma Geological Survey. pp.21118. [Author]
Powers, D.W., R.M. Holt, R.L. Beauheim, and R.G. Richardson. 2006. Advances in depositional models of the Permian Rustler Formation, southeastern New Mexico. In Caves & Karst of Southeastern New Mexico. Fiftyseventh Annual Field Conference Guidebook ed. NM Geological Society. pp.26776. [Author]
Powers, D.W. and D. Owsley. 2003. A field survey of evaporite karst along NM 128 realighment routes. In K.S. Johnson and J.T. Neal, eds. Evaporite karst and engineering / environmental problems in the United States. 109th ed. Oklahoma Geological Survey. pp.23340. [Author]
Ramsey, J.L. and M. Wallace. 1996. Analysis Package for the Culebra Flow and Transport Calculations (Task 3) of the Performance Assessment Analyses Supporting the Compliance Certification Application. Carlsbad, NM: Sandia National Laboratories. ERMS 240516. [PDF / Author]
Roberts, R.M. 2006. Analysis Report for AP070, Analysis of Culebra Pumping Tests Performed Between December 2003 and August 2005. Carlsbad, NM: Sandia National Laboratories. ERMS 543901. [PDF / Author]
Roberts, R.M. 2007. Analysis Report for AP070, Analysis of Culebra Hydraulic Tests Performed Between June 2006 and September 2007. Carlsbad, NM: Sandia National Laboratories. ERMS 547418. [PDF / Author]
Rudeen, D.K. 2003. User's Manual for DTRKMF Version 1.00. Carlsbad, NM: Sandia National Laboratories. ERMS 523246. [PDF / Author]
Toll, N.J. and P.B. Johnson. 2006a. Routine Calculations Report In Support of Task 6 of AP114, SNL14 August 2005 Pumping Test Observation Well Data Processing, Summary of Files. Carlsbad, NM: Sandia National Laboratories. ERMS 543371. [PDF / Author]
Toll, N.J. and P.B. Johnson. 2006b. Routine Calculations Report In Support of Task 6 of AP114, WIPP11 February 2005 Pumping Test Observation Well Data Processing  Summary of Files. Carlsbad, NM: Sandia National Laboratories. ERMS 543651. [PDF / Author]
U.S. Department of Energy (DOE). 1996. Title 40 CFR Part 191 Compliance Certification Application for the Waste Isolation Pilot. Carlsbad, NM: U.S. Department of Energy Waste Isolation Pilot Plant, Carlsbad Area Office. DOE/CAO19962184. [Author]
U.S. Department of Energy (DOE). 2004. Title 40 CFR Part 191 Subparts B and C Compliance Recertification Application for the Waste Isolation Pilot Plant. Carlsbad, NM: US Department of Energy Carlsbad Field Office. DOE/WIPP 20043231. [Author]
U.S. Department of Energy (DOE). 2009. Title 40 CFR Part 191 Subparts B and C Compliance Recertification Application for the Waste Isolation Pilot Plant. Carlsbad, NM: US Department of Energy Carlsbad Field Office. DOE/WIPP 20093424.* [Author]
U.S. Environmental Protection Agency (EPA). 1996. 40 CFR Part 194: Criteria for the Certification and Recertification of the Waste Isolation Pilot Plant's Compliance with 40 CFR Part 191 Disposal Regulations; Final Rule., 1996. Federal Register Vol 61, 52235245. [PDF / Author]