NUREG/CR-6767

PNNL-13832 

Evaluation of Hydrologic Uncertainty Assessments for 

Decommissioning Sites Using Complex and Simplified 

Models

Prepared by

P. D. Meyer

S. Orr

T. J. Nicholson, NRC Project Manager

Pacific Northwest National Laboratory

Richland, WA 99352

March 2002

Prepared for

Division of Systems Analysis and Regulatory Effectiveness

Office of Nuclear Regulatory Research

U. S. Nuclear Regulatory Commission

Washington, DC 20555-0001

NRC JCN W6933

iii

Abstract

This report is the third (and final) report in a series that 

addresses issues related to hydrologic uncertainty assessment at decommissioning sites. Analyses in the first two 

reports in this series emphasized the application of relatively simplified models of subsurface flow and transport. 

Because of their relative computational speed, such simplified models are particularly attractive when the impact of 

uncertainty in flow and transport needs to be evaluated. 

These same simplifications, however, have the potential to 

provide unrepresentative estimates of dose and its uncertainty. Such misrepresentation may have important consequences for decisions based on the dose assessments. The 

significance of this concern was evaluated by comparing 

results from uncertainty assessments conducted on a test 

case using a simplified modeling approach and a more complex/realistic modeling approach.

The test case used a three-dimensional domain with a U-234 

source in the near surface, a 5-m-thick aquifer 7 m below 

the surface, and a small well pumping directly downstream, 

on the boundary of the contaminated zone. Exposure was 

assumed to occur through the drinking water pathway only, 

with all drinking water originating from the pumped well. A 

series of Monte Carlo simulations of flow and transport 

were performed using STOMP as the complex model and 

RESRAD as the simplified model. Hydraulic conductivity 

and air entry were modeled as random fields in STOMP 

using geostatistics from the Las Cruces Trench site. Random distributions for hydraulic parameters in the RESRAD 

simulations were based on both site-specific data and 

generic distributions. Sensitivity analyses were conducted 

on both codes using a combination of Monte Carlo simulation and single parameter variation. 

Peak doses predicted by the simplified model were several 

times higher than peak doses predicted by the complex 

model. This difference was attributed to the lack of dispersion in RESRAD and differences in aquifer mixing. The 

RESRAD concentration breakthrough curves exhibited a 

sharp peak with essentially no contaminant in the well until 

the time of the peak while STOMP predicted a much earlier 

arrival of contaminants in the well. Which code provided 

conservative results thus depended on whether the RESRAD peak occurred before the 1000-year regulatory criterion.

The random field characterization of the subsurface for the 

complex model used all available site data. Uncertainty in 

predicted dose was correspondingly small, with the peak 

dose coefficient of variation being 30%. When the variances 

of parameters in the simplified model were based on a 

generic dataset, the uncertainty in predicted peak dose was 

much larger; the coefficient of variation was 52% in this 

case. When the variances of parameters in the simplified 

model were based on the extensive site-specific data, the 

coefficient of variation for the peak dose was reduced to 

22%. In this case, however, the mean peak dose was actually less similar to the STOMP results than the generic case. 

For the RESRAD Monte Carlo simulations involving random soil hydraulic properties, the variability of peak dose 

was entirely attributed to variability in the aquifer hydraulic 

conductivity. Sensitivity to other parameters was examined 

by varying one parameter at a time. These results indicated 

that the recharge rate, the aquifer gradient, and the depth of 

penetration of the well were significant contributors to 

uncertainty in peak dose and the time of the peak dose. 

RESRAD predicted peak dose was more sensitive to the 

parameter values than was the STOMP predicted peak dose 

for the aquifer hydraulic conductivity, aquifer gradient, and 

depth of well penetration. Sensitivity to the recharge rate 

appeared to be comparable for the two codes.

Stochastic predictions of mean dose over time for the complex model were relatively insensitive to the geostatistical 

parameters. Ensemble mean peak dose predicted by the 

complex model was sensitive to the ensemble mean hydraulic conductivity. For a given ensemble mean hydraulic conductivity, however, the complex model showed no 

correlation between the spatial geometric mean aquifer conductivity of individual realizations and the resulting peak 

dose. This was in contrast to the simplified model in which 

the dose from individual realizations was strongly correlated with the aquifer conductivity. This result has implications for the value of hydraulic conductivity data. Adopting 

the homogeneous parameterization of RESRAD leads to a 

conclusion that reducing the uncertainty in the value of the 

aquifer hydraulic conductivity parameter will have a significant impact on the uncertainty in peak dose. Looking at the 

individual realization results from the STOMP model suggests that characterization of the average aquifer hydraulic 

conductivity is relatively unimportant in reducing uncertainty in peak dose. Characterizing the pattern of aquifer 

heterogeneity is likely to be more important than obtaining a 

value for the spatial mean hydraulic conductivity.

v

Contents

Abstract ...................................................................................................................................................................................... iii

Executive Summary .................................................................................................................................................................... xi

Foreword .................................................................................................................................................................................. xiii

Acknowledgments...................................................................................................................................................................... xv

1 Introduction ............................................................................................................................................................................... 1

1.1 Summary of NUREG/CR-6656........................................................................................................................................ 1

1.2 Summary of NUREG/CR-6695........................................................................................................................................ 4

1.3 Objectives and Methodology of the Current Report ........................................................................................................ 6

2 Test Case Scenario and Hydrogeologic Setting ........................................................................................................................ 9

2.1 Complex Representation .................................................................................................................................................. 9

2.1.1 Las Cruces Trench Data, Domain, and Geostatistics .............................................................................................. 9

2.1.2 Choice of Random Field Generator, Resulting Geostatistics of α and Ks and Spatial Distribution..................... 10

2.1.3 Flow and Transport Parameters ............................................................................................................................. 11

2.1.4 Numerical Modeling Considerations with STOMP: Grid Resolution, Time Step, and Solver............................. 12

2.1.5 Sample Results from Single Realizations.............................................................................................................. 16

2.2 Simplified Representation .............................................................................................................................................. 16

3 Base Case Monte Carlo Simulation Results............................................................................................................................ 25

3.1 STOMP Ensemble Mean Results ................................................................................................................................... 25

3.2 Dose Distributions and Comparison............................................................................................................................... 26

4 Sensitivity Cases ..................................................................................................................................................................... 31

4.1 Sensitivity to Geostatistical Parameters ......................................................................................................................... 31

4.1.1 Variance of Ks and ψe .......................................................................................................................................... 31

4.1.2 Correlation Scale ................................................................................................................................................... 32

4.2 Sensitivity to Other Hydrologic Parameters................................................................................................................... 32

4.2.1 Recharge Rate........................................................................................................................................................ 32

4.2.2 Mean Air Entry Pressure ....................................................................................................................................... 35

4.2.3 Mean Hydraulic Conductivity ............................................................................................................................... 35

4.2.4 Aquifer Hydraulic Gradient................................................................................................................................... 37

4.3 Sensitivity to Scenario Parameters ................................................................................................................................. 38

4.3.1 Depth of Well Penetration into the Aquifer........................................................................................................... 39

5 Conclusions ............................................................................................................................................................................. 43

6 References ............................................................................................................................................................................... 45

A Recommended Soil Parameter Distributions ....................................................................................................................... A-1

A.1 The Normal Distribution ............................................................................................................................................. A-1

A.2 The Lognormal distribution......................................................................................................................................... A-1
vi

A.3 The Beta Distribution .................................................................................................................................................. A-1

A.4 Recommended Probability Distributions for Soil Hydraulic Parameters by Soil Texture.......................................... A-2

A.5 References ................................................................................................................................................................... A-8

B Soil Parameter Correlation Coefficients............................................................................................................................... B-1

B.1 References.................................................................................................................................................................... B-7

C Recommended Soil Bulk Density Distributions................................................................................................................... C-1

C.1 References.................................................................................................................................................................... C-1

D Summary of Water Retention and Conductivity Models ..................................................................................................... D-1

D.1 Van Genuchten Model................................................................................................................................................. D-1

D.2 Brooks-Corey Model ................................................................................................................................................... D-1

D.3 Campbell Model .......................................................................................................................................................... D-1

D.3.1 Calculation of Campbell’s b Parameter .............................................................................................................. D-2

D.4 Additional Parameters ................................................................................................................................................. D-2

D.5 References ................................................................................................................................................................... D-2

Figures

1-1 NRC decommissioning and license termination decision-making framework (from NRC, 1998) ................................... 1

1-2 Probabilistic results from the RESRAD in situ test case simulations of NUREG/CR-6695 ............................................. 5

1-3 Prior and updated probabilistic results from the RESRAD excavation test case simulations of NUREG/CR-6695 ........ 6

2-1 The domain for the test case simulations ........................................................................................................................... 9

2-2 Contaminant source (U-234) concentrations in the solid phase (left) and the liquid phase (right) ................................. 10

2-3 Variograms (top) and histogram (bottom) of log(α) for a single realization with the analytical directional variograms 

shown for comparison (H=horizontal, V=vertical).......................................................................................................... 11

2-4 Variograms (top) and histogram (bottom) of log(Ks) (m/day) (single realization) ......................................................... 12

2-5 Permeability distribution of a single realization (Las Cruces geostatistics) .................................................................... 13

2-6 Histogram of vertical flux in a low-variance, random realization ................................................................................... 14

2-7 Comparison of concentrations at three locations in the domain using three different maximum time step sizes (300, 900, 

and 2600 years) for a random realization (#81) ............................................................................................................... 15

2-8 Concentration depth profiles at the (x, y) location of the well using a 1-m grid and a 0.5-m grid for a uniform domain after

500 years ...........................................................................................................................................................................16

2-9 Comparison between different weighting schemes for transport simulations at three locations in the domain .............. 17

2-10 Saturation (top) and pressure (bottom) distributions in a nonuniform (random) realization (Var[ln(Ks)] = 1.5) ........... 18

2-11 Total head distribution for the same realization as Figure 2-10....................................................................................... 19

2-12 Snapshots of the simulated U-234 plume including transport from the source, extraction from the well, and radioactive 

decay ................................................................................................................................................................................ 20
vii

2-13 Breakthrough curves from a single random realization at four locations in the domain and for the average of the five well 

nodes ................................................................................................................................................................................ 21

2-14 Vertical profiles of concentration at the (x, y) location of the well as calculated by STOMP for a randomly heterogeneous 

realization......................................................................................................................................................................... 21

2-15 Concentration in the well for the RESRAD Mass Balance (left) and Nondispersion (right) models.............................. 23

3-1 Top: Concentrations at the well with plus/minus one standard deviation, maximum well concentration, and concentration 

in the unsaturated zone immediately above the well, Bottom: Concentration on a log scale for the five nodes comprising 

the well (numbered from bottom to top of the aquifer), the node just above the well in the unsaturated zone, a downstream 

node, and the mean well concentration (bottom) ............................................................................................................. 26

3-2 Variance of well concentration at 579 years after start of release, as a function of the number of simulations.............. 27

3-3 Ensemble mean concentration from the STOMP Base Case Monte Carlo simulation compared to the deterministic concentration using equivalent uniform (scaled mean) parameter values ............................................................................. 27

3-4 Empirical cumulative distribution functions for the peak dose (top) and time of peak (bottom) over 300 simulations for 

the complex model using STOMP and the simplified model using RESRAD................................................................ 28

3-5 Mean dose as a function of time for the complex model using STOMP and the simplified model using RESRAD...... 30

4-1 Ensemble mean concentration in the well (solid lines) and plus/minus one standard deviation about the mean well concentration (dashed lines) for the base case and a Monte Carlo simulation consisting of 300 realizations with double the 

variance of log[Ks] and ψe............................................................................................................................................... 31

4-2 Evolution of the variance of the well concentration at 579 years with the number of Monte Carlo realizations - variance 

of log[Ks] two times the base case value ......................................................................................................................... 32

4-3 Effect of correlation scale on breakthrough concentrations at the pumping well............................................................ 33

4-4 Effect of recharge on breakthrough concentrations at the pumping well ........................................................................ 34

4-5 Effect of the recharge rate on peak dose and on the time of the peak dose for the simplified (RESRAD) and the complex 

(STOMP) models ............................................................................................................................................................. 35

4-6 Effect of air-entry value on breakthrough curves at the unsaturated node above the well, at the top well node, and downstream mid-depth in the aquifer ....................................................................................................................................... 36

4-7 Effect of mean saturated conductivity on the ensemble mean well concentration for 200 realization Monte Carlo simulations .................................................................................................................................................................................. 36

4-8 RESRAD sensitivity to aquifer hydraulic conductivity and unsaturated zone bulk density from the base case simulation 

(top) and a comparison of the RESRAD (simplified) and STOMP (complex) results for peak dose sensitivity to aquifer 

hydraulic conductivity (STOMP results from Figure 4-7)............................................................................................... 38

4-9 STOMP sensitivity to the spatially averaged (geometric mean) aquifer hydraulic conductivity using the high recharge 

case................................................................................................................................................................................... 39

4-10 Well concentration over time for aquifer hydraulic gradients of 0.002 to 0.02 from the complex model (STOMP) using 

equivalent uniform hydraulic properties (top) and a comparison of peak dose for the complex and simplified models as a 

function of aquifer hydraulic gradient using equivalent uniform properties in both cases (bottom) .............................. 40

4-11 STOMP model base case results (top) with a fully penetrating well (5-m well depth) and results from a simulation with 

a 3-m deep well (bottom) ................................................................................................................................................. 41

4-12 Well concentration over time for well depths of 1 m to 5 m from the complex model (STOMP) using equivalent uniform 

hydraulic properties (top) and a comparison of peak dose for the complex and simplified models as a function of well 

depth using equivalent uniform properties in both cases (bottom) .................................................................................. 42
viii

Tables

1-1 Summary of radiological criteria for license termination (10 CFR Part 20 Subpart E)..................................................... 2

1-2 Generalized relative importance of dose assessment code parameters in an uncertainty analysis ................................... 3

1-3 Summary of test case deterministic base case simulation results from NUREG/CR-6695 ............................................... 4

2-1 Statistics of log(α) for a single realization ...................................................................................................................... 13

2-2 Statistics of log(Ks) and Ks (m/day) for a single realization ........................................................................................... 13

2-3 Model parameters used in the STOMP simulations......................................................................................................... 14

2-4 Statistics of vertical fluxes (cm/hr) in uniform and random realizations (Las Cruces statistics)..................................... 14

2-5 Generic stochastic parameter inputs used in RESRAD simulations - default distributions from NUREG/CR-6565 for a 

sand soil............................................................................................................................................................................ 22

2-6 Site-specific stochastic parameter distributions used in the RESRAD simulations - based on the Las Cruces Trench 

dataset where applicable .................................................................................................................................................. 23

3-1 Statistics of peak dose and the time of the peak dose ...................................................................................................... 29

A-1 Recommended distributions for Sand ............................................................................................................................ A-2

A-2 Recommended distributions for Loamy Sand................................................................................................................ A-2

A-3 Recommended distributions for Sandy Loam................................................................................................................ A-3

A-4 Recommended distributions for Sandy Clay Loam ....................................................................................................... A-3

A-5 Recommended distributions for Loam........................................................................................................................... A-4

A-6 Recommended distributions for Silt Loam .................................................................................................................... A-4

A-7 Recommended distributions for Silt............................................................................................................................... A-5

A-8 Recommended distributions for Clay Loam .................................................................................................................. A-5

A-9 Recommended distributions for Silty Clay Loam.......................................................................................................... A-6

A-10 Recommended distributions for Sandy Clay.................................................................................................................. A-6

A-11 Recommended distributions for Silty Clay .................................................................................................................... A-7

A-12 Recommended distributions for Clay............................................................................................................................. A-7

B-1 Correlation coefficients for Sand ................................................................................................................................... B-1

B-2 Correlation coefficients for Loamy Sand ....................................................................................................................... B-2

B-3 Correlation coefficients for Sandy Loam ....................................................................................................................... B-2

B-4 Correlation coefficients for Sandy Clay Loam............................................................................................................... B-3

B-5 Correlation coefficients for Loam .................................................................................................................................. B-3

B-6 Correlation coefficients for Silt Loam............................................................................................................................ B-4

B-7 Correlation coefficients for Silt...................................................................................................................................... B-4

B-8 Correlation coefficients for Clay Loam.......................................................................................................................... B-5

B-9 Correlation coefficients for Silty Clay Loam ................................................................................................................. B-5

B-10 Correlation coefficients for Sandy Clay......................................................................................................................... B-6
ix

B-11 Correlation coefficients for Silty Clay ........................................................................................................................... B-6

B-12 Correlation coefficients for Clay.................................................................................................................................... B-7

C-1 Recommended parameters of normal distributions for bulk density ............................................................................. C-1

xi

Executive Summary

Assessments of the dose from contaminated sites and waste 

disposal facilities may rely on the use of relatively simplified models of subsurface flow and transport. Common simplifications include steady state, one-dimensional flow; 

homogeneous and isotropic transport medium properties; 

and unit hydraulic gradient or piston-flow displacement in 

the unsaturated zone. Because of their relative computational speed, such simplified models are particularly attractive when the impact of uncertainty in flow and transport 

needs to be evaluated. Simplifications in the representation 

of flow and transport have the potential to provide unrepresentative estimates of dose and its uncertainty. “Unrepresentative” is used here to describe a dose or uncertainty 

estimate that significantly misrepresents the actual dose or 

the uncertainty associated with it. Such misrepresentation 

may have important consequences for decisions based on 

the dose assessments. The significance of this concern was 

evaluated here by comparing results from uncertainty 

assessments conducted on a test case using a simplified 

modeling approach and a more complex/realistic modeling 

approach. The test case followed the U.S. Nuclear Regulatory Commission’s framework for site decommissioning 

analyses. Subsurface properties were derived from data 

obtained in the Las Cruces Trench experiments with source 

term data reflecting an actual decommissioning case. Comparisons between the two approaches included the probability distribution of peak dose, the relative importance of 

parameters, and the value of site-specific data in reducing 

uncertainty.

Characteristics of the test case were set to eliminate many 

differences between the codes and to highlight the differences due to the underlying conceptualization and parameterization of subsurface flow and transport. We considered a 

three-dimensional domain that resembled the subsurface at 

the Las Cruces Trench site – a site with abundant information on soil properties and information on their geostatistics. 

We assumed a decaying source of U-234 in the near surface, 

a 5-m-thick aquifer 7 m below the surface, and a small well 

pumping directly downstream on the boundary of the contaminated zone. Exposure was assumed to occur through the 

drinking water pathway only, with all drinking water originating from the pumped well. In order to evaluate doses at 

the pumping well over 1000 years, a series of Monte Carlo 

simulations of flow and transport were performed using the 

STOMP code for the complex model and RESRAD for the 

simplified model. The generated random fields for the 

STOMP model were based on the geostatistics of Las 

Cruces. Distributions of random parameters in the RESRAD simulations were based on both site-specific data and 

generic distributions.

Peak doses predicted by the simplified model (RESRAD) 

were several times higher than peak doses predicted by the 

complex model (STOMP). This difference was attributed to 

the lack of dispersion in the RESRAD model and differences in the mixing and dilution in the aquifer. For this test 

case, RESRAD provided conservative results in the sense 

that the peak dose was higher than that predicted by 

STOMP. The RESRAD concentration breakthrough curves 

exhibited a sharp peak, however, with essentially no contaminant in the well until the time of the peak. This could be 

interpreted as nonconservative behavior since the STOMP 

simulations predicted a much earlier arrival of contaminants 

and, prior to the arrival of the RESRAD peak, the well concentrations from the STOMP model were larger than those 

predicted by RESRAD. This point is particularly relevant 

because the regulatory criteria include a time limit: 1000 

years from the time of decommissioning. Although the 

RESRAD predicted peak occurred prior to this limit, had 

the peak occurred after 1000 years, the RESRAD results 

would have been nonconservative in every sense of the 

word.

The complex model assumed the availability of an extensive 

dataset on which to base the random field characterization 

of the subsurface. Uncertainty in predicted dose was correspondingly small, with the peak dose coefficient of variation being 30%. When the variances of parameters in the 

simplified model were based on a generic dataset, the uncertainty in predicted peak dose was much larger; the coefficient of variation was 52% in this case. When the variances 

of parameters in the simplified model were based on the 

available site-specific data, the coefficient of variation for 

the peak dose was reduced to 22%. In this case, however, 

the mean peak dose was actually less similar to the STOMP 

results than the generic case. This unexpected result 

occurred because the site-specific mean aquifer hydraulic 

conductivity was significantly smaller than the generic 

value and because of the strong dependence of the RESRAD peak dose on aquifer hydraulic conductivity. The lognormal distribution of the hydraulic conductivity 

accentuated the effect.

Sensitivity analyses were conducted on both codes using a 

combination of Monte Carlo simulation and single parameter variation. The goals of the sensitivity studies were to 

determine the parameters critical to an understanding of the 

models’ behavior. We also wanted to see whether the simplified and complex models were sensitive to the same 

parameters, thus indicating whether they represent the modeled flow and transport processes in a physically similar 

way.

For the RESRAD Monte Carlo simulations involving random soil hydraulic properties, the variability of peak dose 

was entirely attributed to variability in the aquifer hydraulic 

conductivity. Sensitivity to other parameters was examined 

by varying one parameter at a time. These results indicated 

that the recharge rate, the aquifer gradient, and the depth of 

penetration of the well were significant contributors to 

uncertainty in peak dose and the time of the peak dose. 

Sensitivity analyses of the complex model showed that stochastic predictions of mean dose over time exhibited low 
xii

sensitivity to parameter variability (expressed by their variances), including to the variance of log hydraulic conductivity. This may have been due to the use of a relatively small 

simulation domain. Uncertainty in dose predictions, however, was directly (linearly) dependent on the variance of 

hydraulic conductivity. Ensemble mean dose was insensitive to the correlation scales of the hydraulic properties. 

Uncertainty in dose, however, increased somewhat with an 

increase in the correlation scale. Overall, the complex 

model showed robustness with respect to geostatistical 

parameters.

Ensemble mean peak dose predicted by the complex model 

was strongly sensitive to the ensemble mean hydraulic conductivity. Reducing the ensemble mean hydraulic conductivity by an order of magnitude resulted in more than a fivefold increase in the peak dose (due primarily to reduced 

mixing in the aquifer). Increasing the ensemble mean conductivity by an order of magnitude resulted in nearly a 90% 

reduction in peak dose. It was noted, however, that for a 

given ensemble mean hydraulic conductivity, the complex 

model showed no correlation between the spatial geometric 

mean aquifer conductivity of individual realizations and the 

resulting peak dose. This was in contrast to the simplified 

model in which the dose from individual realizations was 

strongly (negatively) correlated with the aquifer conductivity. This result has implications for the value of hydraulic 

conductivity data. Adopting the homogeneous parameterization of RESRAD leads to a conclusion that reducing the 

uncertainty in the value of the aquifer hydraulic conductivity parameter will have a significant impact on the uncertainty in peak dose. The individual realization results from 

the STOMP model suggest that characterization of the average aquifer hydraulic conductivity is relatively unimportant 

in reducing uncertainty in peak dose. Characterizing the pattern of aquifer heterogeneity is likely to be more important 

than obtaining a value for the spatial mean hydraulic conductivity.

RESRAD predicted peak dose was more sensitive to the 

parameter values than was the STOMP predicted peak dose 

for the aquifer hydraulic conductivity, aquifer gradient, and 

depth of well penetration. Sensitivity to the recharge rate 

appeared to be comparable for the two codes, perhaps a 

result of the boundary condition constraints imposed on the 

STOMP model. Differences in model sensitivities may lead 

to different conclusions about which parameters are critical 

to the analysis and about the importance of site-specific 

data.

Finally, we make some observations on computational 

issues. Flow and transport in single realizations of the complex model showed strong dependence on the spatial structure of generated hydraulic properties. These large 

variations in results, however, diminished sharply when 

averaged over the ensemble, as the number of realizations 

exceeded 200. A minimum of 200 realizations were thus 

conducted for each of the STOMP Monte Carlo simulations. 

STOMP model results showed some sensitivity to mesh resolution, particularly near the water table, where concentration gradients were extremely high. The relatively coarse 

mesh size (1 m) was a calculated compromise without 

which the Monte Carlo simulations could not be performed. 

The variety of numerical considerations required by the 

application of a complex, three-dimensional model was 

clearly illustrated. No such considerations were required to 

execute RESRAD as its internal conceptual-mathematical 

model is fixed. The complex model required many hours of 

computational time to execute a Monte Carlo simulation. 

RESRAD, on the other hand, required no more than a few 

minutes.
xiii

Foreword

This technical contractor report was prepared by Pacific Northwest National Laboratory1 (PNNL) under their U. S. Department of Energy Interagency Work Order (JCN W6933) with the Radiation Protection, Environmental Risk and Waste Management Branch, Division of Systems Analysis and Regulatory Effectiveness, Office of Nuclear Regulatory Research. This final 

research report documents the testing and completion of PNNL's uncertainty assessment methodology. The testing utilized 

field datasets from complex ground-water flow and tracer experiments for comparison of alternative models and parameter 

representations. The field datasets were chosen in consultation with the NRC Project Manager and NRC licensing staff. Earlier 

application of the methodology involved a hypothetical test case provided by NRC licensing staff (i.e., NUREG/CR-6695). 

The PNNL methodology provides a strategy for identifying the critical hydrologic parameters and evaluating their contribution to uncertainty in dose calculations. Results from this work point to the importance of both parameter uncertainty and conceptual model uncertainty. The report's appendices provide a listing of the data distributions used in the testing. Further 

detailed information on the PNNL research study and a listing of PNNL's Simplified Transient Estimation of the Water Budget 

(STEWB) Model is available on the PNNL website: http://nrc-hydro-uncert.pnl.gov/.

The PNNL research study was undertaken in response to an Office of Nuclear Material Safety and Safeguards “user need” to 

support licensing needs for estimating and reviewing hydrologic parameter distributions and their attendant uncertainties for 

site-specific dose assessment modeling as outlined in NUREG-1549. As discussed in NUREG-1549, the NRC staff reviews 

the scenario, conceptual model, and input parameters. All three components may contribute to uncertainty. The PNNL 

research focuses on hydrologic parameter uncertainties in the context of dose assessments for decommissioning sites. The 

information provided in this and earlier reports (i.e., NUREG/CR's 6656 and 6695) supports the NRC staff's efforts in providing dose modeling guidance. Specifically, the report documents the methodology, its relationship to conceptual model uncertainty, and demonstrates its capabilities through the testing exercises which utilized site-specific data.

This final report is the third report in a series of three contractor reports documenting PNNL's uncertainty assessment methodology, its testing and applications to decommissioning sites. This report, as with the earlier ones, is not a substitute for NRC 

regulations, and compliance is not required. The approaches and/or methods described in this NUREG/CR are provided for 

information only. Publication of this report does not necessarily constitute NRC approval or agreement with the information 

contained herein. Use of product or trade names is for identification purposes only and does not constitute endorsement by the 

NRC or Pacific Northwest National Laboratory.

Cheryl A. Trottier, Chief

Radiation Protection, Environmental Risk and Waste Management Branch

Division of Systems Analysis and Regulatory Effectiveness

Office of Nuclear Regulatory Research

1 Pacific Northwest National Laboratory is operated for the U.S. Department of Energy by Battelle Memorial Institute under Contract DE-AC06-76RLO 1830

xv

Acknowledgments

The authors gratefully acknowledge the financial assistance 

of the Office of Nuclear Regulatory Research of the U.S. 

Nuclear Regulatory Commission (NRC) and the guidance 

of the NRC Project Manager, Tom Nicholson. The staff of 

the NRC have also contributed in reviewing and commenting on the activities described in this report.

Joel Malard and Mark White of PNNL provided significant 

help in performing the STOMP simulations. We also appreciate the comments of Dongxiao Zhang of Los Alamos 

National Laboratory who provided a review of this report.

1

1 Introduction

This report is the third report in a series that addresses 

issues related to hydrologic uncertainty assessment at 

decommissioning sites. The analyses described in these 

reports have been conducted under the general framework 

for license termination described in NUREG-1549 (NRC, 

1998). This framework provides for an iterative approach to 

decommissioning analyses that moves from screening-type 

analyses using pre-defined models and generic screening 

parameters to site-specific analyses using more complex, 

realistic models and parameters based on site-specific data. 

A depiction of the framework taken from NUREG-1549 is 

shown in Figure 1-1. 

A primary component of the decommissioning framework 

is a dose assessment (box 4 of Figure 1-1), used to provide 

future estimates of dose that can be compared to the regulatory criteria for license termination (10 CFR Part 20 Subpart 

E). The regulatory criteria are summarized in Table 1-1 and 

involve predictions of contaminant transport and exposure 

1000 years from the time of decommissioning. This long 

time frame, coupled with complex contaminant exposure 

pathways and limitations on site-specific characterization 

data produce an inherent uncertainty in the results of a dose 

assessment. 

Each of the reports in this series addresses a particular 

aspect of dose assessment uncertainty analysis. To provide a 

complete picture of the analyses completed, the first two 

reports are summarized here.

1.1 Summary of NUREG/CR-6656

The first report in this series, NUREG/CR-6656 (Meyer and 

Gee, 1999), described the basic conceptual models and 

mathematical implementations of three dose assessment 

codes (DandD, RESRAD, and MEPAS1

). Uncertainty in the 

hydrologic parameters of these codes was discussed and 

Figure 1-1. NRC decommissioning and license termination decision-making framework (from NRC, 1998)

1 DandD conceptual model and mathematical formulations are described 

in Kennedy and Strenge (1992)

RESRAD web page: http://web.ead.anl.gov/resrad/home2

MEPAS web page: http://mepas.pnl.gov:2080/earth/earth.htm
Introduction

2

data sources for defining best-estimate parameter values and 

parameter uncertainty characterizations were identified. The 

three codes were selected because they were in use by NRC 

staff as part of their decommissioning analyses. The NRC 

does not require the use of these codes, however. 

Each of the codes uses a fairly simplistic representation of 

flow and transport in the subsurface. The common simplifying conditions used in DandD, RESRAD, and MEPAS 

related to flow and transport in the subsurface were summarized as follows.

• Each code uses a relatively simple model for the nearsurface water budget to determine the net infiltration 

rate, i.e., the flow of water through the system that drives 

the contaminant transport. In DandD the net infiltration 

is known, in RESRAD the water budget is based on 

average annual components (precipitation, evapotranspiration, and runoff), and for MEPAS the water budget is 

based on monthly average components.

• Steady-state flow throughout the system is imposed in 

each code. The net infiltration required by each code is a 

constant, average annual value. 

• Flow throughout the system is one-dimensional. In addition, advective transport is also one-dimensional in all 

cases. (MEPAS includes the effect of dispersion, which 

can be three-dimensional in the saturated zone.)

• Only a small number of layers are allowed in the system 

and the physical, hydraulic, and chemical properties of 

each layer must be uniform. No other form of heterogeneity is allowed.

• A fairly simple model of mixing in the aquifer is used. In 

DandD mixing occurs instantaneously throughout the 

fixed aquifer volume. For RESRAD and MEPAS the 

aquifer has infinite lateral extent and a finite and constant thickness.

As a consequence of these simplifications, there are sitespecific conditions under which the codes may provide 

inaccurate, potentially nonconservative results. These conditions include the following.

• The presence of significant preferential flow in the near 

surface could lead to greater net infiltration than might 

be anticipated when relying on the simple water budget 

calculations of the dose codes. 

• The imposed condition of steady-state flow will be violated to some degree at all sites. Net infiltration and percolation to the saturated zone vary in time as a result of 

intermittent precipitation and seasonal changes in evapotranspiration and precipitation. The contaminant flux 

may be affected as a result of temporal variation in net 

infiltration and water content. In addition, transient flow 

conditions may affect transport in the unsaturated zone. 

• Significant heterogeneity in physical, hydraulic, or 

chemical properties may result in two- or three-dimensional flow patterns. Contrasting properties that result in 

focused flow and fast transport paths may produce conditions under which the dose codes predict nonconservative results.

• Subsurface conditions under which flow and transport 

are significantly influenced by fractures or karst formations clearly contradict the simplified conceptual models 

of the dose assessment codes. Because of the potential 

for fast transport of contaminants under these conditions, the codes may produce nonconservative results.

• The presence of unmodeled processes such as colloidal 

transport, nonequilibrium adsorption, and vapor phase 

transport may render the model predictions nonconservative.

The hydrologic parameters of the three codes were listed 

and compared. The relative importance of the parameters to 

the analysis of uncertainty in predicted dose was summarized qualitatively as shown here in Table 1-2. The relative 

importance is composed of the sensitivity of the model to a 

parameter’s value and the uncertainty in the parameter 

value. Sources of uncertainty were discussed and alternative 

measures of sensitivity were described. 

Data sources for best-estimate parameter values and parameter uncertainty information were reviewed. The following 

sources of information were identified. 

Table 1-1. Summary of radiological criteria for license termination (10 CFR Part 20 Subpart E)

Unrestricted Releasea Restricted Release

Dose Criterion 25 mrem TEDE per year peak 

annual dose to the average 

member of the critical group

25 mrem TEDE per year peak 

annual dose to the average 

member of the critical group 

while controls are in place

100 or 500 mrem TEDE per year 

peak annual dose to the average 

member of the critical group 

upon failure of the controls

Time Frame 1000 years 1000 years 1000 years

Other 

Requirements

ALARA ALARA, financial assurance, 

public participation

ALARA, financial assurance, 

public participation

a. TEDE: total effective dose equivalent; ALARA: as low as reasonably achievable
Introduction

3

• UNSODA (Leij et al., 1996) is a database of soil physical and hydraulic property measurements compiled from 

a variety of international sources. UNSODA is available 

from the U. S. Department of Agriculture-Agricultural 

Resarch Service (USDA-ARS) Salinity Laboratory 

(http://www.ussl.ars.usda.gov/). UNSODA contains soil 

property measurements, but some analysis may be 

required to obtain soil parameters for use in dose assessments. Two related databases, compiled from U.S. 

sources, are summarized in Schaap and Leij (1998). 

• In NUREG/CR-6565 (Meyer et al., 1997) generic probability distributions for soil parameters were presented for 

use directly in RESRAD and MEPAS. These parameter 

distributions were based on the work of Carsel and Parrish (1988) who used data from the National Soil Characterization Database. These distributions are included 

in this report as Appendices A and B. Soil bulk density 

distributions obtained using a similar analysis are 

included here as Appendix C.

• The National Soil Characterization Database is a large 

collection of soil properties available from the USDA 

Natural Resources Conservation Service (NRCS) (see 

http://www.nrcs.usda.gov/). This database consists of 

standard soil survey measurements on small-scale soil 

pedons generally obtained near the surface. Some analysis may be required to derive parameter values from the 

measured soil properties. 

• SSURGO is a database of the original NRCS soil surveys digitized for use with Geographic Information System (GIS) software. These are the most detailed soil 

survey data available from the NRCS and are currently 

available for limited areas of the United States. 

SSURGO data are at a sufficiently small scale to be useful for parameter estimation in decommissioning analyses. Additional analyses may be required to derive 

parameter values from the soil property measurements in 

the database. More information and data are available at 

http://www.nrcs.usda.gov/.

• The STATSGO database is an NRCS product that generalizes the detailed soil survey data to a larger scale. This 

database is also intended for use with GIS software. The 

scale of STATSGO data is generally too large to be of 

significant use in estimating parameter values for 

decommissioning analyses, although it may be useful in 

characterizing uncertainty. 

• The National Climatic Data Center (NCDC) provides 

meteorological data measured at locations throughout 

the United States. Many of these data are available at 

http://lwf.ncdc.noaa.gov/oa/ncdc.html, although it may 

require some analysis to derive the parameter values 

required for dose assessment analyses. The use of this 

data in evaluating uncertainty in net infiltration estimates was described by Meyer et al. (1997).

• The NWIS is a database of surface water, groundwater, 

and water quality information available from the U.S. 

Geological Survey (USGS) at 

http://water.usgs.gov/nwis/. The groundwater information is primarily limited to water table levels in selected 

wells.

• Extension Service agents; state, county, and municipal 

staff; and university or industry personnel with experience in local conditions may be excellent sources of 

information and data related to a site. Information 

related to specific soil types, land use patterns, and climatic data are often readily available from extension 

offices and other local sources. Water table depths and 

results of percolation tests (which give some indication 

of soil hydraulic properties) can be obtained in areas 

where septic tank systems are used. Local county or city 

files are often the source for such information.

Additional information on data sources for the RESRAD 

and MEPAS dose assessment codes can be found in Yu et al. 

(1993) and Buck et al. (1995), respectively.

Taken as a whole, NUREG/CR-6656 outlines a methodology for parameter uncertainty assessment that considers the 

potential data limitations and modeling needs of decommisTable 1-2. Generalized relative importance of dose assessment code parameters in an uncertainty analysis from 

NUREG/CR-6656. Results for a specific application may vary.

Uncertainty Due to Variability and/or Lack of Knowledge

Low Medium High

Model Sensitivity

High UZa Thickness Distribution Coefficients

Net Infiltration Rate

Medium Effective Porosity

Bulk Density

Darcy Velocity

Unsaturated Water Content

SZa

 Exposure Parameters

SZ Hydraulic Conductivity

Low Porosity Soil-Type Exponent

Field Capacity

UZ Sat. Hydraulic Conductivity

Dispersivity

a. UZ = unsaturated zone, SZ = saturated zone
Introduction

4

sioning analyses. This methodology uses generic parameter 

distributions based on national or regional databases, sensitivity analysis, probabilistic modeling, and Bayesian updating to incorporate site-specific information. The uncertainty 

assessment methodology described in NUREG/CR-6656 is 

consistent with the technical basis for dose modeling evaluations prepared by NRC staff and presented in NUREG1727 (NRC, 2000).

1.2 Summary of NUREG/CR-6695

The second report in the series, NUREG/CR-6695 (Meyer 

and Taira, 2001) used hypothetical decommissioning test 

cases to illustrate the application of the uncertainty assessment methodology presented in NUREG/CR-6656. The test 

cases were based on source term and scenario information 

provided by NRC staff and on the physical setting of a site 

in Arizona at which NRC-sponsored field studies have been 

carried out. Basic soil and climate information provided by 

University of Arizona staff were used in the application. 

The site consisted of a layered sand/loam soil with a water 

table 13 m below the ground surface. The climate was arid, 

but it was assumed that substantial irrigation occurred. 

Other regional information was obtained from electronic 

sources. For those aspects of the site without reliable data 

sources, national databases were used to estimate site characteristics. 

A series of deterministic simulations were carried out using 

the codes DandD v. 1.0 and RESRAD v. 6.0. Simplifications to the conceptual model of the site were made to 

match the conceptual models embodied in the simulation 

codes. Following the framework described in NUREG1549, a DandD screening simulation was executed with the 

test case source term and all default parameter values. The 

DandD code was subsequently run with parameter values 

more representative of the test case site. The RESRAD simulations used primarily site-specific parameter values for 

the subsurface physical, hydrologic, and chemical parameters. Two source scenarios were modeled with RESRAD. 

An in situ case left the waste in its original buried location 

and assumed that a cover was in place. An excavation case 

was also simulated in which the waste was assumed to have 

been excavated for construction of a house, mixed with 

clean soil from the excavation, and widely spread about in a 

surface layer. This case more closely resembled the DandD 

conceptual model. The resulting peak doses over a 1000-

year period for these deterministic DandD and RESRAD 

simulations are shown in Table 1-3. 

In all cases, the primary contributors to the peak dose were 

the uranium isotopes. For the RESRAD simulations, the 

uranium release occurred over a relatively short period of 

time producing a correspondingly sharp peak in dose. 

Because of this sharp peak, lengthening the travel time 

slightly (less than 10% in the excavation test case) would 

produce a result that satisfies the regulatory criterion by 

moving the peak dose pulse to a point beyond 1000 years. 

Such a result demonstrates the importance of considering 

the uncertainty in parameter values. This result also emphasizes the importance of considering uncertainty in the time 

at which the peak dose occurs, in addition to uncertainty in 

the magnitude of the peak. 

Deterministic sensitivity analyses applied to the RESRAD 

in situ test case included use of a simplified model implemented in a spreadsheet and one-parameter sensitivity calculations applied to the base case parameter values and to a 

set of conservative parameter values. The various sensitivity 

measures were largely in agreement. For the in situ case, 

these analyses indicated that the evapotranspiration coefficient, the uranium distribution coefficients, and the well 

pumping rate were the most important parameters contributing to the uncertainty in peak dose. Soil hydraulic parameters were much less important for this case. Deterministic 

sensitivity analyses were not carried out for the excavation 

case.

Probabilistic analyses were carried out for the in situ and 

excavation cases using the Monte Carlo simulation capability of RESRAD. Histograms and cumulative distributions 

for the peak total dose and the time of the peak dose were 

derived. Figure 1-2(A) shows the empirical cumulative disTable 1-3. Summary of test case deterministic base case simulation results from NUREG/CR-6695

Code Description

Peak Dose 

(mrem/yr)

 Peak Dose 

Time (yr)

DandD Default Parameter Values 829 4

DandD Site-Specific Physical & Hydrologic Parameters 285 17

DandD Site-Specific Distribution Coefficients 198 10

DandD All Site-Specific Parameters 70 98

RESRAD In Situ Case 115 939

RESRAD Excavation Case 16 973
Introduction

5

tribution function for the RESRAD in situ test case. Statistics for the total dose as a function of time, used to obtain 

the peak of the mean dose, were also presented. A sample of 

these results is shown in Figure 1-2(B) for the RESRAD in 

situ test case. The probabilistic results can be used to compare the estimated site performance to the regulatory measures with consideration of parameter uncertainties. 

Probabilistic measures of sensitivity presented were scatter 

plots of peak dose versus parameter values, statistical sensitivity measures calculated by RESRAD (such as the partial 

correlation coefficients), and single-parameter Monte Carlo 

simulations used to clarify the relationships between dose 

and critical parameter values. No single measure was a reliable indicator of the relative importance of the parameters. 

For the in situ case, the evapotranspiration coefficient, the 

well pumping rate, and the uranium distribution coefficients 

were of greatest importance. These results were consistent 

with the deterministic results. For the excavation case, the 

saturated zone hydraulic conductivity was the most important parameter followed by the uranium distribution coefficients. The applications illustrate the value of applying a 

variety of uncertainty analysis methods and understanding 

the behavior of the simulation code. 

There are multiple measures of dose that could be used for 

comparison to the regulatory standard when probabilistic 

simulations are used in decommissioning analyses. One 

could use a statistic (such as the mean or the 50th percentile) of the distribution of peak dose. Peak dose as used here 

refers to the total (all pathways) dose occurring in the first 

1000 years after decommissioning. This is the quantity used 

in Figure 1-2(A) (albeit without the 1000 year limit). An 

alternative measure of dose discussed in NRC (2000) is the 

peak of the mean dose, where the mean dose is calculated as 

a function of time by averaging over all Monte Carlo realizations for each year of the simulation as shown in 

Figure 1-2(B). The median and 90th percentile doses are 

shown in this figure as well as the deterministic base case 

dose. For the in situ test case, the mean dose at any time 

during the simulation was significantly less than the base 

case peak dose. This result occurred because of the sharp 

peak in the dose as a function of time for any given realization and the fact that the parameter variation resulted in the 

time of the peak varying significantly. At any particular 

time, in fact, the fraction of realizations that had a dose 

greater than the 25 mrem/yr limit was quite small. The differences in the statistical quantities that might be used to 

compare with the regulatory criteria require that these quantities be precisely defined.

A method to update parameter probability distributions 

(based on Meyer et al., 1997) was also applied to the excavation test case using the saturated zone hydraulic conductivity. Because site-specific measurements were not 

available for this parameter, data were generated using four 

measurements of physical properties from the deepest samples available at the test case site. Updating the saturated 

zone hydraulic conductivity distribution with this generated 

data had its greatest effect on the standard deviation of the 

peak dose, which was significantly reduced. The empirical 

cumulative distribution functions for the prior and updated 

peak dose (for the excavation case) are shown in Figure 1-3. 

The percentage of realizations exceeding 25 mrem/yr was 

reduced from 13% to 2% for this test case. 

1.0

0.8

0.6

0.4

0.2

Empirical cdf

0 100 200 300 400

Peak Dose (mrem/yr)

N = 100

Min = 0

Max = 488.7

Mean = 120.48

StDev = 78.317

Median = 115.35

100

80

60

40

20

0

Total Dose (mrem/yr)

0 400 800 1200

Time (yr)

 Basecase

 90th Percentile

 Mean

 Median

(A) (B)

Figure 1-2. Probabilistic results from the RESRAD in situ test case simulations of NUREG/CR-6695
Introduction

6

1.3 Objectives and Methodology of the 

Current Report

Analyses in the first two reports in this series emphasized 

the application of relatively simplified models of subsurface 

flow and transport. Common simplifications of the codes 

examined in the reports include steady-state, one-dimensional flow, homogeneous and isotropic transport medium 

properties, and unit hydraulic gradient in the unsaturated 

zone. While recognizing the potential impact of these simplifications, the simplified codes are nonetheless attractive 

because of their relative computational speed. This speed is 

particularly important when the impact of uncertainty in 

flow and transport needs to be evaluated. 

Simplifications in the representation of flow and transport 

have the potential to result in not only unrepresentative estimates of dose, but also unrepresentative estimates of uncertainty. “Unrepresentative” is used here to describe an 

estimate of expected values (or mean values) and their 

uncertainty that significantly misrepresent the “actual” 

mean and uncertainty. Such misrepresentation may have 

important consequences for decisions based on the dose 

assessment, including conclusions drawn about the relative 

importance of parameters and the value of site-specific data 

in reducing uncertainty. The significance of this concern is 

evaluated in this report by comparing test case results from 

uncertainty assessments conducted using a simplified, onedimensional, statistical modeling approach and a complex, 

three-dimensional, stochastic modeling approach. We use 

the results from the complex, three-dimensional stochastic 

model as our best estimate of the “actual” mean and uncertainty. Comparisons between the two modeling approaches 

include the probability distribution of peak dose, the relative 

importance of parameters, and the value of site-specific data 

in reducing uncertainty.

The test case described in the current report used subsurface 

properties derived from data obtained in the Las Cruces 

Trench experiments (Wierenga et al., 1989; 1990) with a 

source term and exposure scenario loosely based on that 

used in NUREG/CR-6695. The Las Cruces Trench site was 

chosen because the characterization of the unsaturated zone 

at the site was excellent, tracer experimental data and modeling results were available, and the NRC sponsored the 

experiments. 

The simplified modeling approach used the RESRAD code 

(v. 6.1). Although our conclusions are specific to the RESRAD code, we believe similar analyses using other simplified codes (such as DandD and MEPAS) would highlight 

analogous concerns. RESRAD was chosen because it is a 

familiar code to NRC staff and NRC licensees. The more 

complex modeling approach used the STOMP1

 code in a 

stochastic framework of Monte Carlo simulations. STOMP 

was chosen because of its flexibility and extensive feature 

set, its low price (free), and because the authors have experience with it. For comparisons between the two approaches, 

concentrations calculated by STOMP were converted to 

dose consistent with the procedures used in RESRAD.

Monte Carlo simulation was used in this analysis because of 

the complexity of the scenario, which included coupled 

unsaturated-saturated flow and transport and a pumped 

1.0

0.8

0.6

0.4

0.2

0

Empirical cdf

10 20 30 40

Peak Dose (mrem/yr)

Prior

N = 100

Min = 8.9

Max = 49.07

Mean = 15.467

StDev = 7.9565

Median = 11.755

 Prior

 Updated

Updated

N = 100

Min = 8.9

Max = 37.3

Mean = 13.077

StDev = 4.265

Median = 11.76

Figure 1-3. Prior and updated probabilistic results from the RESRAD excavation test case simulations of 

NUREG/CR-6695

1 STOMP web page: http://www.pnl.gov/etd/stomp
Introduction

7

well. Stochastic moment methods (Gelhar, 1993; Dagan and 

Neuman, 1997; Zhang, 2001) provide an alternative to 

Monte Carlo simulation for estimating flow and transport 

uncertainty. These methods have the potential to significantly reduce the computational requirements of the analyses discussed in this report. For our application, however, 

the available stochastic moment methods could not be used 

without unacceptable simplification of the modeled scenario. Recent improvements to address the effects of a coupled unsaturated-saturated system (Sun and Zhang, 2000; 

Zhang and Lu, 2002) suggest that stochastic moment methods may be a viable alternative in the near future.

The following chapter describes the setting of the test case 

and the exposure scenario used to provide relative equivalence in the simplified and complex approaches. Probabilistic simulation results are then presented followed by results 

from a limited sensitivity analysis.

The emphasis in this report series, including the current 

report, is on prediction uncertainty resulting from uncertainty in parameter values. This emphasis is not intended to 

downplay the importance of uncertainty arising from other 

sources such as model structure and future scenarios. The 

importance of conceptual model uncertainty was recognized 

in the first report of this series and played a role in the analysis of the test case results discussed in the second report. 

Conceptual model uncertainty is the subject of an NRC 

research project currently being conducted at the University 

of Arizona. Our test case analyses indirectly include model 

uncertainty since we compare results from a three-dimensional model and a one-dimensional model, a significant 

difference in problem conceptualization. Other uncertainties 

in model structure investigated in this project (through sensitivity analyses) include variations in pumping well location and discharge rates, boundary conditions (affecting 

ambient mean gradient), and variations in geostatistical 

parameters. The NRC staff have also recognized the importance of scenario uncertainty. Although this topic is not 

addressed in detail in this report series, it is touched on in 

this report as we looked at the sensitivity of model results to 

parameters related to the future exposure scenario. 

In addition to discussing primarily parameter uncertainty, 

the analyses in this report series are limited to the hydrologic aspects of decommissioning cases. For buried contaminants, this means that the analysis is limited to the exposure 

pathway involving infiltration of water at the ground surface, leaching of contaminants, and transport of contaminants through the groundwater to a point of exposure. While 

other exposure pathways may be considered in a decommissioning analysis, subsurface water transport is likely to be a 

significant pathway for buried contaminants.

9

2 Test Case Scenario and Hydrogeologic Setting

This chapter describes the simulation domain, the exposure 

scenario, and the hydrogeologic setting of the test case. In 

all cases, decisions about the characteristics of the test case 

were made to result in RESRAD and STOMP models that 

looked as identical as possible, but at the same time allowed 

us to model the test case as realistically as feasible. The 

complex representation, which used the STOMP code, is 

discussed first. The simplified representation, which used 

the RESRAD code, is discussed at the end of this chapter.

2.1 Complex Representation

The test case was based on the hydrogeology at the site of 

the Las Cruces Trench Experiments, New Mexico 

(Wierenga, et al., 1989, 1990). We did not adhere strictly to 

the regional hydrogeology, however, but enforced a water 

table seven meters below the ground surface. In addition, 

we assumed a five-meter-thick aquifer. The region modeled 

in STOMP was a three-dimensional domain of 20 m by 20 

m by 12m deep, as shown in Figure 2-1. In our coordinate 

system, groundwater flow occurred along the x-axis. These 

assumptions about the hydrogeology were done to make the 

test case more comparable to that in NUREG/CR-6695. 

These assumed conditions likely also reduced the differences between the simplified and complex modeling 

approaches. 

To fit the conceptual model embodied in the RESRAD 

code, it was assumed that a source existed as a contaminant 

buried in the top 1 m of soil and that this contaminant would 

leach down with the percolating water when exposed to 

infiltrating precipitation. The near-surface contamination 

covered an area 8 m by 8 m and was located 4 m from the 

upstream boundary along the X axis (centered with respect 

to the Y axis). The contaminated zone was 1-m thick. For 

simplicity, the source consisted of U-234 only with an initial 

concentration of 122 pCi/g in the solid phase. A uniform 

distribution coefficient of Kd = 15 cm3

/g was used. This is 

the same value used in NUREG/CR-6695 for this soil type. 

The source was leached over time according to the sourceterm model used in RESRAD. The resulting solid and liquid 

phase concentrations of the source as a function of time are 

shown in Figure 2-2. Following the RESRAD manual, we 

assigned a half-life time of 2.45 x 105

 years for U-234 in the 

STOMP simulations. 

To conform with RESRAD it was assumed that a well 

tapped the aquifer on the downstream boundary of the contaminated zone. This well was assumed to supply domestic 

water needs only and the sole route of exposure was 

assumed to be through the drinking water pathway. It was 

also assumed that 100% of the drinking water was derived 

from the well and that drinking water ingestion was 2 l/day. 

The pumping rate from the well was 118 m3/yr, the approximate median of the distribution for this quantity given in the 

DandD documentation (Beyeler, et al., 1999). 

2.1.1 Las Cruces Trench Data, Domain, and 

Geostatistics

We started with data and the geostatistics of soil hydraulic 

properties at the Las Cruces site, as presented by Rockhold 

et al. (1996). Based on their findings, we assumed the existence of Miller’s similarity principle, particularly with the 

assumption of identical scaling factors for saturated hydraulic conductivities and for pressure heads. We adopted the 

scaling factor statistics described in Table 1, Figure 6, and 

Equation 18 of Rockhold et al. (1996). We used the Brooks 

and Corey (B-C) water retention model with the Burdine 

hydraulic conductivity model (Eqs. 1 and 2 in Rockhold et 

al., 1996). We used the mean-scaled parameters as determined by Rockhold et al. (1996), that is, Ks

* = 4.16 m/day 

(mean saturated hydraulic conductivity), ψe* = 5 cm (mean 

air-entry value), θs = 0.32 (saturated water content or porosity), θr = 0 (residual water content), b = 0.267 (B-C exponent parameter). The last three parameters were assumed 

uniform over the whole domain (as done in Rockhold et al., 

1996). The scaling theory implies the following relationships:

ψ = ψ*

/α, hence, ψe = ψe

*

/α (2-1)

K = α2

 K*, hence, Ks = α2 Ks

* (2-2)

where α is the scaling factor. α does not depend on saturation, but it is heterogeneous or random in a similar manner 

to K or Ks. Hence, rather than generating random fields of 

Ks and ψe, we generated a random field of α based on the 

geostatistics presented in Rockhold et al. (1996), and then 

Figure 2-1. The domain for the test case simulations

Vadose zone

20 m

20 m

Aquifer

Well

8 m

8 m 

U-234

GW 

Flow

7 m

5 m
Test Case Scenario and Hydrogeologic Setting

10

translated it to fields of Ks and ψe (the necessary input 

parameters) using the above relationships (Eqs. 2-1 and 

2-2). 

The particle density was also assumed uniform at 2.65 

g/cm3, which implied a bulk density of 1.8 g/cm3

. 

2.1.2 Choice of Random Field Generator, 

Resulting Geostatistics of α and Ks and Spatial 

Distribution of Permeabilities

There are several methods to generate correlated random 

fields. The most commonly used are methods based on 

Turning Bands, LU decomposition, fast Fourier transform 

(FFT), and sequential simulations (also called the method of 

conditional distribution). We chose SGSIM, a sequential 

Gaussian random field generator from GSLIB (Deutsch and 

Journel, 1997) because it allows nested, anisotropic, conditional simulations. At the core of each random field generator there is a random number generator, which generates 

independent pseudo random numbers based on three deterministic parameters, including a seed – the only random 

number generator parameter that can be changed by the 

user. A good discussion on random number generators is 

presented in Press et al. (1986). The random field generator 

uses the uncorrelated numbers generated by the random 

number generator and transforms them into correlated values in space (or time), based on a correlation structure (in 

the form of a covariance function or a variogram) provided 

by the user. Descriptions of the different random field generators can be found in Deutsch and Journel (1997), and an 

extensive comparison between different random field generators and random number generators can be found in Orr 

(1993). 

We used SGSIM to generate random realizations of log-α

using the normalized variogram model defined by Eq. 18 of 

Rockhold et al. (1996) (after correction), that is: 

(2-3)

where

and h0 and h90 are lag distances in the horizontal and vertical directions, respectively. Simulations were carried out on 

a 1-m grid over the domain of Figure 2-1. (See Section 2.1.4 

for a discussion of the numerical grid.)

The resulting variograms are presented in Figure 6 of Rockhold et al. (1996) and Figure 2-3 below (after translating the 

integral scales into effective ranges by multiplying them by 

3). Eq. 2-3 implies a nested variogram with statistical 

anisotropy expressed by principal vertical and horizontal 

directions, with statistical isotropy in the horizontal plane. 

In Figure 2-3, note the relatively good fit between experimental and theoretical variograms for a single realization; 

however, the experimental variogram is presumed isotropic, 

horizontally, while in reality, generated single realizations 

whose size is less than 40 correlation scales almost always 

exhibit anisotropy (apparently, due to imperfection of all 

random field generators). Note also that the vertical variogram shows practically no correlation between block values 

in the vertical direction. Sample statistics of a single realization of α are given in Table 2-1.

Next, we used Eqs. 2-1 and 2-2 (above) to generate realizations of Ks and ψe, a different value for each block (or element). Figure 2-4 shows the geostatistics of log(Ks). As 

expected, the normalized variograms almost resemble those 

of log(α). Sample statistics of log(Ks) and Ks for a single 

realization are given in Table 2-2.

0

20

40

60

80

100

120

140

0 100 200 300 400 500

Time (yr)

Solid Phase Source Conc. (pCi/g)

0.0

1000.0

2000.0

3000.0

4000.0

5000.0

6000.0

7000.0

8000.0

9000.0

0 100 200 300 400 500

Time (yr)

Liquid Phase Source Conc (pCi/l)

Figure 2-2. Contaminant source (U-234) concentrations in the solid phase (left) and the liquid phase (right)

γ( ) h 0.132 0.543 1 e h1 – ( ) – 0.325 1 e h2 – = + + ( ) –

h1 h0 ( ) ⁄ 6.347 2 h90 ( ) ⁄ 0.039 2 [ ] +

1 2⁄ =

h2 h0 ( ) ⁄ 1.089 2 h90 ( ) ⁄ 1.129 2 [ ] +

1 2⁄ =
Test Case Scenario and Hydrogeologic Setting

11

Typically, in Monte Carlo simulations every element is considered a zone. Consequently, each element (or basic block) 

possesses its own hydraulic and transport properties. 

Figure 2-5 shows one realization of the permeability field. 

These values, as well as values of ψe and the other (uniform) parameters, were written into STOMP’s input file in a 

generic way that enables us to regenerate the input automatically (4,800 values for each parameter) for each new realization. Our sample statistics started with 300 simulations 

for each case we analyzed. After examining the results, particularly the (sample) ensemble variances vs. number of 

simulations, we concluded that the number of simulations 

could be reduced to 200 with relatively insignificant reduction in accuracy (of ensemble mean and variance) and significant savings in computer time. We checked the output at 

different locations and particularly at the location of the 

pumping well.

2.1.3 Flow and Transport Parameters

Table 2-3 presents the hydraulic and transport parameters 

used in our models, following the parameters presented by 

Rockhold et al. (1996) for the Las Cruces Trench site. For 

each realization, different Ks and ψe values were assigned to 

each element, while other flow parameters, such as porosity, 

were assigned uniformly. This is due to the dominating variability of the saturated hydraulic conductivity as indicated 

by the statistics presented by Rockhold et al. (1996). The 

coefficient of variation of porosity, for example, is smaller 

than the coefficient of variation of Ks by an order of magnitude, implying negligible variations in porosity. Other transport parameters included longitudinal and transverse 

dispersivities of 0.5 m and 0.1 m, respectively.

A recharge rate of 25.4 cm/yr (or 10 inch/yr) was imposed 

as a Neumann boundary condition at the top of the domain. 

0

50

100

150

200

250

-3.87

-3.35

-2.83

-2.31

-1.79

-1.27

-0.75

-0.23

0.29

0.81

1.33

1.85

2.37

2.89

Normalized values 

Count

Figure 2-3. Variograms (top) and histogram (bottom) of log(α) for a single realization with the analytical 

directional variograms shown for comparison (H=horizontal, V=vertical)

0

0.2

0.4

0.6

0.8

1

1.2

0 5 10 15 20

Lag (m)

γ

X-Vario

Y-Vario

Z-Vario

Analytic H

Analytic V
Test Case Scenario and Hydrogeologic Setting

12

Initial conditions assumed gravitational/static conditions in 

the vadose zone and steady state flow in the aquifer. Boundary conditions included prescribed heads on aquifer boundaries, creating a gradient of 0.02. All other boundaries were 

considered no-flow. In terms of transport, we assumed 

boundary conditions of zero concentration gradient on all 

the vertical side surfaces, implying no dispersive transport 

(only advective transport). A time-varying concentration 

Dirichlet condition was imposed over the 8 m x 8 m contaminated zone. Liquid phase concentration as a function of 

time was presented above.

We ran flow simulations with STOMP until steady state was 

reached (in less than 1000 years), and then used the steady 

velocity fields as the basis for transport simulations, using 

STOMP. This was done in order to save computation time, 

keeping in mind (a) this is not too far from reality (natural 

recharge had been there long before contamination started), 

and (b) the goal of finding long-term effects. This is not a 

recommended shortcut in cases where short-term predictions are desired and variations in recharge events dominate 

the flow.

2.1.4 Numerical Modeling Considerations with 

STOMP: Grid Resolution, Time Step, and 

Solver

The basic block (element) size for the STOMP modeling 

was a 1-m cube. With a domain size of 20 m x 20 m x 12 m, 

this translates to 4800 blocks with different hydraulic properties (Ks and ψe). These generated properties were fed into 

STOMP for flow and transport simulations. The choices of 

domain size and block size were a compromise between a 

desire for high resolution of the randomly varying parameters (including adequate representation of the geostatistics) 

and the need to be able to conduct Monte Carlo simulations 

in a reasonable amount of time on the available computers. 

The chosen mesh size imposes certain numerical considerations as discussed below.

0

50

100

150

200

250

-3.65

-3.02

-2.39

-1.75

-1.12

-0.48

0.15

0.79

1.42

2.05

2.69

3.32

3.96

4.59

Log(Ks)

Count

Figure 2-4. Variograms (top) and histogram (bottom) of log(Ks) (m/day) (single realization)

0

0.2

0.4

0.6

0.8

1

1.2

0 5 10 15 20

Lag (m)

γ

K - X Dir

K - Y Dir

K - Z Dir
Test Case Scenario and Hydrogeologic Setting

13

Peclet Number (Pe): According to the scheme and solver 

used in STOMP, the maximum Peclet number allowed is 

Pe = 2. Since Pe = v dx / D = dx / αL = ratio of grid spacing 

to longitudinal dispersivity, the minimum allowed dispersivity is half the cell size. Following Bear (1979, p. 261-

263), the numerical dispersivity equals half the mesh size, 

which also corresponds to Pe = 2. Given the chosen block 

Table 2-1. Statistics of log(α) for a single 

realization

Mean 0.0341

Median 0.0344

Standard Deviation 1.01

Kurtosis -0.0595

Skewness 0.000947

Range 7.17

Minimum -3.87

Maximum 3.31

Count 4800

0

5

10

Z-Direct oi n Node Positions, m

0

5

10

15

20

X or R-Direction Node Positions, m

0 5 10 15 20 Y or Theta-Direction Node Positions, m

X

Y

Z

Intrinsic Permeability, m2

1.81E-10

1.76E-10

1.71E-10

1.67E-10

1.62E-10

1.57E-10

1.53E-10

1.48E-10

1.43E-10

1.39E-10

1.34E-10

1.29E-10

1.25E-10

1.20E-10

1.15E-10

1.11E-10

1.06E-10

1.01E-10

9.66E-11

9.20E-11

8.73E-11

8.26E-11

7.80E-11

7.33E-11

6.86E-11

6.39E-11

5.93E-11

5.46E-11

4.99E-11

4.53E-11

1.00E-11

5.00E-12

1.00E-12

5.00E-13

1.00E-13

5.00E-14

Frame 001  19 Sep 2001  permeability distribution Las Cruces, # 0001

Figure 2-5. Permeability distribution of a single realization (Las Cruces geostatistics)

Table 2-2. Statistics of log(Ks) and Ks (m/day) for a 

single realization

log(Ks) Ks

Mean 1.10 6.38

Median 1.10 3.02

Standard Deviation 1.23 10.7

Kurtosis -0.06 47.2

Skewness 0.00 5.47

Range 8.75 163.6

Minimum -3.65 0.03

Maximum 5.10 163.59

Count 4800 4800
Test Case Scenario and Hydrogeologic Setting

14

size, this Peclet number constraint determined the longitudinal dispersivity of αL = 0.5.

Courant Number (Co): The Courant number condition 

imposes a limit on the maximum size of the time step in a 

numerical simulation. Co = v (dt/dx) = (q/θ)(dt/dx) < 4 was 

suggested by El-Kadi and Ling (1993) and El-Kadi (1995). 

For the STOMP code, Co = q (dt/dx) < 1 has been suggested.1

 The largest flux (q) we encountered in our simulations is about 4 cm/hr in the saturated zone (θ = 0.32). 

Following El-Kadi’s formulation, since dz = dx = 1 m = 

100 cm, dt < (0.32/4)x100x4 = 32 hours. However, such 

high flux rates are limited to the vicinity of the pumping 

well, while fluxes in the vadose zone are dominated by the 

recharge (the boundary condition), which is q ≅ 

0.0029 cm/hr. Statistics of the vertical flux are given in 

Table 2-4 for a uniform domain and two different random 

realizations, one with a low variance (of α) and one with a 

high variance. A histogram of the vertical flux for the lowvariance random realization is shown in Figure 2-6. The 

mean and median values given in Table 2-4 are smaller than 

the local maximum rate around the well by three orders of 

magnitude, implying much higher allowable time steps. 

Table 2-3. Model parameters used in the STOMP 

simulations

Parameter 

Value (deterministic or 

stochastic)

Block size 1 m (X, Y, and Z)

Recharge 254 mm/yr (10 inch/yr)

GW gradient 0.02

Pumping rate 118 m3

/yr

Log α - mean 0

Log α - Std. Dev. 0.61

θs 0.32

θr 0

β (Brooks-Corey) 0.267

Initial source concentration 122 pCi/g (solid) = 8133 

pCi/l (liquid) 

U-234 Half-life 2.45 x 105 years

Kd 15 cm3

/g

Mean-scaled Ks* (geometric mean)

4.16 m/day = 4.8E-3 cm/s 

Mean air entry ψe* 5 cm

Particle density 2.65 g/cm3

Longitudinal dispersivity 0.5 m

Transverse dispersivity 0.1 m

1 M.D. White, personal communication, May 2001.

Table 2-4. Statistics of vertical fluxes (cm/hr) in uniform 

and random realizations (Las Cruces 

statistics)

 Uniform

Low 

Variance

High 

Variance

Mean 0.00208 -0.00662 -0.00498

Median 0.00288 -0.00297 -0.00277

Std. Dev. 0.00107 0.0667 0.0577

Kurtosis -0.872 32.4 63.8

Skewness -0.878 -1.40 1.26

Range 0.00299 1.74 1.78

Minimum 0 -0.933 -0.805

Maximum 0.00299 0.812 0.973

Count 4800 4400 4800

1

10

100

1000

10000

-0.933-0.801-0.669-0.536-0.404-0.272-0.140-0.007

0.1250.2570.3890.5220.6540.786

Vz (cm/hr)

Count

Figure 2-6. Histogram of vertical flux in a low-variance, random realization
Test Case Scenario and Hydrogeologic Setting

15

We experimented with transport simulations, incrementing 

time steps from dt = 30 up to dt = 5000 hours. A few results 

using an extreme random realization are shown in 

Figure 2-7. These are concentrations over time at the top 

node of the well screen, at the node just above the water 

table at the well location, and at a node downstream from 

the well, 2.5-m from the aquifer bottom, half way between 

the contaminated zone and the outflow boundary. The 

results indicate practically identical transport at least up to 

dt = 2600 hours. We found similar results from simulations 

of a more variable random field (with twice the standard 

deviation). Nevertheless, we selected the value of 600 hours 

as the maximum time step allowed in most transport simulations.

As was stated earlier, the choice of grid resolution of 1 m3

was a compromise. High resolution is desirable for two reasons: (a) to avoid upscaling of parameters obtained from 

core samples, and (b) to include as many elements as possible within an integral scale, such that the spatial variability 

is captured and is represented correctly by the spatial variogram. Several studies (e.g., Orr, 1993) have shown the need 

for at least 5 elements, preferably 10, per integral scale. A 

block size of 1 m implies upscaling of hydraulic properties 

from the smaller core size. The 1-m block size still provides 

a fair representation of lateral correlation, with an integral 

scale of about 6.4 m (from Eq. 2-3), but misses the vertical 

correlation, with almost one block per integral scale. This is 

not too disturbing because the short correlation scale 

implies almost no correlation vertically anyway. Note that 

within blocks there is perfect vertical correlation. Another 

important factor is the effect of boundary conditions governed by the distance of stresses and responses in the aquifer 

from the boundaries. A distance of at least three to five integral scales is needed in order to ignore boundary effects 

(Orr, 1993). Our random fields are generally below this 

threshold with about three correlation scales across the 

domain. 

Due to the limited size of the domain, we do not expect it to 

be ergodic; this means that we cannot imply ensemble statistics from spatial statistics and vice versa. Nevertheless, 

the univariate statistics and spatial (experimental) variograms were shown (in Section 2.1.1) to be close to the theoretical (ensemble) values. The Monte Carlo simulation 

results presented in Chapters 3 and 4 are all sample-ensemble quantities (mean, variance, covariance).

In order to have an idea about possible deviations from a 

more rigorous resolution, we ran a few simulations with 

double the grid resolution, i.e., 40 x 40 x 24 blocks (block 

size of 0.5-m on each side). Figure 2-8 shows results as concentration profiles with depth at the (x, y) location of the 

well for the two resolutions using a uniform domain. The 

bottom 5 m in each plot represents the well (and the saturated zone). The similarity between the results leads us to 

the conclusion that the 1-m block size is acceptable. The 

difference in the location of the peak is due to the effective 

difference in the location of the bottom of the contaminated 

zone: 1.0 m for the low resolution grid and 0.5 m for the 

high resolution grid. Note that this comparison does not 

address the question of domain size relative to the correlation scale and the overall effects of boundary conditions.

The matrix solver of preference in STOMP for large and 

complex domains is SPLIB – a preconditioned conjugate 

gradient solver. STOMP also includes different weighting 

60

50

40

30

20

10

0

Concentration (pCi/L)

0 200 400 600 800 1000

Time (yr)

Node above well Top node of well

Downstream in aquifer

Figure 2-7. Comparison of concentrations at three locations in the domain using three different maximum time step 

sizes (300, 900, and 2600 years) for a random realization (#81). Note, concentration differences are 

negligible.
Test Case Scenario and Hydrogeologic Setting

16

schemes for transport calculations, namely Patankar (the 

default), TVD, and Superbee. We experimented and compared the performances of these schemes in an attempt to 

maximize time step size, minimize dispersion, and optimize 

grid resolution. Results from our experiments are shown in 

Figure 2-9. The figure shows concentrations of U-234 over 

time (breakthrough curves) at three different locations – the 

top node of the well screen, one node just above the well, 

and one point downstream, half way between the well and 

the outflow boundary. A close look at the breakthrough 

curves, and a comparison with breakthrough curves generated by simulations with different maximum time steps 

revealed high sensitivity of the TVD and Superbee schemes 

to the maximum time step, while the Patankar scheme 

showed robustness and consistency over a large range of 

maximum time steps. Therefore, we used the Patankar 

scheme in all our STOMP simulations.

2.1.5 Sample Results from Single Realizations 

This section shows a sample of some of the results obtained 

from single realization simulations with STOMP. 

Figures 2-10 and 2-11 show distributions of saturation, 

pressure head, and total head resulting from steady state 

flow (imposed by steady recharge) in a nonuniform realization (of hydraulic conductivities and air-entry values) from 

the random field described above. The heterogeneous nature 

of the domain is clearly illustrated along with the location of 

the water table. Figure 2-12 shows the evolution of the U234 plume within the nonuniform domain over 900 years in 

100-year increments. For this particular plume realization, it 

can be seen that there is a significant concentration gradient 

with depth in the aquifer and that a portion of the contaminant is not captured by the pumping well. In addition, concentrations are highest in the unsaturated zone. 

Breakthrough curves resulting from flow and transport simulations in a single random realization are shown in 

Figure 2-13. Concentrations at four locations in the domain 

are shown: in the unsaturated zone directly above the well, 

at the top and bottom nodes of the well, and downstream at 

the top of the aquifer, half-way between the source and the 

boundary. The average concentration in the five nodes comprising the well is also shown. The behavior observed in 

Figure 2-12 is seen here in detail. The concentration in the 

unsaturated zone is significantly larger than in the aquifer. 

In addition, the concentration at the top of the well is much 

larger than at the bottom of the well. The peak average well 

concentration for this realization was about 30 pCi/l. 

To provide additional detail on the complex model results, 

Figure 2-14 shows a series of concentration profiles with 

depth at the location of the well as obtained from a STOMP 

simulation for a single random realization. These profiles 

illustrate the limited vertical mixing in the aquifer and the 

significant dilution in the aquifer due to horizontal flow. 

2.2 Simplified Representation

The simplified representation of the test case site was made 

as close to the complex representation discussed above as 

12

10

8

6

4

2

0

Vertical Distance (m)

0 200 400 600 800 1000 1200

Concentration (pCi/L)

 0.5-m Grid

 1.0-m Grid

Figure 2-8. Concentration depth profiles at the (x, y) location of the well using a 1-m grid and a 0.5-m grid for a 

uniform domain after 500 years
Test Case Scenario and Hydrogeologic Setting

17

possible. The primary differences were required by the conceptual model embodied in RESRAD and are summarized 

here.

• One-dimensional flow and transport. RESRAD is limited to a one-dimensional analysis.

• Uniform domain. RESRAD allows for multiple layers in 

the domain, but since there were no distinct layers 

observed at the Las Cruces Trench site, the test case was 

modeled as statistically homogeneous in the complex 

representation and as homogeneous in the simplified 

representation.

• Unit hydraulic gradient in the unsaturated zone. RESRAD makes the unit gradient assumption to simplify the 

unsaturated flow solution. This assumption is not 

required by STOMP, which solves the Richards equation 

for flow in the unsaturated zone.

• Simplified aquifer mixing model. RESRAD uses a simplified analysis to estimate mixing and dilution in the 

aquifer. The primary reason a three-dimensional analysis 

was adopted for the complex representation was to more 

accurately represent the effect of the pumping well.

The RESRAD domain consisted of a 1-m thick contaminated zone above a 6-m unsaturated zone (7 m total of 

unsaturated soil) and a 5-m saturated zone. There was no 

cover over the contaminated zone, which was 8 m x 8 m in 

area. Based on the characterization of the Las Cruces 

Trench site soils, all soils in the test case (contaminated 

zone, unsaturated zone, and aquifer) were classified as a 

sand. Precipitation was set at 25.4 cm/y with irrigation, 

evapotranspiration, and runoff set to zero, resulting in the 

required recharge rate. As discussed above, the source term 

was calculated using RESRAD.

In the RESRAD Monte Carlo simulations, the hydraulic 

properties of the unsaturated zone and the aquifer were 

modeled as random parameters. This is the corollary of the 

three-dimensional Monte Carlo simulations carried out with 

STOMP. Note that statistically, however, the Monte Carlo 

simulations have different meanings for the two representations. The RESRAD model (1-D) assumes univariate statistics, i.e., either multiple scenarios (an ensemble) of uniform 

domains, or multiple parallel non-interacting columns averaged in space. The STOMP model (3-D) assumes a multivariate joint distribution, i.e., multiple scenarios of spatially 

correlated random realizations (a spatially correlated random field).

Two sets of distributions of the soil hydraulic parameters 

were used for the RESRAD simulations. One of these was a 

set of generic distributions taken from NUREG/CR-6565 

(Meyer et al., 1997). Distributions for the sand classification 

were selected. For the saturated hydraulic conductivity, lognormal distributions with the mean and standard deviation 

given in NUREG/CR-6565 were used instead of the specified beta distributions. These parameter distributions are 

given in Tables 2-5. Correlations between parameters 

greater than 0.6 as given in NUREG/CR-6565 were 

included in the simulations.

The second set of distributions used with RESRAD was 

based on data from the Las Cruces Trench site. These 

parameter distributions are given in Table 2-6. Compared to 

the generic distributions, the following changes were made. 

Mean bulk density was set to the value implied by the 

300

250

200

150

100

50

0

Concentration (pCi/L)

0 200 400 600 800 1000

Time (yr)

 Patankar

 TVD

 Superbee

Node above well

Top node of well

Downstream in aquifer

Figure 2-9. Comparison between different weighting schemes for transport simulations at three locations in the 

domain
Test Case Scenario and Hydrogeologic Setting

18

0

5

10

0

Z-Direction Node Positions, m

5

10

15

20 X or R-Direction Node Positions, m

0

5

10

15

20

Y or Theta-Direction Node Positions, m

X Y

Z

Aqueous Pressure(Pa), pa

143849.88

142328.77

140807.65

139286.53

137765.42

136244.30

134723.18

133202.06

131680.95

130159.83

128638.71

127117.60

125596.48

124075.36

122554.25

121033.13

119512.01

117990.89

116469.78

114948.66

113427.54

111906.43

110385.31

108864.19

107343.08

105821.96

104300.84

102779.72

101258.61

99737.49

98216.37

96695.26

95174.14

93653.02

92131.91

90610.79

89089.67

87568.55

86047.44

84526.32

Frame 001 26 May 2001 Saturations, pressure and hydraulic heads for realization #200

0

5

10

0

Z-Direction Node Positions, m

5

10

15

20 X or R-Direction Node Positions, m

0

5

10

15

20

Y or Theta-Direction Node Positions, m

X Y

Z

Aqueous-Phase Saturation

0.98

0.97

0.95

0.94

0.92

0.90

0.89

0.87

0.85

0.84

0.82

0.81

0.79

0.77

0.76

0.74

0.72

0.71

0.69

0.68

0.66

0.64

0.63

0.61

0.60

0.58

0.56

0.55

0.53

0.51

0.50

0.48

0.47

0.45

0.43

0.42

0.40

0.38

0.37

0.35

Frame 001 26 May 2001 Saturations, pressure and hydraulic heads for realization #200

Figure 2-10. Saturation (top) and pressure (bottom) distributions in a nonuniform (random) realization 

(Var[ln(Ks)] = 1.5)
Test Case Scenario and Hydrogeologic Setting

19

0

5

10

0

Z-Direction Node Positions, m

5

10

15

20 X or R-Direction Node Positions, m

0

5

10

15

20

Y or Theta-Direction Node Positions, m

X Y

Z queous Hydraulic Head (Fresh Water),

11.21

11.04

10.88

10.71

10.55

10.38

10.22

10.05

9.89

9.72

9.56

9.39

9.22

9.06

8.89

8.73

8.56

8.40

8.23

8.07

7.90

7.74

7.57

7.41

7.24

7.08

6.91

6.75

6.58

6.42

6.25

6.09

5.92

5.76

5.59

5.43

5.26

5.10

4.93

4.76

Frame 001  26 May 2001  Saturations, pressure and hydraulic heads for realization #200

Figure 2-11. Total head distribution for the same realization as Figure 2-10
Test Case Scenario and Hydrogeologic Setting

20

0

5

10

0

Z-Direction Node Positions, m

5

10

15

20 X or R-Direction Node Positions, m

0

5

10

15

20

Y or Theta-Direction Node Positions, m

X Y

Aqueous-Phase u234 Concentration, 1/liter Z

3788

3452

3145

2866

2611

2379

2168

1975

1800

1640

1494

1361

1240

1130

1030

938

855

779

710

647

589

537

489

446

406

370

337

307

280

255

232

212

193

176

160

146

133

121

110

101

92

84

76

69

63

58

52

48

44

40

Figure 2-12. Snapshots of the simulated U-234 plume including transport from the source, extraction from the well, 

and radioactive decay. The snapshots are in 100-year increments, from 100 years to 900 years after the 

initial release.

0

5

10

0

Z-Direction Node Positions, m

5

10

15

20 X or R-Direction Node Positions, m

0

5

10

15

20

Y or Theta-Direction Node Positions, m

X Y

entration, 1/liter Z

0

5

1

0

5

10

15

20 X or R-Direction Node Positions, m

0

5

10

15

20

Y or Theta-Direction Node Positions, m

X

Aqueous-Phase u234 Concentration, 1/liter

3968

3616

3295

3002

2735

2492

2271

2069

1885

1718

1565

1426

1300

1184

1079

983

896

816

744

678

617

563

513

467

426

388

353

322

293

267

244

222

202

184

168

153

139

127

116

105

96

88

80

73

66

60

55

50

46

42

20

5

0

5

10

0

Z-Direction Node Positions, m

5

10

15

20 X or R-Direction Node Positions, m

0

5

10

15

20

Y or Theta-Direction Node Positions, m

X Y

Aqueous-Phase u234 Concentration, 1/liter Z 3968

3616

3295

3002

2735

2492

2271

2069

1885

1718

1565

1426

1300

1184

1079

983

896

816

744

678

617

563

513

467

426

388

353

322

293

267

244

222

202

184

168

153

139

127

116

105

96

88

80

73

66

60

55

50

46

42

20

5

Frame 001  28 Apr 2001  500 years

0

5

1

0

5

10

15

20 X or R-Direction Node Positions, m

0

5

10

15

20

Y or Theta-Direction Node Positions, m

X

Aqueous-Phase u234 Concentration, 1/liter

3968

3616

3295

3002

2735

2492

2271

2069

1885

1718

1565

1426

1300

1184

1079

983

896

816

744

678

617

563

513

467

426

388

353

322

293

267

244

222

202

184

168

153

139

127

116

105

96

88

80

73

66

60

55

50

46

42

20

5

Frame 001  28 Apr 2001  600 years

0

5

10

0

Z-Direction Node Positions, m

5

10

15

20 X or R-Direction Node Positions, m

0

5

10

15

20

Y or Theta-Direction Node Positions, m

X Y

entration, 1/liter Z

0

5

10

0

Z-Direction Node Positions, m

5

10

15

20 X or R-Direction Node Positions, m

0

5

10

15

20

Y or Theta-Direction Node Positions, m

X Y

Aqueous-Phase u234 Concentration, 1/liter Z 3968

3616

3295

3002

2735

2492

2271

2069

1885

1718

1565

1426

1300

1184

1079

983

896

816

744

678

617

563

513

467

426

388

353

322

293

267

244

222

202

184

168

153

139

127

116

105

96

88

80

73

66

60

55

50

46

42

20

5

Frame 001  28 Apr 2001  800 years

0

5

1

0

5

10

15

20 X or R-Direction Node Positions, m

0

5

10

15

20

Y or Theta-Direction Node Positions, m

X

Aqueous-Phase u234 Concentration, 1/liter

3968

3616

3295

3002

2735

2492

2271

2069

1885

1718

1565

1426

1300

1184

1079

983

896

816

744

678

617

563

513

467

426

388

353

322

293

267

244

222

202

184

168

153

139

127

116

105

96

88

80

73

66

60

55

50

46

42

20

5

Frame 001  28 Apr 2001  900 years

0

5

10

0

Z-Direction Node Positions, m

5

10

15

20 X or R-Direction Node Positions, m

0

5

10

15

20

Y or Theta-Direction Node Positions, m

X Y

ntration, 1/liter

ZZ

X Y

U-234 Concentration Increasing

100 yr. 200 yr. 300 yr.

400 yr. 500 yr. 600 yr.

700 yr. 800 yr. 900 yr.
Test Case Scenario and Hydrogeologic Setting

21

400

300

200

100

0

Concentration (pCi/l)

0 200 400 600 800 1000

Time (yr)

 Node Above Well

 Top Node of Well

 Well Average

 Downstream Node

 Bottom Node of Well

Figure 2-13. Breakthrough curves from a single random realization at four locations in the domain and for 

the average of the five well nodes

0 500 1000 1500 2000 2500 3000 3500

Concentration (pCi/L)

12

10

8

6

4

2

0

Vertical Distance (m)

 100 yrs

 200 yrs

 300 yrs

 400 yrs

 500 yrs

 600 yrs

 700 yrs

 800 yrs

 900 yrs

 1000 yrs

Figure 2-14. Vertical profiles of concentration at the (x, y) location of the well as calculated by STOMP for a 

randomly heterogeneous realization
Test Case Scenario and Hydrogeologic Setting

22

measured average porosity of 0.32 and a particle density of 

1.8 g/cm3. The Las Cruces Trench sample mean and 

variance of the saturated water content were used for the 

porosity moments. The effective porosity distribution was 

set equal to the total porosity distribution. The distribution 

of the unsaturated zone soil-type exponent was based on the 

van Genuchten water retention parameters derived from the 

Las Cruces Trench data. The soil-type exponent (b) was 

calculated from the van Genuchten parameters using the 

relationship presented in NUREG/CR-6565 (Appendix C). 

Field-measured values of saturated hydraulic conductivity 

were used as the basis for the Ks distributions. Lognormal 

distributions were assumed for both b and Ks. The sample 

mean values of ln(b) and ln(Ks) were taken as the geometric 

mean values of the parameter distributions. The variance of 

the ln(b) and ln(Ks) parameter distributions were 

determined from a relationship presented in Gilbert (1987; 

Eq. 4.14) for the variance of the mean. 

(2-4)

where

=the variance of the estimated mean site value 

for ln(b) or ln(Ks)

σ2 = 0.27682 and 1.21342, the sample variances of 

the ln(b) and ln(Ks) measurements, respectively

n = 448, the number of measurements

ρc = 0.033, the average cross-correlation between 

the measurement locations

This relationship considers the multitude of samples at the 

Las Cruces Trench site and the spatial correlation of those 

samples. The value of ρc was derived from the spatial covariance of Ks (discussed above) and the known measurement 

locations from the Las Cruces Trench. The same value of ρc 

was used for ln(b). The means and standard deviations listed 

in Table 2-6 for the soil-type exponent and the saturated 

hydraulic conductivity parameter distributions were calculated from the estimated mean and variance of the mean for 

ln(b) and ln(Ks) (see, for example, NUREG/CR-6565, 

Appendix A). The same correlations between parameters 

used with the generic distributions were applied to the sitespecific distributions. 

RESRAD has two options for modeling mixing and dilution 

in the aquifer, referred to as the Mass Balance and the Nondispersion models. Based on the size of the contaminated 

zone (8 m x 8 m) and the guidance presented in the RESRAD documentation, the Mass Balance model would be 

recommended for this test case. Deterministic results using 

the site-specific mean parameter values are shown for both 

models in Figure 2-15. (Note that the scaled mean values for 

Ks and b=1/β from Table 2-3 were used in the deterministic 

RESRAD simulations, not the mean values from Table 2-6.) 

This figure shows well concentrations as a function of time 

and illustrates a number of points. Well concentrations from 

RESRAD reflect the shape of the source concentration 

(Figure 2-2) for a contaminant that experiences little decay 

over the time period of the simulation. In addition, the Mass 

Balance model results in a significantly larger peak concentration in the well, ten times larger in this case. By comparing the RESRAD results to the STOMP results in 

Figure 2-13, it can be seen that both RESRAD models predict a higher peak well concentration. The peak of the well 

concentration for this single STOMP realization was 

approximately 30 pCi/l, about one-third of the RESRAD 

Nondispersion model result. (As will be shown in 

Chapter 3, the ensemble average peak well concentration 

Var x( ) σ

2

n

----- 1 ρc = { } + ( ) n – 1

Var x( )

Table 2-5. Generic stochastic parameter inputs used in RESRAD simulations - default distributions from NUREG/CR6565 for a sand soil

Parameter Distribution Mean Std. Dev. Min. Max.

Unsat. Zone Bulk Density (g/cm^3) normal 1.58 0.158 1.092 2.068

Unsat. Zone Total Porosity normal 0.43 0.06 0.245 0.615

Unsat. Zone Effective Porosity normal 0.383 0.061 0.195 0.571

Unsat. Zone Field Capacity normal 0.061 0.015 0.015 0.107

Unsat. Zone Sat. Hyd. Conductivity (m/yr) lognormal 2604.5 1388.0 490.3 10775.2

Unsat. Zone Soil-Type Exponent lognormal 0.998 0.218 0.500 1.901

Saturated Zone Bulk Density (g/cm^3) normal 1.58 0.158 1.092 2.068

Saturated Zone Total Porosity normal 0.43 0.06 0.245 0.615

Saturated Zone Effective Porosity normal 0.383 0.061 0.195 0.571

Saturated Zone Sat. Hyd. Conductivity (m/yr) lognormal 2604.5 1388.0 490.3 10775.2
Test Case Scenario and Hydrogeologic Setting

23

was approximately 25 pCi/l.) It can also be observed, however, that with RESRAD there is an inherent delay in arrival 

of contaminant at the well until the peak arrives. As a result 

of dispersion, STOMP predicts a much earlier first arrival of 

contaminant at the well.

The RESRAD Mass Balance model assumes instantaneous 

distribution of contaminant throughout the depth of the 

aquifer and that all contaminant is withdrawn through the 

well. Figures 2-12 to 2-14 illustrate that these conditions are 

not reflected in the complex representation results. Contaminant can be seen being transported out of the domain in this 

figure. Only limited vertical mixing in the aquifer is indicated by the results of the complex representation. Based on 

the observations from the single realization STOMP simulations, it was decided that the RESRAD Nondispersion 

model was more representative of the STOMP results. Only 

results from the Nondispersion model are reported in the 

remainder of this report. 

Table 2-6. Site-specific stochastic parameter distributions used in the RESRAD simulations - based on the Las Cruces 

Trench dataset where applicable (compare to Table 2-5)

Parameter Distribution Mean Std. Dev. Min. Max.

Unsat. Zone Bulk Density (g/cm^3) normal 1.8 0.158 1.312 2.288

Unsat. Zone Total Porosity normal 0.32 0.034 0.215 0.425

Unsat. Zone Effective Porosity normal 0.32 0.034 0.215 0.425

Unsat. Zone Field Capacity normal 0.061 0.015 0.015 0.107

Unsat. Zone Sat. Hyd. Conductivity (m/yr) lognormal 1559.3 360.2 751.1 3073.3

Unsat. Zone Soil-Type Exponent lognormal 3.239 0.169 2.755 3.799

Saturated Zone Bulk Density (g/cm^3) normal 1.8 0.158 1.312 2.288

Saturated Zone Total Porosity normal 0.32 0.034 0.215 0.425

Saturated Zone Effective Porosity normal 0.32 0.034 0.215 0.425

Saturated Zone Sat. Hyd. Conductivity (m/yr) lognormal 1559.3 360.2 751.1 3073.3

0

20

40

60

80

100

120

600 700 800 900 1000

Time (yr)

Well Conc. (pCi/L)

0

200

400

600

800

1000

1200

600 700 800 900 1000

Time (yr)

Well Conc. (pCi/L) Figure 2-15. Concentration in the well for the RESRAD Mass Balance (left) and Nondispersion (right) models

25

3 Base Case Monte Carlo Simulation Results

This chapter presents the results of Monte Carlo simulations 

for the complex representation of the test case (using 

STOMP) and the simplified representation (using RESRAD). RESRAD v. 6.0 has a built-in Monte Carlo simulation capability. STOMP, however, required some additional 

work to allow it to be used in a Monte Carlo mode with 

minimal user interaction. Several programs and scripts in 

Fortran and C-Shell were written in order to perform a variety of tasks: translating generated random fields of scaling 

factors to other parameters (Ks and ψe); writing these to 

valid STOMP input files; extracting the required results 

from the STOMP output files; and calculating the ensemble 

statistics. STOMP simulations were run on a number of 

machines. Sequential runs were carried out on SGI and 

SUN Unix machines. In addition, a “task farm” script and a 

C++ program were written for parallel processing on one of 

the supercomputers at PNNL. Because of the independence 

of each realization, one of the available programs that utilize free processor time on a network of PCs or Unix 

machines could also be used. 

3.1 STOMP Ensemble Mean Results

We started with the base case of mild variability, with variance of log(Ks) equal to 1.5. We generated 300 realizations 

of Ks and ψe and ran 300 simulations of flow and transport 

using STOMP. The resulting sample-ensemble mean concentrations of U-234 at different locations are shown in 

Figure 3-1. Average concentration at the five nodes comprising the well is shown along with plus and minus one 

standard deviation from this average. Maximum well concentrations are also shown. In addition, ensemble mean concentrations at the node immediately above the well and at a 

downstream node at the top of the aquifer are shown. Average well concentration and concentrations at individual 

nodes are also shown in the figure on a log scale to accentuate the differences, particularly between well nodes. 

In Figure 3-1, well concentration is calculated as the arithmetic weighted average of the five well nodes (assuming 

mean uniform mixing within the pumping well). Note that 

ensemble mean concentrations 1 m above the well are 10 

times higher than mean well concentrations. This is because 

this point is in the vadose zone, just above the capillary 

fringe, with a lower water content and no dilution from 

aquifer effects. Note also that the confidence intervals 

(defined here by ± one standard deviation of well concentration) seem relatively narrow due to relatively low concentration variance, implying overall low sensitivity of well 

doses to domain heterogeneity. The maximum well concentration (over 300 simulations) is about double the average 

well concentration. 

Figure 3-2 shows the development of ensemble variance of 

well concentrations at a point in time (and space) as a function of the number of simulations. This is a measure of sample robustness, i.e., it shows how many simulations are 

needed to reach a stable, representative sample variance. 

The figure implies stabilization after 200 simulations. This 

prompted us to reduce the number of simulations in a later 

stage of our numerical experiments (a significant savings of 

computer time).

Figure 3-3 compares the well concentration obtained as the 

ensemble mean of the 300 realizations to the well concentration obtained from a single (deterministic) simulation 

using the scaled mean parameters (Ks* and ψe*), where Ks* 

also equals Kg, the geometric mean. The modest difference 

in peak concentration of only 20% is remarkable. If this difference was in hydraulic heads (rather than concentrations), 

it would be expected to shrink to zero under two-dimensional mean uniform, horizontal flow. In our case, despite 

the three-dimensional, mixed vertical, unsaturated flow and 

transport in the vadose zone and (essentially) horizontal, 

radial flow and transport in the aquifer, the geometric mean 

of saturated hydraulic conductivity is an acceptable rough 

estimate of the effective conductivity of the whole complex 

domain for transport calculations. The scaled mean parameters represent approximate equivalent uniform parameters in 

the sense that the well concentrations in the uniform case 

are an approximation of the ensemble average well concentrations. Calculations with the scaled mean parameters were 

done either as initial tests, or for comparison of numerical 

schemes (e.g., between solvers), or for sensitivity analyses 

(wherever permissible). 

A rigorous analysis of equivalent parameters for this particular transport problem was beyond the scope of this project. 

However, since (a) the parameter distributions were based 

on the scaling assumption, (b) the geometric mean has been 

shown to be the expected effective hydraulic conductivity in 

two-dimensional domains and in slightly heterogeneous 

three-dimensional domains (with relatively small variance) 

(Orr, 1993; Neuman and Orr, 1993)1

, and (c) Figure 3-3 

(and other results) show a fair similarity between mean 

breakthrough curves and breakthrough curves obtained with 

the geometric mean Ks, for the purpose of this study we 

consider the scaled mean parameter values to be acceptable 

equivalent uniform parameter values. Figure 3-3 demonstrates, however, that the use of the scaled mean parameters 

in a deterministic analysis will not completely match the 

ensemble mean well concentration (missing it by 20%) and 

that the deterministic results may not be conservative.

1 Note that the formulations discussed in these publications are based on 

hydraulic responses in fully saturated domains; our case focuses on concentration/doses in a complex unsaturated-saturated domain. 
Base Case Monte Carlo Simulation Results

26

3.2 Dose Distributions and Comparison

The base case RESRAD Monte Carlo simulation also consisted of 300 realizations. As discussed in the previous 

chapter, two sets of soil hydraulic parameter distributions 

were used in the RESRAD simulations. One was generic 

distributions and the other set had several mean and variance values based on the data from the Las Cruces Trench 

site. Because well concentrations are not available from the 

results of a RESRAD Monte Carlo simulation, comparisons 

between the two models were made using dose. Well concentrations from STOMP were converted to total effective 

0

50

100

150

200

250

300

350

0 100 200 300 400 500 600 700 800 900 1000

Time (yr)

Concentration (pCi/L)

c-well

c-above

c+stdev

c-stdev

c-max

0.1

1

10

100

1000

0 100 200 300 400 500 600 700 800 900 1000

Time (yr)

Conc. (pCi/L)

1 node above well

well node 5

well node 4

well node 1

well node 2

well node 3

top of aquifer, downstream

well mean

Figure 3-1. Top: Concentrations at the well with plus/minus one standard deviation, maximum well 

concentration, and concentration in the unsaturated zone immediately above the well, 

Bottom: Concentration on a log scale for the five nodes comprising the well (numbered from 

bottom to top of the aquifer), the node just above the well in the unsaturated zone, a 

downstream node, and the mean well concentration (bottom). Results are ensemble averages 

(except for c-max, c-stdev, and c+stdev) for the base case using 300 realizations.
Base Case Monte Carlo Simulation Results

27

dose equivalent by applying the same drinking water ingestion rate and dose conversion factor for U-234 that were 

used in RESRAD. Limiting the test case simulation to the 

drinking water pathway simplified this comparison. 

Figure 3-4 shows the cumulative distribution functions for 

peak dose and the time of the peak dose for the STOMP 

base case simulation and the two RESRAD simulations. 

Statistics of the peak dose and the time of the peak dose are 

given in Table 3-1 for the three cases. The STOMP simulations result in a significantly smaller variance, and hence, 

smaller uncertainty, than the RESRAD simulations for both 

peak dose and time to peak. This result is not unexpected for 

the generic case since the variance of the RESRAD parame0

10

20

30

40

50

60

70

80

0 50 100 150 200 250 300

Simulations

Var(c-well)

Figure 3-2. Variance of well concentration at 579 years after start of release, as a function of the 

number of simulations

0

5

10

15

20

25

30

0 200 400 600 800 1000

Time (yr)

Well Conc. (pCi/L)

Mean of 300 Realizations

Equivalent Uniform

Figure 3-3. Ensemble mean concentration from the STOMP Base Case Monte Carlo simulation 

compared to the deterministic concentration using equivalent uniform (scaled mean) 

parameter values. The concentration shown is the arithmetic average of the five well nodes.
Base Case Monte Carlo Simulation Results

28

ters represents the variability over a large number of sites 

across the United States, while the covariances of the 

STOMP parameters (Ks and ψe) and resulting concentrations/doses represent the variability at a single site. The 

larger dose variance from the RESRAD simulation reflects 

the intended conservativeness of the generic parameter distributions. For the site-specific RESRAD simulation the 

parameter variances are a mix of generic values and values 

specific to the Las Cruces Trench site. The reduction in the 

site-specific peak dose variance relative to the generic case 

reflects the site-specific Ks data. The lack of a similar 

reduction in the variance of the time of the peak dose primarily reflects the use of generic parameter values (for the 

unsaturated zone bulk density in particular). 

The mean dose predicted from the RESRAD simulations is 

larger than the mean dose from STOMP, reflecting the 

deterministic results presented in Chapter 2. The larger 

mean dose predicted by RESRAD would be expected if the 

RESRAD code represented a conservative analysis of the 

test case. 

1.0

0.8

0.6

0.4

0.2

0

Prob {D<d}

10 20 30 40 50

Peak Dose (mrem/yr)

 Simplified Model/Generic Distributions

 Simplified Model/Site-Specific Distributions

 Complex Model/Spatially Variable Properties

1.0

0.8

0.6

0.4

0.2

0

Prob {T<t}

300 400 500 600 700 800 900 1000

Time of Peak Dose (yr)

 Simplified Model/Generic Distributions

 Simplified Model/Site-Specific Distributions

 Complex Model/Spatially Variable Properties

Figure 3-4. Empirical cumulative distribution functions for the peak dose (top) and time of peak (bottom) over 

300 simulations for the complex model using STOMP and the simplified model using RESRAD. 

Error bars represent the mean value and plus/minus one standard deviation.
Base Case Monte Carlo Simulation Results

29

An unexpected result of the simulations is that the use of the 

site-specific parameter distributions results in a mean dose 

(and a mean time to peak dose) from RESRAD that is less 

similar to the STOMP results than when using generic 

parameter distributions with RESRAD. An explanation for 

this behavior requires a detailed look at the models (see the 

following chapter). In the simplest terms, however, the peak 

dose predicted by RESRAD is primarily a function of the 

hydraulic conductivity in the saturated zone and the time of 

the peak dose is primarily a function of the unsaturated zone 

bulk density (considering only the soil hydraulic parameters 

listed in Tables 2-5 and 2-6). The differences between the 

generic and site-specific mean values for these two parameters explain the differences in generic and site-specific mean 

peak dose and time of peak dose. The relative accuracy of 

the generic results is thus fortuitous. The behavior displayed 

in Figure 3-4 reflects the significant structural differences in 

the simplified and complex models. 

The RESRAD distributions of peak dose are highly skewed 

whereas the peak dose distribution from the STOMP simulation appears closer to normal. The large degree of skewness in the RESRAD results indicates a potential sensitivity 

to parameter values, particularly to the 1D hydraulic conductivity that constitutes the flow path in the simplified 

model, an issue that will be discussed in the following chapter.

The mean dose as a function of time is shown in Figure 3-5 

for the two RESRAD cases (simplified model) and for the 

STOMP base case (complex model). The peaks of these 

mean doses are similar, particularly in magnitude. This 

result contrasts with the means of the peak dose distributions shown in Figure 3-4, which are significantly different 

for the three simulations. For the STOMP model, the peak 

of the mean dose and the mean of the peak dose distribution 

are almost the same because of the smoothness of the breakthrough curve(s). The RESRAD results for the same quantities are very different for this test case because of the 

sharpness of the breakthrough curve, which is a result of the 

lack of dispersion in the model and the form of the contaminant source release function. Differences between the two 

RESRAD cases shown in Figure 3-5 reflect the higher peak 

dose and longer time to peak displayed for these same cases 

in Figure 3-4. Note that doses are underestimated by RESRAD prior to 550-600 years. Note also that similarity of the 

simulated peak mean and mean peak doses implies robustness and reliability of the complex model, in contrast to the 

simplified approach. These results are specific to the test 

case considered and the relationship between the two models’ results for the peak of the mean and the mean of the 

peak doses can be expected to vary with the problem specification, particularly when additional exposure pathways 

are considered. 

Table 3-1. Statistics of peak dose and the time of the peak dose for the simplified model with generic parameter 

distributions, the simplified model with site-specific parameter distributions, and the complex model 

with spatially varying properties

Peak Dose (mrem/yr) Time of Peak Dose (yr)

Generic Site-Specific Complex Generic Site-Specific Complex

Number 300 300 300 300 300 300

Minimum 3.2 12.1 1.73 333 429 481

Maximum 56.5 38.3 10.5 729 1003 591

Mean 15.9 21.6 5.3 513 664 523

Std. Deviation 8.2 4.7 1.6 67 103 18

Median 14.2 21.0 5.2 515 659 520
Base Case Monte Carlo Simulation Results

30

8

6

4

2

0

Mean Dose (mrem/yr)

0 200 400 600 800 1000 1200

Time (yr)

 Simplified Model/Generic Distributions

 Simplified Model/Site-Specific Distributions

 Complex Model/Spatially Variable Properties

Figure 3-5. Mean dose as a function of time for the complex model using STOMP and 

the simplified model using RESRAD
31

4 Sensitivity Cases

An important part of the evaluation of uncertainty estimates 

as predicted by the simplified and complex models is a comparison of the sensitivity to parameter values. Models that 

are highly sensitive to certain parameters are worrisome 

because even small uncertainty in these parameters translates to large uncertainty in predictions. Too much robustness (very low sensitivity) may indicate model weakness or 

a deficiency in model structure, such as over- or underparameterization or strong boundary effects. Since the complex model is more representative of reality, ideally we 

would also like to see that the simplified and complex models are sensitive to the same parameters, thus providing evidence that they represent the flow and transport processes in 

a physically similar way. 

This chapter describes the sensitivity analyses that were carried out with the STOMP and RESRAD models and draws 

some conclusions from the results. 

4.1 Sensitivity to Geostatistical 

Parameters

The geostatistical characterization of the soil hydraulic 

properties used in the complex representation requires 

parameters that have no corollary in the simplified representation. The sensitivity of the STOMP model to these parameters is discussed in this section.

4.1.1 Variance of Ks and ψe

Sensitivity of the STOMP model to the degree of heterogeneity was evaluated by conducting a Monte Carlo simulation for a case with double the variance of log conductivity 

and air entry pressure. This simulation was carried out using 

the same realizations of the scaling factor, α, but translating 

the values to double the variance of Ks and ψe. This translates to fluctuations in conductivity at least an order of magnitude larger than in the base case. Figure 4-1 shows the 

resulting ensemble mean breakthrough curves. The figure 

shows that mean concentrations in the well are slightly 

smaller than the concentrations simulated in the base case. 

The variance of the well concentration increased for the 

high variance case as reflected in the plots of well concentration plus and minus one standard deviation, shown in 

Figure 4-1. The magnitude of the concentration variance 

increase was less significant than anticipated, however. 

Figure 4-2 shows the variance of well concentration at a 

point in time as a function of the number of realizations for 

the high variance case. The well variance is much less stable 

than for the base case due to the higher variability (compare 

to Figure 3-2). Here, 200 simulations are not sufficient. 

Variability of the peak dose and the time of the peak dose 

did not increase significantly, however. These results indicate relatively low sensitivity of the dose to the degree of 

heterogeneity in the subsurface. This conclusion may be the 

result of the small domain used in the simulation, a possibility that was not explored. 

30

25

20

15

10

5

0

Concentration (pCi/l)

0 200 400 600 800 1000

Time (yr)

 Base Case

 Double Variance

Figure 4-1. Ensemble mean concentration in the well (solid lines) and plus/minus one standard 

deviation about the mean well concentration (dashed lines) for the base case and a Monte 

Carlo simulation consisting of 300 realizations with double the variance of log[Ks] and ψe
Sensitivity Cases

32

4.1.2 Correlation Scale

Monte Carlo simulations were also carried out to evaluate 

the sensitivity to correlation scale. Two simulations consisting of 200 realization each were run, one with the correlation scale doubled and one with it cut in half, relative to the 

base case value. Figure 4-3 shows the effect of the horizontal correlation scale on ensemble mean concentrations at the 

pumping well, plus/minus one standard deviation about this 

mean, and the maximum well concentration. Differences in 

the mean well concentration are small. The figure shows, 

however, an increase in the variability/uncertainty of the 

well concentrations as the correlation scale increases, as 

evidenced by the increasing difference between the 

plus/minus one standard deviation curves and between the 

mean and the maximum curves. This result suggests that the 

presence of larger structures within the domain tends to 

increase the well concentration variability.

4.2 Sensitivity to Other Hydrologic 

Parameters

Sensitivity analyses were carried out for a number of parameters that have a representation in both the complex and the 

simplified models. The results for the hydrologic parameters are presented here, with comparisons between the 

STOMP and RESRAD sensitivities. One parameter not considered in the sensitivity analyses was the contaminant distribution coefficient. Previous studies have demonstrated 

the importance of the contaminant distribution coefficient in 

simplified and complex representations. 

For the sake of reducing computer time, a number of the 

STOMP sensitivity analyses presented in this section were 

performed on “equivalent uniform” domains. That is, a 

small number of deterministic simulations were carried out 

with Ks and ψe assigned their scaled mean values and with a 

single parameter value varied incrementally. Based on earlier results (e.g., Figure 3-3), we concluded that the differences between variations in ensemble mean responses and 

variations in responses from an equivalent uniform domain 

would be insignificant for sensitivity analysis purposes, 

with the exception of variations in geostatistical parameters 

discussed in the previous section. 

4.2.1 Recharge Rate

Additional Monte Carlo simulations were carried out to 

evaluate the sensitivity of the STOMP model to the recharge 

rate. Two simulations consisting of 300 realizations each 

were run, one with the recharge rate increased by 40% and 

the other with recharge cut to one-half the base case value. 

Figure 4-4 shows the effect of the higher and lower recharge 

rates on the mean breakthrough curves. As expected, high 

recharge translates to higher concentrations at the well and a 

shorter time to peak. The 40% increase in recharge resulted 

in about a 30% shorter time to peak and about a 20% higher 

peak concentration. Reducing recharge by half (50%) 

resulted in almost doubling the time to peak and about a 

45% reduction in peak concentration. This sensitivity to the 

recharge rate is on a par with that observed in RESRAD 

simulations and reflects the fact that the recharge rate controls the contaminant flux through the unsaturated zone. 

Results from recharge sensitivity calculations for the complex and simplified models are compared in Figure 4-5. The 

complex model results were obtained from the ensemble 

mean well concentrations shown in Figure 4-4. The simplified model results were obtained by running RESRAD with 

the recharge rate as the only variable parameter (other 

parameters were set to the site-specific mean values). The 

effects of recharge on the peak dose and the time of the peak 

0

20

40

60

80

100

120

140

0 50 100 150 200 250 300

Simulations

Well concentration variance

(pCi/L)^2

Figure 4-2. Evolution of the variance of the well concentration at 579 years with the number of Monte Carlo 

realizations - variance of log[Ks] two times the base case value
Sensitivity Cases

33

80

60

40

20

0

Concentration (pCi/L)

0 200 400 600 800 1000

Time (yr)

 Well Concentration

 Well Conc +/- One Std. Dev.

 Maximum Well Concentration

Base Case

Figure 4-3. Effect of correlation scale on breakthrough concentrations at the pumping well. 

Mean concentration, mean +/- one standard deviation, and maximum 

concentration are shown for 200 realizations.

80

60

40

20

0

Concentration (pCi/L)

0 200 400 600 800 1000

Time (yr)

100% Increase in Correlation Scale

80

60

40

20

0

Concentration (pCi/L)

0 200 400 600 800 1000

Time (yr)

50% Decrease in Correlation Scale
Sensitivity Cases

34

60

50

40

30

20

10

0

Concentration (pCi/l)

0 200 400 600 800 1000

Time (yr)

Base Case

60

50

40

30

20

10

0

Concentration (pCi/l)

0 200 400 600 800 1000

Time (yr)

40% Increase in Recharge

60

50

40

30

20

10

0

Concentration (pCi/l)

0 200 400 600 800 1000 1200

Time (yr)

 Well Concentration

 Well Conc. +/- One Std. Dev.

 Maximum Well Concentration

50% Decrease in Recharge

Figure 4-4. Effect of recharge on breakthrough concentrations at the pumping well. Mean 

concentration, mean +/- one standard deviation, and maximum concentration are 

shown for 300 realizations.
Sensitivity Cases

35

dose are shown in Figure 4-5. Other than the previously 

noted differences in magnitudes, the primary difference is 

the greater sensitivity of the RESRAD peak dose results to 

the recharge rate. The slope of the peak dose relationship for 

RESRAD is more than three times the slope for the STOMP 

results. This stems from differences in model structure (one 

vs. three dimensions).

4.2.2 Mean Air Entry Pressure

Sensitivity to the air entry pressure was evaluated by conducting two simulations on uniform domains using air entry 

pressures that were twice and one-half the base case value. 

Concentration results for the two sensitivity cases are shown 

in Figure 4-6. There is little sensitivity to air-entry values. 

This is, again, due to the imposed boundary conditions on 

the top surface. Air-entry variations translate to variations in 

pressure head in the unsaturated zone, but mass flux in the 

vadose zone is fixed by the specified source concentration 

and recharge rate. The air-entry value does not apply to saturated flow in the aquifer. 

4.2.3 Mean Hydraulic Conductivity

Two additional Monte Carlo simulations, each consisting of 

200 realizations, were carried out to examine the sensitivity 

to the mean hydraulic conductivity. One simulation used a 

mean Ks an order of magnitude larger than the base case 

value and one used a mean Ks an order of magnitude 

smaller. Figure 4-7 shows the well concentration results for 

the two simulations with the base case results also shown 

for comparison. Reducing the mean hydraulic conductivity 

by an order of magnitude resulted in more than a five-fold 

increase in the peak concentration and a 20% increase in the 

time of the peak. Increasing the mean conductivity by an 

order of magnitude resulted in nearly a 90% reduction in 

well concentration and a 10% decrease in the time to peak. 

At first glance, the variations in concentrations caused by 

varying mean conductivities seem counter-intuitive, and the 

variations in peak time are much smaller than expected. 

This is apparently due to the boundary conditions imposed 

on the system. In particular, the prescribed water flux and 

source concentrations at the surface imply a steady flux of 

solute mass in the vadose zone regardless of mean Ks. Since 

the majority of the travel time is in the unsaturated zone, 

these boundary conditions limit the sensitivity of the time of 

the peak dose to the mean Ks. On the other hand, with prescribed heads on upstream and downstream boundaries of 

the aquifer, a higher mean conductivity in the aquifer entails 

higher flow rates. This implies greater mixing and flushing 

in the aquifer, which leads to lower concentrations at the 

well. Similarly, lower flow rates in the aquifer result in less 

mixing and hence higher concentrations at the well. 

Note that if the contaminant source boundary condition on 

top was defined as solute mass rather than solute concentrations, the results would be different. In that case, since a 

lower mean Ks implies an apparent increase in water content in the unsaturated zone (in order to maintain the same 

water flux)1

, lower source concentrations would result. By 

the same reasoning, a finite concentration at the surface 

30

25

20

15

10

5

0

Peak Dose (mrem/yr)

0.15 0.20 0.25 0.30 0.35

Recharge (m/yr)

1200

1000

800

600

400

Time of Peak Dose (yr)

 Peak Dose, Simplified Model

 Peak Dose, Complex Model

 Time of Peak, Simplified Model

 Time of Peak, Complex Model

Figure 4-5. Effect of the recharge rate on peak dose and on the time of the peak dose for the simplified (RESRAD) 

and the complex (STOMP) models 
Sensitivity Cases

36

implies a higher solute mass in the vadose zone, which augments the increase in well concentrations caused by the 

lower aquifer flow. Careful examination of the results from 

the mean Ks sensitivity simulations (not shown) indicates 

that concentration variations in the unsaturated zone are 

more pronounced than the concentration variations within 

1 The constant recharge at the surface entails mean steady gravitational 

flow (in all simulations), implying uniform pressure head (plus random 

variations in the random fields) and unsaturated hydraulic conductivity 

approximately equal to the imposed flux.

200

150

100

50

0

Concentration (pCi/L)

0 200 400 600 800 1000

Time (yr)

 Mean Uniform Air Entry

 2x Mean Air Entry

 1/2 Mean Air Entry

Node above well

Top node of well

Downstream in aquifer

Figure 4-6. Effect of air-entry value on breakthrough curves at the unsaturated node above the well, at the top well 

node, and downstream mid-depth in the aquifer. All simulations are single deterministic runs with 

uniform properties.

140

120

100

80

60

40

20

0

Mean Well Concentration (pCi/l)

0 200 400 600 800 1000

Time (yr)

 0.1*Ks

 Base Case

 10*Ks

Figure 4-7. Effect of mean saturated conductivity on the ensemble mean well concentration for 200 

realization Monte Carlo simulations
Sensitivity Cases

37

the nodes comprising the well. These are important differences in model structure (of otherwise the same model) that 

one has to keep in mind when comparing different models.

For the simplified representation, the hydraulic conductivity 

in the aquifer dominated the variability in peak dose. This 

can be seen most clearly by looking at scatter plots derived 

from the RESRAD Monte Carlo simulation results. 

Figure 4-8 (upper plot) shows the peak dose as a function of 

the values of aquifer hydraulic conductivity and unsaturated 

zone bulk density used in the RESRAD simulation. The 

lack of correlation between the dose and the bulk density is 

representative of the scatter plots for the remaining random 

parameters. In fact, virtually all the variability in the peak 

dose was due to the variability in the aquifer hydraulic conductivity. The lower plot in Figure 4-8 shows the same 

RESRAD results for the peak dose as a function of the aquifer hydraulic conductivity and includes the STOMP results 

for comparison. The STOMP results were taken from the 

breakthrough curves presented in Figure 4-7 for the three 

values of mean hydraulic conductivity. Although the range 

of hydraulic conductivities used in the STOMP simulations 

was significantly larger than that used in RESRAD, the 

overall sensitivity of peak dose to hydraulic conductivity 

appears similar for the two models. 

The quantities compared in Figure 4-8 (lower plot) are not 

strictly the same. Each point of the complex model results 

represents an ensemble mean peak dose and an ensemble 

mean hydraulic conductivity. Each point of the simplified 

model curve represents a peak dose and a spatial average 

hydraulic conductivity for a particular realization. Averaging all the results for the simplified model would give a single point equivalent to the central point (of three) from the 

complex model. In order to provide the most direct comparison to the RESRAD sensitivity results shown in Figure 4-8, 

a similar plot was prepared from the results of one of the 

STOMP simulations (the high recharge simulation). The 

spatial geometric mean hydraulic conductivity in the aquifer 

was calculated for each realization. This was then plotted 

versus the peak dose for each realization (calculated from 

average well concentrations). The resulting scatter plot is 

shown in Figure 4-9 and demonstrates that there is no correlation between the spatial geometric mean conductivity of a 

single realization and peak concentration at the pumping 

well as simulated by STOMP. This result is in sharp contrast 

to the RESRAD results. Note that in Figure 4-9 the spatial 

geometric mean hydraulic conductivity varies over approximately the same range as the three ensemble mean values 

used to obtain the STOMP results shown in Figures 4-7 and 

4-8. 

It’s possible that the effect of variation in other parameters 

obscures the relationship between peak dose predicted by 

STOMP and the spatial average aquifer hydraulic conductivity. More likely, in single bounded realizations the generated parameter values (Ks and ψe) tend to be statistically 

anisotropic, implying preferred continuity in certain directions, and hence, preferential flow and transport. Rather 

than an isotropic spatial average of hydraulic conductivities, 

the patterns in which these parameters are distributed in 

each realization determines the transport both in the vadose 

zone and in the aquifer. Orr (1993) showed how these patterns of hydraulic conductivities in single bounded realizations are reflected by anisotropic spatial covariances and 

variograms, implying a non-ergodic behavior and preferred 

directions. Orr (1993) also showed how these irregularities 

are themselves random, disappearing as the number of realizations increases and generated parameter values are averaged over the ensemble of realizations. Indeed, this is why 

generating and simulating a sufficient number of realizations in Monte Carlo simulations is essential for reliable 

analyses. This is consistent with the ensemble averaging we 

perform in order to obtain unbiased estimates of predicted 

values (i.e., the expected values). 

The results shown in Figures 4-8 and 4-9 demonstrate that 

although the complex model ensemble results and the simplified model results behave similarly, the individual realizations from the complex model behave quite differently. 

The results have implications for the value of hydraulic conductivity data. Adopting the homogeneous parameterization 

of RESRAD leads to a conclusion that reducing the uncertainty in the value of the aquifer hydraulic conductivity 

parameter will have a significant impact on the uncertainty 

in peak dose. Looking at the individual realization results 

from the STOMP model suggests that characterization of 

the average aquifer hydraulic conductivity is relatively 

unimportant in reducing uncertainty in peak dose. Sensitivity of the ensemble mean peak dose to the ensemble mean 

hydraulic conductivity is not significant in this case because 

the actual aquifer is a single realization. Characterizing the 

pattern of aquifer heterogeneity is likely to be more important than simply obtaining a value for the spatial mean 

hydraulic conductivity.

4.2.4 Aquifer Hydraulic Gradient

Sensitivity to the aquifer hydraulic gradient was evaluated 

by conducting several additional STOMP simulations using 

the equivalent uniform parameters (scaled mean parameter 

values for Ks and ψe) with the aquifer gradient reduced 

from the base case value of 0.02. Well concentration results 

are shown in the upper plot of Figure 4-10. Peak dose is 

inversely related to the aquifer gradient, increasing approximately seven times with a ten-fold decrease in the aquifer 

gradient. A smaller gradient results in less mixing of uncontaminated groundwater with the contaminated water 

recharging the aquifer from the unsaturated zone. Peak dose 

calculated from these results of the complex model are compared to the equivalent simplified model results in the lower 

plot of Figure 4-10. The relationship between peak dose and 

the aquifer gradient is similar for the two models. Although 
Sensitivity Cases

38

peak dose for the simplified model is significantly higher, 

both models produce about a seven-fold increase in the peak 

dose over the range of aquifer gradients considered. This 

similarity occurs in spite of the fact that RESRAD uses a 

relatively simplified aquifer dilution model.

4.3 Sensitivity to Scenario Parameters

There are a number of parameters related to the exposure 

scenario that might contribute to the uncertainty in the dose. 

A comprehensive study of these parameters deserves further 

consideration, but was not undertaken for this report. We 

examined the sensitivity of the two models to only the depth 

of the well penetration. Our choice of three-dimensional 

modeling for the complex representation was made primarily to allow for realistic representation of the pumping well. 

We thus had an opportunity to compare the effect of well 

depth to the simpler representation used in RESRAD.

35

30

25

20

15

Peak Dose (mrem/yr)

1000 1500 2000 2500

Aquifer Saturated Hydraulic Conductivity (m/yr)

1.4 1.6 1.8 2.0 2.2

Unsaturated Zone Bulk Density (g/cm3)

 Aquifer Ks

 Density

40

30

20

10

0

Peak Dose (mrem/yr)

14x103 2 4 6 8 10 12

Aquifer Saturated Hydraulic Conductivity (m/yr)

 Simplified Model

 Complex Model

Figure 4-8. RESRAD sensitivity to aquifer hydraulic conductivity and unsaturated zone bulk density from 

the base case simulation (top) and a comparison of the RESRAD (simplified) and STOMP 

(complex) results for peak dose sensitivity to aquifer hydraulic conductivity (STOMP results 

from Figure 4-7)
Sensitivity Cases

39

4.3.1 Depth of Well Penetration into the 

Aquifer

The base case results for the STOMP model showed a sharp 

decrease in concentration with depth within the well (see 

Figure 3-1, well nodes 1 to 5). This suggests that the average well concentrations may be sensitive to the depth of 

well penetration. In order to explore this issue, we ran an 

additional Monte Carlo simulation with the STOMP model 

using only the three uppermost well nodes as an active well 

screen, with the same total discharge (or pumping rate). 

Figure 4-11 shows the results with the base case (five well 

nodes) shown for comparison. The ensemble mean well 

concentration is actually slightly lower for the three-node 

well, while the variance of the well concentration increased 

somewhat. Overall, however, the differences in well concentrations are insignificant. This result is not entirely unexpected because ultimately, well concentrations are 

dominated by pumping and ambient flow, two quantities 

that remained the same in the two simulations. 

Results from deterministic STOMP simulations using the 

equivalent uniform parameter values are shown in 

Figure 4-12 (upper plot). Well concentrations are shown for 

well depths of one to five meters (one to five nodes). Sensitivity to the well depth is significant for these simulations, 

with the peak well concentration increasing from 20 pCi/l to 

over 30 pCi/l as the well depth is reduced from 5 to 3 m and 

increasing to more than 50 pCi/l for a well depth of 1 m. 

This sensitivity to the well depth is in contrast to the Monte 

Carlo ensemble mean results presented in Figure 4-11. Peak 

dose results for the deterministic STOMP simulations are 

compared to equivalent results from RESRAD in the lower 

plot of Figure 4-12. The peak dose predicted by RESRAD 

appears to be more sensitive to the depth of the well penetration, particularly for well depths less than 3 m. 

12

10

8

6

4

Peak Dose (mrem/yr)

102 103 104

Geometric Mean Aquifer Hydraulic Conductivity (m/yr)

Figure 4-9. STOMP sensitivity to the spatially averaged (geometric mean) aquifer hydraulic conductivity using 

the high recharge case
Sensitivity Cases

40

120

100

80

60

40

20

0

Well Concentration (pCi/L)

0 200 400 600 800 1000

Time (yr)

Aquifer Gradient

 0.002

 0.008

 0.014

 0.02

120

100

80

60

40

20

0

Peak Dose (mrem/yr)

20x10-3 5 10 15

Aquifer Gradient (m/m)

 Simplified Model

 Complex Model

Figure 4-10. Well concentration over time for aquifer hydraulic gradients of 0.002 to 0.02 from the complex model 

(STOMP) using equivalent uniform hydraulic properties (top) and a comparison of peak dose for the 

complex and simplified models as a function of aquifer hydraulic gradient using equivalent uniform 

properties in both cases (bottom).
Sensitivity Cases

41

50

40

30

20

10

0

Concentration (pCi/L)

0 200 400 600 800 1000

Time (yr)

 Well Mean

 Well + Stdv

 Well - Stdv

 Well Maximum

 Downstream

50

40

30

20

10

0

Concentration (pCi/L)

0 200 400 600 800 1000

Time (yr)

 Well Mean

 Well + Stdv

 Well - Stdv

 Well Maximum

 Downstream

Figure 4-11. STOMP model base case results (top) with a fully penetrating well (5-m well depth) and results from a 

simulation with a 3-m deep well (bottom). Results are shown for the ensemble mean well 

concentration, plus/minus one standard deviation, the maximum well concentration, and a 

downstream aquifer node for a 200 realization Monte Carlo simulation.
Sensitivity Cases

42

50

40

30

20

10

0

Well Concentration (pCi/L)

0 200 400 600 800 1000

Time (yr)

Well Depth

 1 m

 2 m

 3 m

 4 m

 5 m

80

60

40

20

0

Peak Dose (mrem/yr)

1 2 3 4 5

Depth of Well Penetration (m)

 Simplified Model

 Complex Model

Figure 4-12. Well concentration over time for well depths of 1 m to 5 m from the complex model (STOMP) using 

equivalent uniform hydraulic properties (top) and a comparison of peak dose for the complex and 

simplified models as a function of well depth using equivalent uniform properties in both cases 

(bottom).
43

5 Conclusions

This report has evaluated an uncertainty assessment as carried out using RESRAD by comparing the results to similar 

analyses conducted using a more complex code, STOMP. A 

test case was developed for the evaluation using data from 

the Las Cruces Trench site. Characteristics of the test case 

were set to eliminate many differences between the codes 

and to highlight the differences due to the underlying conceptualization and parameterization of subsurface flow and 

transport.

Mean peak dose predicted by the simplified model (RESRAD) for a base case Monte Carlo simulation was several 

times higher than the mean peak dose predicted by the complex model (STOMP). This difference was attributed to the 

lack of dispersion in the RESRAD model and differences in 

the mixing and dilution in the aquifer. For this test case, 

RESRAD provided conservative results in the sense that the 

peak dose was higher than that predicted by STOMP. The 

RESRAD concentration breakthrough curves exhibited a 

sharp peak, however, with essentially no contaminant in the 

well until the time of the peak. This could be interpreted as 

nonconservative behavior since the STOMP simulations 

predicted a much earlier arrival of contaminants and, prior 

to the arrival of the RESRAD peak, the well concentrations 

from the STOMP model were larger than those predicted by 

RESRAD. This point is particularly relevant because the 

regulatory criteria include a time limit: 1000 years from the 

time of decommissioning. If the RESRAD predicted peak 

had occurred after 1000 years, the RESRAD results would 

have been nonconservative in every sense of the word.

The complex model assumed the availability of an extensive 

dataset on which to base the random field characterization 

of the subsurface. Uncertainty in predicted dose was correspondingly small, with the peak dose coefficient of variation being 30%. When the variances of parameters in the 

simplified model were based on a generic dataset, the uncertainty in predicted peak dose was much larger; the coefficient of variation was 52% in this case. When the variances 

of parameters in the simplified model were based on the 

available site-specific data, the coefficient of variation for 

the peak dose was reduced to 22%. In this case, however, 

the mean peak dose was actually less similar to the STOMP 

results than the generic case. This unexpected result 

occurred because the site-specific mean aquifer hydraulic 

conductivity was significantly smaller than the generic 

value and because of the strong dependence of the RESRAD peak dose on aquifer hydraulic conductivity. The lognormal distribution of the hydraulic conductivity 

accentuated the effect.

For the STOMP model, the peak of the mean dose as a function of time and the mean of the peak dose distribution were 

almost the same because of the smoothness of the breakthrough curve(s). The RESRAD results for the same quantities were very different for this test case because of the 

sharpness of the breakthrough curves, which was a result of 

the lack of dispersion in the model. Similarity of the simulated peak mean and mean peak doses implies robustness 

and reliability of the complex model, in contrast to the simplified approach

Sensitivity analyses illustrated the relative importance of 

parameters. For the RESRAD Monte Carlo simulations 

involving random soil hydraulic properties, the variability 

of peak dose was entirely attributed to variability in the 

aquifer hydraulic conductivity. Sensitivity to other parameters was examined by varying one parameter at a time. 

These results indicated that the recharge rate, the aquifer 

gradient, and the depth of penetration of the well would be 

significant contributors to uncertainty in peak dose and also 

the time of the peak dose. 

Sensitivity analyses of the complex model showed that stochastic predictions of mean dose over time exhibited low 

sensitivity to parameter variability (expressed by their variances), including to the variance of log hydraulic conductivity. This may have been due to the use of a relatively small 

simulation domain. Uncertainty in dose predictions, however, was directly (linearly) dependent on the variance of 

hydraulic conductivity. Ensemble mean dose was insensitive to the correlation scales of the hydraulic properties. 

Uncertainty in dose, however, increased somewhat with an 

increase in the correlation scale. Overall, the complex 

model showed robustness with respect to geostatistical 

parameters.

Ensemble mean peak dose predicted by the complex model 

was strongly sensitive to the ensemble mean hydraulic conductivity. Reducing the ensemble mean hydraulic conductivity by an order of magnitude resulted in more than a fivefold increase in the peak dose. Increasing the ensemble 

mean conductivity by an order of magnitude resulted in 

nearly a 90% reduction in peak dose. It was noted, however, 

that for a single ensemble mean hydraulic conductivity, the 

complex model showed no correlation between the spatial 

geometric mean aquifer conductivity of individual realizations and the resulting peak dose. This was in contrast to the 

simplified model in which the dose from individual realizations was strongly correlated with the aquifer conductivity. 

This result has implications for the value of hydraulic conductivity data. Adopting the homogeneous parameterization 

of RESRAD leads to a conclusion that reducing the uncertainty in the value of the aquifer hydraulic conductivity 

parameter will have a significant impact on the uncertainty 

in peak dose. Looking at the individual realization results 

from the STOMP model suggests that characterization of 

the average aquifer hydraulic conductivity is relatively 

unimportant in reducing uncertainty in peak dose. Characterizing the pattern of aquifer heterogeneity is likely to be 

more important than obtaining a value for the spatial mean 

hydraulic conductivity.

RESRAD predicted peak dose was more sensitive to the 

parameter values than was the STOMP predicted peak dose 
Conclusions

44

for the aquifer hydraulic conductivity, aquifer gradient, and 

depth of well penetration. Sensitivity to the recharge rate 

appeared to be comparable for the two codes, perhaps a 

result of the boundary condition constraints imposed on the 

STOMP model. Differences in model sensitivities may lead 

to different conclusions about which parameters are critical 

to the analysis and about the importance of site-specific 

data.

Finally, some observations on computational issues. Flow 

and transport in single realizations of the complex model 

showed strong dependence on the spatial structure of generated hydraulic properties. These large variations in results, 

however, diminished sharply when averaged over the 

ensemble, as the number of realizations exceeded 200. A 

minimum of 200 realizations were thus conducted for each 

of the STOMP Monte Carlo simulations. STOMP model 

results showed some sensitivity to mesh resolution, particularly near the water table, where concentration gradients 

were extremely high. The relatively coarse mesh size was a 

calculated compromise without which the Monte Carlo simulations could not be performed. The variety of numerical 

consideration required by the application of a complex, 

three-dimensional model was clearly illustrated. No such 

considerations were required to execute RESRAD as its 

internal conceptual-mathematical model is fixed. The complex model required many hours of computational time to 

execute a Monte Carlo simulation. RESRAD, on the other 

hand required no more than a few minutes.

We note here that our results are specific to the application 

of RESRAD to the test case we simulated. It may be valuable to conduct similar studies with other multimedia models (e.g. MEPAS/FRAMES, DandD) and for other test 

cases. The complex representation of the test case was less 

complex than reality at many sites. The impact of such processes as preferential flow or chemical heterogeneity were 

not investigated. Our application also did not make use of 

monitoring data and calibration/inversion. A similar test 

case comparison using such data would be instructive.
45

6 References

Bear, J., Hydrauics of Ground Water, McGraw-Hill Book 

Co., New York, NY, 569 pp., 1979.

Beyeler, W. E., P. A. Davis, F. A Duran, M. C. Daily, and T. 

A. Feeney, “Residual Radioactive Contamination from 

Decommissioning: Parameter Analysis,” Draft

NUREG/CR-5512, Vol. 3, Sandia National Laboratories, 

Albuquerque, New Mexico, October, 1999.

Buck, J. W., G. Whelan, J. G. Droppo Jr., D. L. Strenge, K.J. 

Castleton, J. P. McDonald, C. Sato, G. P. Streile, “Multimedia Environmental Pollutant Assessment System (MEPAS) 

Application Guidance,” PNL-10395, Pacific Northwest 

National Laboratory, Richland, WA, 1995.

Carsel, R. F., and R. S. Parrish, “Developing joint probability distributions of soil water retention characteristics,” Wat. 

Resour. Res., 24(5):755-770, 1988.

Dagan, G., and S. P. Neuman (eds.), Subsurface Flow and 

Transport: A Stochastic Approach, Cambridge University 

Press, New York, 1997.

Deutsch, C.V., and A.G. Journel, GSLIB – Geostatistical 

Software Library and User’s Guide, Second Ed., Applied 

Geostatistical Series, Oxford University Press, New York – 

Oxford, 369 pp., 1997.

El-Kadi, A.I (ed.), Groundwater Models for Resources 

Analysis and Management, CRC/Lewis Pub., Boca Raton, 

Florida. 384pp., 1995.

El-Kadi, A.I. and G. Ling, “The Courant and Peclet criteria 

for numerical solution of the Richards equation,” Water 

Resources Research, 29(10):3485-3494, 1993.

Fayer, M.J., J.B. Sisson, W.A. Jordan, A.H. Lu, and P.R. 

Heller, “Subsurface Injection of Radioactive Tracers – Field 

experiment for model validation testing,” NUREG/CR-5996

(PNL-8499), U.S. Nuclear Regulatory Commission, Washington, DC, 1993.

Gelhar, L. W., Stochastic Subsurface Hydrology, PrenticeHall, Englewood Cliffs, New Jersey, 1993.

Gilbert, R. O. Statistical Methods for Environmental Pollution Monitoring, Van Nostrand Reinhold, New York, NY, 

320 pp., 1987.

Kennedy, J.D., and D.L. Strenge, “Residual Radioactive 

Contamination from Decommissioning: Technical Basis for 

Translating Contamination Levels to Annual Total Effective 

Dose Equivalent,” NUREG/CR-5512, (PNL-7994), U.S. 

Nuclear Regulatory Commission, Washington, D.C., 1992.

Kincaid, C.T., M.P. Bergeron, C.R. Cole, M.D. Freshley, 

N.L. Hassig, V.G. Johnson, D.I. Kaplan, R.J. Serne, G.P. 

Streile, D.L. Strenge, P.D. Thorne, L.W. Vail, G.A. Whyatt, 

S.K. Wurstner, “Composite Analysis for Low-Level Waste 

Disposal in the 200 Area Plateau of the Hanford Site,” 

PNNL-11800, Pacific Northwest National Laboratory, Richland, WA, 1998.

Khaleel, R., and E.J. Freeman, “Variability and Scaling of 

Hydraulic Properties for 200 Area Soils, Hanford Site,” 

WHC-EP-0883, UC-900, Westinghouse Hanford Company, 

Richland, WA, 1995.

Leij, F. J., W. J. Alves, M. Th. van Genuchten, and J. R. 

Williams, “The UNSODA Unsaturated Soil Hydraulic 

Database User’s Manual Version 1.0,” EPA/600/R-96/095,

U.S. Environmental Protection Agency, Washington, D.C., 

1996.

Meyer, P.D., M.L. Rockhold, and G.W. Gee, “Uncertainty 

Analyses of Infiltration and Subsurface Flow and Transport 

for SDMP Sites,” NUREG/CR-6565, U.S. Nuclear Regulatory Commission, Washington, DC, 1997.

Meyer, P. D. and G. W. Gee, “Information on Hydrologic 

Conceptual Models, Parameters, Uncertainty Analysis, and 

Data Sources for Dose Assessments at Decommissioning 

Sites,” NUREG/CR-6656, U.S. Nuclear Regulatory Commission, Washington, D.C., 1999.

Meyer, P. D. and R. W. Taira, “Hydrologic Uncertainty 

Assessment for Decommissioning Sites: Hypothetical Test 

Case Applications,” NUREG/CR-6695, U.S. Nuclear Regulatory Commission, Washington, D.C., 2001.

Neuman, S.P., and S. Orr, "Prediction of steady state flow in 

nonuniform geologic media by conditional moments: exact 

nonlocal formalism, effective conductivities, and weak 

approximation'', Water Resources Research, 29(2):341-364, 

1993.

Neuman, S.P., S. Orr, O. Levin, and E. Paleologos, “Theory 

and High Resolution Finite Element Analysis of Two 

Dimensional and Three Dimensional Effective Permeability 

in Strongly Heterogeneous Porous Media,” in Computational Methods in Water Resources IX, Vol. 2: Mathematical 

Modeling in Water Resources, edited by T.F. Russell, R.E. 

Ewing, C.A. Brebbia, W.G. Gray, and G.F. Pinder, Elsevier, 

New York, 1992.

NRC, “Decision Methods for Dose Assessments to Comply 

with Radiological Criteria for License Termination,” Draft 

NUREG-1549, U.S. Nuclear Regulatory Commission, 

Washington, DC, February 3, 1998

NRC, “NMSS Decommissioning Standard Review Plan,” 

NUREG-1727, U.S. Nuclear Regulatory Commission, 

Washington, DC, September, 2000.

Orr, S., Stochastic Approach to Steady State Flow in Nonuniform Geologic Media, Ph.D. Dissertation, Department of 

Hydrology and Water Resources, University of Arizona, 

Tucson, AZ, 356 pp., 1993.
References

46

Press, W.H., B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, Numerical Recipes – The Art of Scientific Computing, 

Cambridge University Press, Cambridge, 818 pp., 1986.

Rockhold, M.L., R.E. Rossi, and R.G. Hills, “Applications 

of similar media scaling and conditional simulations for 

modeling water flow and tritium transport at the Las Cruces 

Trench Site,” Wat. Resour. Res., Vol. 32, No. 3, pp 595-609, 

1996.

Schaap, M. G. and F. J. Leij, “Database-related accuracy and 

uncertainty of pedotransfer functions,” Soil Science,

163(10):765-779, 1998.

Sisson, J.B., and A. Lu, “Field Calibration of Computer 

Models for Applications to Buried Liquid Discharges: A 

Status Report,” RHO-ST-46P, Rockwell Hanford Operations, Richland, WA, 1984.

Smoot, J.L., and R.E. Williams, “A Geostatistical Methodology to Assess the Accuracy of Unsaturated Flow Models,” NUREG/CR-6411 (PNL-10866), Richland, WA, 1996.

Sun, A. Y. and D. Zhang, “Prediction of solute spreading 

during vertical infiltration in unsaturated, bounded heterogeneous porous media,” Water Resour. Res., 36(3 ):715-723, 

2000.

Wierenga, P. J., A. F. Toorman, D. B. Hudson, J. Vinson, M. 

Nash, and R. G. Hills, “Soil Physical Properties at the Las 

Cruces Trench Site,” NUREG/CR-5441, U.S. Nuclear Regulatory Commission, Washington, DC, 1989.

Wierenga, P. J., D. B. Hudson, R. G. Hills, I. Porro, J. Vinson, and M. R. Kirkland, “Flow and Transport Experiments 

at the Las Cruces Trench Site: Experiments 1 and 2,” 

NUREG/CR-5607, U.S. Nuclear Regulatory Commission, 

Washington, DC, 1990.

Yu, C., C. Loureiro, J.-J. Cheng, L. G. Jones, Y. Y. Wang, Y. 

P. Chia, and E. Faillace, “Data Collection Handbook to Support Modeling the Impacts of Radioactive Material in Soil,” 

ANL/EAIS-8, Argonne National Laboratory, Argonne, IL, 

1993.

Zhang, D., Stochastic Methods for Flow in Porous Media: 

Coping with Uncertainties, Academic Press, San Diego, 

California, 2001.

Zhang, D. and Z. Lu, “Stochastic analysis of flow in a heterogeneous unsaturated-saturated system,” Wat. Resour. Res.,

38(2), 10.1029/2001WR000515, 2002.
A-1

Appendix A: Recommended Soil Parameter Distributions

Probability distributions of three types (normal, lognormal, 

and beta) were used to approximate the soil hydraulic 

parameter distributions generated from the Carsel and Parrish (1988) statistics. This appendix provides a summary of 

these distributions and presents tables of recommended distributions for selected soil hydraulic parameters. The information provided in Sections A.1 – A.3 can be found in 

many good probability or statistics textbooks (e.g., Ang and 

Tang, 1975). Definitions of parameters can be found in 

Appendix D.

This appendix was originally presented in NUREG/CR6565 (Meyer et al., 1997).

A.1 The Normal Distribution

The normal distribution has a density function given by

(A-1)

where x is the soil parameter being modeled and and 

are the parameters of the distribution. The mean, μ, and the 

variance, σ2, are related to the parameters of the normal distribution as follows.

(A-2)

(A-3)

Although the normal distribution is unbounded, soil parameters modeled by a normal distribution often have physical 

limits. These limits can be enforced by specifying that the 

soil parameter values fall between given quantiles of the 

normal distribution. In the tables below, the lower (A) and 

upper (B) limits of each normal distribution are the 0.001 

and 0.999 quantiles calculated as follows.

(A-4)

A.2 The Lognormal distribution

The lognormal distribution has a density function given by

(A-5)

where γ and ζ are the parameters of the distribution. The 

mean and variance of the lognormal distribution are related 

to the parameters as follows.

(A-6)

(A-7)

These relationships can also be inverted.

(A-8)

(A-9)

The lognormal distribution is thus completely specified by 

either its parameters or its mean and variance. 

The lognormal distribution is bounded below by zero, but 

has no upper bound. In the tables below, the lower and 

upper bound for the lognormal distributions are the 0.001 

and 0.999 quantiles calculated as follows.

(A-10)

A.3 The Beta Distribution

The beta distribution has a density function given by

(A-11)

where q and r are parameters controlling the shape of the 

distribution and A and B are the lower and upper limits of 

the distribution. β(q,r) is the beta function, calculated 

through its relationship to the gamma function.

(A-12)

where Γ( ) indicates the gamma function.

The mean and variance of the beta distribution are related to 

the parameters as follows.

(A-13)

(A-14)

With some algebraic manipulation, these relationships can 

be inverted to provide the shape parameters as a function of 

the mean, variance, and limits.

(A-15)

(A-16)

The beta distribution can thus be completely specified by its 

lower and upper limits and either its mean and variance or 

its shape parameters.

In the tables below, the lower and upper limits for the beta 

distributions are the actual limits, A and B.

f x( ) 1

2πσ' --------------- 1

2

-- x – μ'

σ' ------------    2

= exp

μ' σ'

μ μ = '

σ

2

σ'

2 =

A = μ' 3.09 – σ' B = μ' 3.09 + σ'

f x( ) 1

2πζx ---------------- 1

2

-- ln( ) x – γ

ζ ---------------------    2

= exp

μ γ 1

2

--ζ

2 = exp( ) +

σ

2

μ

2 ζ2 = [ ] exp 1 ( ) –

γ μ ln 1

2

--ζ

2 = –

ζ σ

2

μ

2

----- 1 +

     = ln

A = exp( ) γ – 3.09ζ B = exp( ) γ + 3.09ζ

f x( ) 1

β( ) q r, ----------------( ) x A– q – 1

( ) B x – r – 1

( ) B A– q r + – 1 = ---------------------------------------------------

β( ) q r, Γ( )Γ q ( )r

Γ( ) q r + = ----------------------

μ A q

q r + = + -----------( ) B A–

σ

2 qr

( ) q r +

2

( ) q r + + 1

--------------------------------------------( ) B A– 2 =

q ( )μ B – μ ( ) – A

σ

2 ----------------------------------- 1 –    μ – A

B A– ------------    =

r q

B – μ

μ – A ------------    =
Recommended Soil Parameter Distributions

A-2

A.4 Recommended Probability 

Distributions for Soil Hydraulic 

Parameters by Soil Texture

Tables A-1 to A-12 contain the recommended distributions 

for the selected soil hydraulic parameters. Each table represents a particular USDA soil textural classification. 

Observed correlations between parameters are given in 

Appendix B. 

 

Table A-1. Recommended distributions for Sand

Parameter Distribution1 Mean Std. Deviation Lower Limit2 Upper Limit2

θs Normal 0.430 0.0600 0.245 0.615

θr LN(-3.09,0.224)* 0.0466 0.0106 0.0228 0.0907

pe Normal 0.383 0.0610 0.195 0.572

fc LN(-2.83,0.241) 0.0607 0.0150 0.0280 0.124

wp LN(-3.09,0.224) 0.0466 0.0106 0.0227 0.0907

awc LN(-4.34,0.387) 0.0141 6.12E-03 3.95E-03 0.0431

α [cm-1] Normal 0.147 0.0255 0.0687 0.226

n LN(0.978,0.0998)* 2.67 0.267 1.95 3.62

hb LN(1.93,0.183) 7.02 1.38 3.92 12.1

λ LN(0.502,0.161) 1.67 0.267 1.00 2.72

b LN(-0.0253,0.216) 0.998 0.226 0.501 1.90

Ks [cm/s] Beta(1.398,1.842) 8.22E-03 4.39E-03 3.50E-04 0.0186

Table A-2. Recommended distributions for Loamy Sand

Parameter Distribution1 Mean Std. Deviation Lower Limit2 Upper Limit2

θs Normal 0.410 0.0900 0.132 0.688

θr Normal 0.0569 0.0145 0.0121 0.102

pe Normal 0.353 0.0913 0.0711 0.635

fc LN(-2.55,0.281) 0.0809 0.0224 0.0327 0.186

wp Normal 0.0570 0.0146 0.0119 0.102

awc LN(-3.85,0.491) 0.0239 0.0125 4.65E-03 0.0966

α [cm-1] Normal* 0.125 0.0404 2.03E-04 0.250

n LN(0.816,0.0910) 2.27 0.209 1.71 3.00

hb LN(2.15,0.401) 9.58 8.59 2.48 29.5

λ LN(0.226,0.164) 1.27 0.209 0.756 2.08

b LN(0.305,0.258) 1.40 0.397 0.610 3.01

Ks [cm/s] Beta(0.7992,1.910) 3.99E-03 3.17E-03 3.90E-05 0.0134
Recommended Soil Parameter Distributions

A-3

Table A-3. Recommended distributions for Sandy Loam

Parameter Distribution1 Mean Std. Deviation Lower Limit2 Upper Limit2

θs Normal 0.410 0.0899 0.132 0.688

θr Beta(2.885,2.304) 0.0644 0.0169 0.0173 0.102

pe Normal 0.346 0.0915 0.0629 0.628

fc LN(-2.21,0.314) 0.116 0.0369 0.0417 0.291

wp Normal 0.0659 0.0179 0.0106 0.121

awc LN(-3.12,0.489) 0.0498 0.0256 9.75E-03 0.200

α [cm-1] Beta(1.816,3.412) 0.0757 0.0368 8.72E-03 0.202

n LN(0.634,0.0818)* 1.89 0.155 1.46 2.43

hb LN(2.71,0.538) 17.7 12.0 2.85 79.4

λ Normal 0.892 0.155 0.412 1.37

b LN(0.632,0.282) 1.96 0.597 0.786 4.50

Ks [cm/s] LN(-7.46,1.33) 1.17E-03 1.37E-03 9.62E-06 0.0347

Table A-4. Recommended distributions for Sandy Clay Loam

Parameter Distribution1 Mean Std. Deviation Lower Limit2 Upper Limit2

θs Normal 0.390 0.0700 0.174 0.606

θr Beta(2.202,2.010) 0.101 6.09E-03 0.0860 0.114

pe Normal 0.289 0.0703 0.0723 0.507

fc LN(-1.59,0.254) 0.212 0.0568 0.0933 0.449

wp LN(-2.14,0.158) 0.120 0.0214 0.0724 0.193

awc Beta(1.890,3.817) 0.0920 0.0393 0.0204 0.237

α [cm-1] LN(-3.04,0.639) 0.0572 0.0337 6.62E-03 0.343

n LN(0.388,0.0858)* 1.48 0.127 1.13 1.92

hb LN(3.04,0.639) 26.2 21.3 2.92 151.

λ Normal 0.479 0.127 0.0865 0.872

b LN(1.41,0.275) 4.27 1.39 1.75 9.57

Ks [cm/s] LN(-9.30,1.75) 3.23E-04 5.98E-04 4.12E-07 0.0202
Recommended Soil Parameter Distributions

A-4

Table A-5. Recommended distributions for Loam

Parameter Distribution1 Mean Std. Deviation Lower Limit2 Upper Limit2

θs Normal 0.430 0.0998 0.122 0.738

θr Beta(3.639,2.652) 0.0776 0.0127 0.0374 0.107

pe Normal 0.352 0.101 0.0414 0.663

fc LN(-1.68,0.300) 0.194 0.0609 0.0735 0.468

wp LN(-2.40,0.250) 0.0935 0.0246 0.0418 0.196

awc LN(-2.40,0.462) 0.101 0.0454 0.0218 0.380

α [cm-1] Beta(1.576,3.625) 0.0367 0.0202 3.51E-03 0.113

n LN(0.442,0.0730) 1.56 0.114 1.24 1.95

hb LN(3.470,0.598) 38.9 29.3 5.05 203.

λ Normal 0.560 0.114 0.209 0.911

b LN(1.08,0.271) 3.07 0.900 1.28 6.82

Ks [cm/s] LN(-9.26,1.66) 2.92E-04 4.91E-04 5.51E-07 0.0159

Table A-6. Recommended distributions for Silt Loam

Parameter Distribution1 Mean Std. Deviation Lower Limit2 Upper Limit2

θs Normal 0.45 0.0800 0.203 0.697

θr Beta(3.349,2.566) 0.0670 0.0142 0.0243 0.0998

pe Normal 0.383 0.0813 0.132 0.634

fc Normal 0.252 0.0776 0.0119 0.491

wp LN(-2.22,0.397) 0.117 0.0471 0.0318 0.368

awc Normal 0.135 0.0402 0.0107 0.259

α [cm-1] LN(-4.10,0.554)* 0.0193 0.0115 2.99E-03 0.0919

n LN(0.343,0.0851) 1.41 0.120 1.08 1.83

hb LN(4.10,0.554) 70.3 41.9 10.9 335.

λ Normal 0.414 0.120 0.0417 0.786

b LN(1.28,0.334) 3.80 1.42 1.28 10.1

Ks [cm/s] LN(-10.4,1.49)* 9.33E-05 2.24E-04 3.12E-07 3.11E-03
Recommended Soil Parameter Distributions

A-5

Table A-7. Recommended distributions for Silt 

Parameter Distribution1 Mean Std. Deviation Lower Limit2 Upper Limit2

θs Normal 0.456 0.110 0.1206 0.799

θr Beta(1.717,1.072) 0.0352 8.97E-03 0.0131 0.0490

pe Normal 0.425 0.110 0.0839 0.766

fc Normal 0.236 0.0578 0.0575 0.415

wp LN(-2.46,0.295) 0.0890 0.0268 0.0342 0.212

awc Normal 0.147 0.0395 0.0252 0.269

α [cm-1] Normal* 0.0178 5.73E-03 3.91E-05 0.0355

n Normal* 1.38 0.0369 1.27 1.49

hb LN(4.10,0.403) 68.1 74.8 17.3 209.

λ Normal 0.380 0.0369 0.266 0.494

b LN(1.16,0.140) 3.21 0.465 2.06 4.89

Ks [cm/s] LN(-10.0,0.475)* 4.89E-05 2.76E-05 9.95E-06 1.87E-04

Table A-8. Recommended distributions for Clay Loam

Parameter Distribution1 Mean Std. Deviation Lower Limit2 Upper Limit2

θs Normal 0.410 0.0900 0.132 0.688

θr Normal 0.0954 9.68E-03 0.0655 0.125

pe Normal 0.315 0.0905 0.0349 0.594

fc LN(-1.27,0.297) 0.292 0.0862 0.112 0.700

wp LN(-1.84,0.257) 0.164 0.0468 0.0714 0.350

awc Beta(2.986,4.318) 0.128 0.0497 9.34E-03 0.301

α [cm-1] LN(-4.22,0.719)* 0.0190 0.0153 1.59E-03 0.136

n Normal 1.32 0.0973 1.02 1.62

hb LN(4.22,0.719) 88.0 71.3 7.37 628.

λ Normal 0.318 0.0973 0.0170 0.618

b LN(1.73,0.323) 5.97 2.37 2.08 15.3

Ks [cm/s] LN(-11.3,2.17) 9.93E-05 2.51E-04 1.42E-08 9.76E-03
Recommended Soil Parameter Distributions

A-6

Table A-9. Recommended distributions for Silty Clay Loam

Parameter Distribution1 Mean Std. Deviation Lower Limit2 Upper Limit2

θs Normal 0.430 0.0699 0.214 0.646

θr Normal* 0.0880 9.00E-03 0.0602 0.116

pe Normal 0.342 0.0705 0.124 0.560

fc Normal 0.347 0.0710 0.127 0.566

wp LN(-1.61,0.233) 0.205 0.0508 0.0970 0.410

awc Normal 0.142 0.0333 0.0387 0.245

α [cm-1] LN(-4.72,0.563) 0.0104 6.08E-03 1.57E-03 0.0508

n Normal* 1.23 0.0610 1.04 1.42

hb LN(4.72,0.563) 132. 81.4 19.7 638.

λ Normal 0.230 0.0610 0.0416 0.418

b LN(1.96,0.265) 7.13 2.34 3.02 15.5

Ks [cm/s] LN(-12.3,1.59) 1.54E-05 3.38E-05 3.44E-08 6.49E-04

Table A-10. Recommended distributions for Sandy Clay

Parameter Distribution1 Mean Std. Deviation Lower Limit2 Upper Limit2

θs Normal 0.380 0.0500 0.226 0.534

θr Beta(4.000,1.487) 0.0993 0.0116 0.0508 0.117

pe Normal 0.281 0.0513 0.122 0.439

fc LN(-1.23,0.210) 0.299 0.0623 0.153 0.559

wp Beta(1.142,4.640) 0.165 0.0344 0.121 0.346

awc Normal 0.134 0.0356 0.0238 0.244

α [cm-1] LN(-3.77,0.562)* 0.0270 0.0164 4.06E-03 0.131

n LN(0.241,0.0653)* 1.28 0.0834 1.04 1.56

hb LN(3.77,0.562) 50.7 30.5 7.64 246.

λ Normal 0.275 0.0834 0.0177 0.533

b LN(1.89,0.260) 6.90 2.27 2.97 14.8

Ks [cm/s] LN(-12.2,2.02)* 3.55E-05 1.48E-04 9.59E-09 2.50E-03
Recommended Soil Parameter Distributions

A-7

Table A-11. Recommended distributions for Silty Clay

Parameter Distribution1 Mean Std. Deviation Lower Limit2 Upper Limit2

θs Normal 0.360 0.0698 0.144 0.576

θr Normal* 0.0706 0.0228 1.47E-04 0.141

pe Normal 0.289 0.0735 0.0623 0.517

fc Normal 0.334 0.0678 0.124 0.543

wp LN(-1.49,0.220) 0.230 0.0512 0.114 0.444

awc Normal 0.103 0.0303 9.63E-03 0.197

α [cm-1] LN(-5.66,0.584)* 4.13E-03 2.60E-03 5.73E-04 0.0211

n LN(0.145,0.0430) 1.16 0.0499 1.01 1.32

hb LN(5.66,0.584) 340. 216. 47.3 1743.

λ Beta(2.591,3.268) 0.157 0.0499 0.0404 0.304

b LN(2.29,0.259) 10.2 2.96 4.43 22.0

Ks [cm/s] LN(-13.9,1.31)* 2.19E-06 4.08E-06 1.64E-08 5.37E-05

Table A-12. Recommended distributions for Clay

Parameter Distribution1 Mean Std. Deviation Lower Limit2 Upper Limit2

θs Normal 0.380 0.0900 0.102 0.658

θr Beta(1.501,1.580) 0.0685 0.0344 8.36E-04 0.140

pe Normal 0.311 0.0963 0.0138 0.609

fc Normal 0.340 0.0893 0.0638 0.615

wp Beta(2.751,4.921) 0.263 0.0770 0.0939 0.567

awc LN(-2.66,0.429) 0.0761 0.0299 0.0186 0.263

α [cm-1] LN(-5.54,0.893) 6.18E-03 7.59E-03 2.50E-04 0.0621

n Beta(0.8857,2.400) 1.13 0.0697 1.04 1.36

hb Beta(0.8002,1.546) 353. 257. 14.1 1007

λ Beta(0.8854,2.399) 0.127 0.0697 0.0397 0.365

b Beta(1.751,11.61) 14.1 6.24 4.93 75.0

Ks [cm/s] LN(-12.36,2.269) 3.65E-05 1.08E-04 3.87E-09 4.76E-03

1. LN(,) = Lognormal(γ,ζ); see Section A.2. Beta(,) =Beta(q,r); see Section A.3.

2. Lower Limit and Upper Limit are 0.001 and 0.999 quantiles for Normal and Lognormal distributions.

* Indicates that the recommended distribution is the same type as used by Carsel and Parrish (1988). This 

applies to the parameters θr

, α, n, and Ks only.
Recommended Soil Parameter Distributions

A-8

A.5 References

Ang, A. H-S. and W. H. Tang, Probability Concepts in 

Engineering Planning and Design, Volume 1, Basic Principles, John Wiley & Sons, New York, 409 pp., 1975.

Carsel, R. F., and R. S. Parrish, “Developing joint probability distributions of soil water retention characteristics,” 

Water Resources Research, 24(5):755-770, 1988.

Meyer, P.D., M.L. Rockhold, and G.W. Gee, “Uncertainty 

Analyses of Infiltration and Subsurface Flow and Transport 

for SDMP Sites,” NUREG/CR-6565, U.S. Nuclear Regulatory Commission, Washington, DC, 1997.
B-1

Appendix B: Soil Parameter Correlation Coefficients

The correlation coefficient is a measure of the strength of a 

linear relationship between two random variables (i.e., soil 

parameters), X and Y. Sample correlation coefficients were 

calculated as follows [e.g., Ang and Tang (1975)].

(B-1)

where

= sample correlation coefficient

xi, yi = sample values for parameters X and Y

= sample mean values calculated as 

(B-2)

sx, sy = sample standard deviations calculated as

(B-3)

N = the number of sample values

Correlations between parameters were induced by applying 

the correlations between θr

, α, n, and Ks given in Carsel and 

Parrish (1988). The rank correlation method of Iman and 

Conover (1982) as embodied in the Latin hypercube sampling code of Iman and Shortencarier (1984) was used. Note 

that the correlations given in the tables below do not necessarily appear to be the same as those of Carsel and Parrish 

(1988) since their correlations were calculated after the 

parameters were transformed to normal distributions. The 

correlations given below were calculated on the untransformed parameters. Definitions of parameters in the tables 

can be found in Appendix D.

This appendix was originally presented in NUREG/CR6565 (Meyer et al., 1997).

 

ρ

ˆ 1

N – 1 ------------

xiyi

i = 1

 – Nxy

sxsy

= ---------------------------------

ρ

ˆ

x y,

x 1

N

--- xi

i = 1

N

= 

sx

2 1

N – 1 ------------ xi ( ) – x 2

i = 1

N

= 

Table B-1. Correlation coefficients for Sand

θs θr pe fc wp awc α n hb λ b Ks

θs 1 -0.01 0.99 0.15 -0.01 0.38 0.00 0.00 -0.02 0.00 -0.29 0.00

θr 1 -0.18 0.94 1 0.59 0.12 -0.84 -0.12 -0.84 0.91 -0.50

pe 1 -0.02 -0.18 0.27 -0.02 0.15 0.01 0.15 -0.44 0.08

fc 1 0.94 0.82 -0.11 -0.91 0.11 -0.91 0.89 -0.67

wp 1 0.59 0.12 -0.84 -0.12 -0.84 0.91 -0.50

awc 1 -0.49 -0.79 0.49 -0.79 0.59 -0.78

α 1 0.29 -0.97 0.29 -0.09 0.73

n 1 -0.28 1 -0.88 0.84

hb 1 -0.28 0.09 -0.68

λ 1 -0.88 0.84

b 1 -0.65

Ks 1
Soil Parameter Correlation Coefficients

B-2

Table B-2. Correlation coefficients for Loamy Sand

θs θr pe fc wp awc α n hb λ b Ks

θs 1 -0.01 0.99 0.27 0.00 0.49 0.00 0.01 -0.03 0.01 -0.50 0.01

θr 1 -0.16 0.85 1 0.34 -0.29 -0.58 0.16 -0.58 0.71 -0.34

pe 1 0.13 -0.16 0.42 0.05 0.10 -0.05 0.10 -0.60 0.07

fc 1 0.85 0.79 -0.53 -0.76 0.33 -0.76 0.57 -0.58

wp 1 0.35 -0.30 -0.58 0.17 -0.58 0.72 -0.35

awc 1 -0.60 -0.68 0.39 -0.68 0.19 -0.63

α 1 0.38 -0.57 0.38 -0.29 0.88

n 1 -0.22 1 -0.64 0.65

hb 1 -0.22 0.17 -0.38

λ 1 -0.64 0.65

b 1 -0.41

Ks 1

Table B-3. Correlation coefficients for Sandy Loam

θs θr pe fc wp awc α n hb λ b Ks

θs 1 0.00 0.98 0.38 0.03 0.54 0.01 0.00 0.01 0.00 -0.44 0.01

θr 1 -0.19 0.72 1 0.34 0.14 -0.79 -0.17 -0.79 0.77 -0.22

pe 1 0.24 -0.16 0.46 -0.02 0.15 0.04 0.15 -0.57 0.05

fc 1 0.78 0.90 -0.35 -0.85 0.35 -0.85 0.51 -0.51

wp 1 0.42 0.08 -0.82 -0.10 -0.82 0.77 -0.25

awc 1 -0.56 -0.65 0.57 -0.65 0.20 -0.56

α 1 0.36 -0.77 0.36 -0.11 0.82

n 1 -0.28 1 -0.78 0.60

hb 1 -0.28 0.05 -0.51

λ 1 -0.78 0.60

b 1 -0.33

Ks 1
Soil Parameter Correlation Coefficients

B-3

Table B-4. Correlation coefficients for Sandy Clay Loam

θs θr pe fc wp awc α n hb λ b Ks

θs 1 0.00 1 0.48 0.21 0.58 -0.01 0.00 -0.02 0.00 -0.43 -0.01

θr 1 -0.09 -0.02 0.23 -0.16 0.37 -0.11 -0.36 -0.11 0.21 0.16

pe 1 0.48 0.19 0.59 -0.04 0.01 0.01 0.01 -0.45 -0.03

fc 1 0.88 0.97 -0.67 -0.81 0.66 -0.81 0.42 -0.50

wp 1 0.73 -0.51 -0.81 0.68 -0.81 0.65 -0.33

awc 1 -0.69 -0.73 0.58 -0.73 0.24 -0.54

α 1 0.77 -0.70 0.77 -0.49 0.82

n 1 -0.65 1 -0.76 0.71

hb 1 -0.65 0.57 -0.39

λ 1 -0.76 0.71

b 1 -0.38

Ks 1

Table B-5. Correlation coefficients for Loam

θs θr pe fc wp awc α n hb λ b Ks

θs 1 0.00 0.99 0.55 0.18 0.63 0.03 0.00 0.04 0.00 -0.46 0.03

θr 1 -0.13 0.29 0.79 -0.03 -0.04 -0.70 0.07 -0.70 0.67 0.14

pe 1 0.50 0.08 0.63 0.03 0.09 0.03 0.09 -0.54 0.01

fc 1 0.75 0.93 -0.63 -0.71 0.70 -0.71 0.28 -0.41

wp 1 0.47 -0.42 -0.87 0.56 -0.87 0.69 -0.16

awc 1 -0.62 -0.49 0.63 -0.49 0.00 -0.46

α 1 0.60 -0.73 0.60 -0.37 0.82

n 1 -0.55 1 -0.79 0.41

hb 1 -0.55 0.39 -0.42

λ 1 -0.79 0.41

b 1 -0.21

Ks 1
Soil Parameter Correlation Coefficients

B-4

Table B-6. Correlation coefficients for Silt Loam

θs θr pe fc wp awc α n hb λ b Ks

θs 1 -0.01 0.98 0.48 0.20 0.70 -0.02 -0.01 -0.02 -0.01 -0.20 -0.02

θr 1 -0.18 0.50 0.66 0.18 -0.29 -0.59 0.27 -0.59 0.63 -0.25

pe 1 0.39 0.08 0.66 0.03 0.10 -0.06 0.10 -0.31 0.03

fc 1 0.91 0.87 -0.72 -0.80 0.73 -0.80 0.63 -0.45

wp 1 0.58 -0.63 -0.89 0.73 -0.89 0.86 -0.36

awc 1 -0.65 -0.50 0.54 -0.50 0.20 -0.44

α 1 0.74 -0.75 0.74 -0.56 0.80

n 1 -0.69 1 -0.88 0.48

hb 1 -0.69 0.68 -0.39

λ 1 -0.88 0.48

b 1 -0.31

Ks 1

Table B-7. Correlation coefficients for Silt

θs θr pe fc wp awc α n hb λ b Ks

θs 1 -0.02 1 0.90 0.52 0.97 0.00 0.00 0.03 0.00 -0.39 0.02

θr 1 -0.10 0.25 0.57 -0.02 -0.19 -0.60 0.04 -0.60 0.70 -0.21

pe 1 0.88 0.48 0.96 0.01 0.05 0.03 0.05 -0.44 0.04

fc 1 0.81 0.92 -0.35 -0.37 0.16 -0.37 -0.03 -0.30

wp 1 0.51 -0.60 -0.72 0.48 -0.72 0.44 -0.45

awc 1 -0.10 -0.06 -0.09 -0.06 -0.35 -0.14

α 1 0.55 -0.49 0.55 -0.41 0.89

n 1 -0.20 1 -0.84 0.44

hb 1 -0.20 0.13 -0.29

λ 1 -0.84 0.44

b 1 -0.34

Ks 1
Soil Parameter Correlation Coefficients

B-5

Table B-8. Correlation coefficients for Clay Loam

θs θr pe fc wp awc α n hb λ b Ks

θs 1 0.00 0.99 0.65 0.38 0.76 -0.02 0.00 -0.04 0.00 -0.40 0.01

θr 1 -0.11 -0.50 -0.46 -0.43 0.73 0.58 -0.74 0.58 -0.35 0.51

pe 1 0.69 0.43 0.80 -0.10 -0.06 0.04 -0.06 -0.36 -0.04

fc 1 0.89 0.90 -0.60 -0.71 0.55 -0.71 0.23 -0.42

wp 1 0.60 -0.55 -0.84 0.75 -0.84 0.57 -0.33

awc 1 -0.52 -0.45 0.26 -0.45 -0.13 -0.42

α 1 0.79 -0.62 0.79 -0.42 0.89

n 1 -0.80 1 -0.73 0.58

hb 1 -0.80 0.70 -0.36

λ 1 -0.73 0.58

b 1 -0.26

Ks 1

Table B-9. Correlation coefficients for Silty Clay Loam

θs θr pe fc wp awc α n hb λ b Ks

θs 1 0.00 0.99 0.73 0.46 0.85 -0.01 0.00 0.00 0.00 -0.19 -0.03

θr 1 -0.13 -0.42 -0.46 -0.21 0.72 0.55 -0.71 0.55 -0.37 0.47

pe 1 0.77 0.51 0.87 -0.10 -0.07 0.10 -0.07 -0.15 -0.09

fc 1 0.90 0.75 -0.62 -0.65 0.58 -0.65 0.35 -0.45

wp 1 0.40 -0.68 -0.84 0.79 -0.84 0.69 -0.42

awc 1 -0.29 -0.11 0.02 -0.11 -0.30 -0.32

α 1 0.86 -0.75 0.86 -0.57 0.83

n 1 -0.84 1 -0.82 0.60

hb 1 -0.84 0.80 -0.41

λ 1 -0.82 0.60

b 1 -0.32

Ks 1
Soil Parameter Correlation Coefficients

B-6

Table B-10. Correlation coefficients for Sandy Clay

θs θr pe fc wp awc α n hb λ b Ks

θs 1 0.00 0.97 0.58 0.35 0.68 0.02 0.00 0.01 0.00 -0.23 0.05

θr 1 -0.23 -0.70 -0.82 -0.42 0.75 0.88 -0.92 0.88 -0.78 0.28

pe 1 0.72 0.53 0.75 -0.15 -0.20 0.22 -0.20 -0.05 -0.02

fc 1 0.89 0.89 -0.70 -0.78 0.68 -0.78 0.44 -0.33

wp 1 0.58 -0.67 -0.85 0.87 -0.85 0.74 -0.24

awc 1 -0.58 -0.55 0.35 -0.55 0.05 -0.35

α 1 0.87 -0.74 0.87 -0.59 0.58

n 1 -0.86 1 -0.79 0.44

hb 1 -0.86 0.83 -0.26

λ 1 -0.79 0.44

b 1 -0.23

Ks 1

Table B-11. Correlation coefficients for Silty Clay

θs θr pe fc wp awc α n hb λ b Ks

θs 1 0.00 0.95 0.94 0.74 0.84 -0.01 -0.01 -0.01 -0.01 -0.31 0.02

θr 1 -0.31 -0.29 -0.49 0.19 0.89 0.79 -0.88 0.79 -0.46 0.64

pe 1 0.98 0.86 0.74 -0.28 -0.25 0.26 -0.25 -0.15 -0.18

fc 1 0.91 0.70 -0.32 -0.33 0.24 -0.33 -0.07 -0.24

wp 1 0.34 -0.50 -0.64 0.52 -0.64 0.32 -0.34

awc 1 0.14 0.33 -0.36 0.33 -0.70 0.03

α 1 0.84 -0.72 0.84 -0.47 0.86

n 1 -0.78 1 -0.77 0.64

hb 1 -0.78 0.63 -0.44

λ 1 -0.77 0.64

b 1 -0.31

Ks 1
Soil Parameter Correlation Coefficients

B-7

B.1 References

Ang, A. H-S. and W. H. Tang, Probability Concepts in 

Engineering Planning and Design, Volume 1, Basic Principles, John Wiley & Sons, New York, 409 pp., 1975.

Carsel, R. F., and R. S. Parrish, “Developing joint probability distributions of soil water retention characteristics,” 

Water Resources Research, 24(5):755-770, 1988.

Iman, R. L. and W. J. Conover, “A distribution-free 

approach to inducing rank correlation among input variables,” Communications in Statistics, B11(3):311-334, 

1982.

Iman, R. L. and M. J. Shortencarier, “A FORTRAN 77 Program and User’s Guide for the Generation of Lain Hypercube and Random Samples for Use with Computer 

Models,” NUREG/CR-3624 (also SAND83-2365), U.S. 

Nuclear Regulatory Commission, Washington, D.C., 1984.

Meyer, P.D., M.L. Rockhold, and G.W. Gee, “Uncertainty 

Analyses of Infiltration and Subsurface Flow and Transport 

for SDMP Sites,” NUREG/CR-6565, U.S. Nuclear Regulatory Commission, Washington, DC, 1997.

Table B-12. Correlation coefficients for Clay

θs θr pe fc wp awc α n hb λ b Ks

θs 1 0.00 0.93 0.88 0.73 0.76 0.00 -0.01 0.00 -0.01 -0.26 -0.01

θr 1 -0.36 -0.38 -0.50 0.13 0.70 0.79 -0.85 0.79 -0.52 0.53

pe 1 0.96 0.85 0.66 -0.25 -0.29 0.31 -0.29 -0.06 -0.20

fc 1 0.95 0.55 -0.38 -0.45 0.36 -0.45 0.08 -0.30

wp 1 0.25 -0.47 -0.63 0.57 -0.63 0.33 -0.32

awc 1 0.09 0.28 -0.38 0.28 -0.63 -0.06

α 1 0.82 -0.61 0.82 -0.46 0.86

n 1 -0.78 1 -0.73 0.64

hb 1 -0.78 0.67 -0.37

λ 1 -0.73 0.64

b 1 -0.30

Ks 1

C-1

Appendix C: Recommended Soil Bulk Density Distributions

Dry soil bulk density data were obtained from the U.S. Natural Resources Conservation Service Soil Characterization 

Database, dated May 1994. The data were divided according to USDA soil textural class based on the sand, silt, and 

clay percentages. The distribution of these data over textural 

classes can be seen in Figure 5-5 of NUREG/CR-6656 

(Meyer and Gee, 1999).

For each textural class, the Kolmogorov-Smirnov D-statistic was calculated using hypothetical normal and lognormal 

distributions. The parameters for the hypothetical normal 

and lognormal distributions were based on the data. In all 

cases the D-statistic from the normal distribution was 

smaller than that from the lognormal distribution. In addition, plots of the bulk density histogram for each soil texture 

were examined and appeared to better fit a normal distribution. For these reasons, the normal distribution is recommended for modeling bulk density. 

The mean and standard deviation of the bulk density data 

are given in the table below. These values reflect the elimination of outliers from the dataset for each textural class. 

Outliers were defined as those points outside the mean plus 

or minus four times the standard deviation, where the mean 

and standard deviation were calculated as the sample mean 

and sample standard deviation with the potential outliers 

eliminated. The upper and lower limits given in the table are 

the actual limits of the data extracted from the Soil Characterization Database. Samples with a bulk density less than 

0.5 g/cm3 or greater than 2.0 g/cm3 were not included in the 

analysis.

Bulk density is highly correlated to saturated water content 

and effective porosity with correlations generally between 

-0.95 and -0.99.

C.1 References

Meyer, P.D. and G.W. Gee, “Information on Hydrologic 

Conceptual Models, Parameters, Uncertainty Analysis, and 

Data Sources for Dose Assessments at Decommissioning 

Sites,” NUREG/CR-6656, U.S. Nuclear Regulatory Commission, Washington, DC, 1999.

Table C-1. Recommended parameters of normal distributions for bulk density

Soil Texture Number of Samples Mean Std. Deviation Lower Limit Upper Limit

Sand 811 1.578 0.158 1.0 1.99

Loamy Sand 1889 1.515 0.262 0.5 2.0

Sandy Loam 7195 1.461 0.268 0.5 2.0

Sandy Clay Loam 2189 1.518 0.186 0.76 2.0

Loam 5198 1.418 0.240 0.54 1.99

Silt Loam 6411 1.366 0.227 0.5 1.99

Silt 195 1.330 0.202 0.63 1.74

Clay Loam 3396 1.410 0.197 0.63 2.0

Silty Clay Loam 3139 1.405 0.148 0.8 1.91

Sandy Clay 386 1.491 0.177 0.87 1.94

Silty Clay 2165 1.37 0.154 0.73 1.82

Clay 4539 1.292 0.177 0.58 1.9

D-1

Appendix D: Summary of Water Retention and Conductivity Models

This appendix defines the parameters appearing in the tables 

of Appendices A and B. The information in this appendix 

was taken from NUREG/CR-6565 (Meyer et al., 1997).

Richards equation (Richards, 1931) forms the basis for most 

process-based descriptions of water movement in the unsaturated zone. Richards equation can be expressed as

(D-1)

where

θ = volumetric water content, or volume of water 

per unit bulk volume of soil,

h = soil water tension, h ≥ 0

z = depth, measured positive downward from the 

soil surface,

t = time, 

K(h) = hydraulic conductivity, and

u = a source or sink term used to account for 

water uptake by plant roots.

To solve Richards equation, constitutive functions relating 

the unsaturated hydraulic conductivity and the water content 

to the pressure head are needed. The most commonly used 

relationships are those of van Genuchten (1980), Brooks 

and Corey (1964), and Campbell (1974), although other 

expressions are available (Mualem, 1992; Rossi and 

Nimmo, 1994; Fayer and Simmons, 1995). 

D.1 Van Genuchten Model

The van Genuchten water retention relationship is

(D-2)

where 

Se = effective saturation = , 0 ≤ Se ≤ 1

α = curve fitting parameter related to air entry 

pressure

n, m = curve fitting parameters related to pore size 

distribution; the relationship, m=1-1/n, is 

often assumed

θr = residual water content

θs = saturated water content 

The van Genuchten hydraulic conductivity relationship, 

based on the hydraulic conductivity model of Mualem 

(1976) is

. (D-3)

or

(D-4)

where

Ks = saturated hydraulic conductivity

D.2 Brooks-Corey Model

The Brooks-Corey water retention relationship is

 for h ≥ hb (D-5)

 otherwise. (D-6)

When combined with the relative permeability model of 

Burdine (1953), Brooks and Corey derived the following 

hydraulic conductivity relationship. 

(D-7)

or

 for h ≥ hb (D-8)

and otherwise. (D-9)

where 

hb = curve fitting parameter related to air entry 

pressure

λ = curve fitting parameter related to pore size 

distribution.

Carsel and Parrish (1988) used the following equivalence 

between the Brooks-Corey and van Genuchten parameters:

hb = α-1 and λ = n - 1.

D.3 Campbell Model

Campbell (1974) adopted a water retention relationship similar to Brooks and Corey’s, but with θr = 0. 

 for h ≥ hb (D-10)

 otherwise. (D-11)

Note that because θr = 0, b ≠ 1/λ. Campbell (1974) derived 

a corresponding hydraulic conductivity relationship. 

(D-12)

∂t

∂θ

∂z

∂ –K( ) h ∂z

∂h = – K h( ) + u

Se( ) h 1 ( ) αh n [ ] + –m =

θ θr –

θs θr – ---------------

K Se ( ) Ks Se 1 1 Se

1 ⁄ m ( ) – m [ ] – 2

=

K h( ) Ks

1 ( ) αh n – 1 1 ( ) αh n [ ] + –m{ } – 2

1 ( ) αh n [ ] +

0.5m = ------------------------------------------------------------------------

Se( ) h

hb

h

----    λ

=

Se( ) h = 1

K Se ( ) Ks Se ( )3 2 + ⁄ λ =

K h( ) Ks

hb

h

----    2 3 + λ

=

K h( ) = Ks

θ

θs

---- hb

h

----    1 ⁄ b

=

θ

θs

---- 1 =

K( ) θ Ks

θ

θs

----    2b + 3

=
Summary of Water Retention and Conductivity Models

D-2

or

 for h ≥ hb (D-13)

and otherwise. (D-14)

where 

b = curve fitting parameter related to pore size 

distribution.

Note that all of these single-valued relationships (Equations 

D-2 through D-14) assume that hysteresis is not important.

D.3.1 Calculation of Campbell’s b Parameter

An expression for b in terms of θs, θr, and λ is derived by 

assuming that the Brooks-Corey model (Equation D-7) and 

the Campbell model (Equation D-12) predict the same 

hydraulic conductivity for a given value of water content. In 

this case, the water content used is that corresponding to an 

effective saturation of 0.5. Assuming Se = 0.5 and using the 

definition of effective saturation given above, it follows that

(D-15)

Substituting this expression in Equation D-12 and equating 

this with Equation D-7 leads to

(D-16)

Equation D-16 can be solved for b,

(D-17)

D.4 Additional Parameters

Several additional soil hydraulic parameters may be 

required by dose assessment codes. These parameters and 

the methods by which they were calculated are discussed 

here.

• Effective porosity, pe = θs - θr

• Field capacity, fc = θ(K = 10-8 cm/s)

Field capacity is generally interpreted as the water content at which drainage from a field soil becomes negligible (see the discussion by Hillel, 1980). Field capacity is 

often calculated as the water content at a specified tension, usually taken to be 340 cm (1/3 bar). Hillel (1980) 

argues, however, that the field capacity should be based 

on the drainage rate considered negligible (which is a 

function of the intended application). Field capacity was 

calculated here as the water content at which the unsaturated hydraulic conductivity equals 10-8 cm/s using the 

K h( ) Ks

hb

h

----    2 3 + ⁄ b

=

K h( ) = Ks

θ

θs

---- 0.5 θs θr ( )θ – + r

θs

--------------------------------------- 0.5 1 θr θs = = ( ) + ⁄

0.5 1 θr θs ( ) ( ) + ⁄ 3 2 + b 0.53 2 + ⁄ λ =

b 0.5 ln( ) 0.5 ( ) 3 2 + ⁄ λ

0.5 1 θr θs ln[ ] ( ) + ⁄

---------------------------------------------- 3 –

  	 
   =

van Genuchten model (10-8 cm/s ≅ 3 mm/yr).The value 

of 10-8 cm/s was chosen because it represents a water 

flux at which contaminant transport is likely to be insignificant. This value also results in somewhat larger field 

capacity values and a more realistic available water 

capacity for very coarse textured soils than using the 

water content at 1/3 bar soil pressure. See Meyer and 

Gee (1999) for a more detailed discussion of field capacity. 

• Wilting point, wp = θ(h = 15,300 cm)

Wilting point is the minimum water content (or maximum tension) at which plants can extract water from the 

soil. Wilting point was calculated as the water content at 

a tension of 15,300 cm (15 bars). 

• Available water capacity, awc = fc - wp

Available water capacity represents the amount of water 

available for plant uptake.

D.5 References

Brooks, R. H., and A. T. Corey, “Hydraulic Properties of 

Porous Media,” Hydrology Paper 3, Colorado State Univ., 

Ft. Collins, CO, 1964.

Burdine, N. T. “Relative permeability calculations from 

pore-size distribution data.” Petr. Trans., Am. Inst. Mining 

Metall. Eng., 198:71-77, 1953.

Campbell, G. S., “A simple method for determining unsaturated conductivity from moisture retention data,” Soil Sci.,

117:311-314, 1974.

Fayer, M. J., and C. S. Simmons, “Modified soil water 

retention functions for all matric suctions,” Water Resources 

Research, 31:1233-1238, 1995.

Meyer, P.D., M.L. Rockhold, and G.W. Gee, “Uncertainty 

Analyses of Infiltration and Subsurface Flow and Transport 

for SDMP Sites,” NUREG/CR-6565, U.S. Nuclear Regulatory Commission, Washington, DC, 1997.

Meyer, P. D. and G. W. Gee, “Flux-based estimation of field 

capacity,” Journal of Geotechnical and Geoenvironmental 

Engineering, 125(7):595-599, 1999.

Mualem, Y., “A new model for predicting the hydraulic 

conductivity of unsaturated porous media,” Water 

Resources Research, 12(3):513-522, 1976.

Mualem, Y., “Modeling the Hydraulic Conductivity of 

Unsaturated Porous Media,” in Indirect Methods for Estimating the Hydraulic Properties of Unsaturated Soils, M. 

Th. van Genuchten, F. J. Leij, and L. J. Lund (eds.), University of California, Riverside, CA, pp. 15-36, 1992.

Richards, L. A., “Capillary conduction of liquids through 

porous mediums,” Physics, 1, 318-333, 1931.
Summary of Water Retention and Conductivity Models

D-3

Rossi, C. and J. R. Nimmo, “Modeling of soil water retention from saturation to oven dryness,” Water Resources 

Research, 30(3):701-708, 1994.

van Genuchten, M. Th., “A closed-form equation for predicting the hydraulic conductivity of unsaturated soils,” Soil 

Sci. Soc. Am. J., 44:892-898, 1980.