Geostatistical Modeling of Sediment Abundance in a 

Heterogeneous Basalt Aquifer at the 

Idaho National Laboratory, Idaho

Scientific Investigations Report 2006-5316

Prepared in cooperation with the U.S. Department of Energy

U.S. Department of the Interior

U.S. Geological Survey

DOE/ID-22201

Geostatistical Modeling of Sediment 

Abundance in a Heterogeneous Basalt 

Aquifer at the Idaho National Laboratory, 

Idaho

By John A. Welhan, Idaho Geological Survey, 

Renee L. Farabaugh, Department of Geosciences, Idaho State University, 

Melissa J. Merrick, Department of Geosciences, Idaho State University, and 

Steven R. Anderson, formerly with U.S. Geological Survey 

Prepared in cooperation with the 

U.S. Department of Energy DOE/ID-22201

Scientific Investigations Report 2006–5316

U.S. Department of the Interior

U.S. Geological Survey
U.S. Department of the Interior

DIRK KEMPTHORNE, Secretary

U.S. Geological Survey

Mark D. Myers, Director

U.S. Geological Survey, Reston, Virginia: 2007

For product and ordering information:

World Wide Web: http://www.usgs.gov/pubprod

Telephone: 1-888-ASK-USGS

For more information on the USGS--the Federal source for science about the Earth, its natural and living resources, 

natural hazards, and the environment:

World Wide Web: http://www.usgs.gov

Telephone: 1-888-ASK-USGS

Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the 

U.S. Government.

Although this report is in the public domain, permission must be secured from the individual copyright owners to 

reproduce any copyrighted materials contained within this report.

Suggested citation:

Welhan, J.A., Farabaugh, R.L., Merrick, M.J., and Anderson, S.R., 2007, Geostatistical modeling of sediment abundance in a heterogeneous basalt aquifer at the Idaho National Laboratory, Idaho: U.S. Geological Survey Scientific 

Investigations Report 2006-5316, 32 p.
iii

Contents

Abstract...........................................................................................................................................................1

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

Purpose and Scope .............................................................................................................................3

Previous Investigations........................................................................................................................3

Interpreted Stratigraphy.............................................................................................................3

Composited Stratigraphy............................................................................................................3

Subsurface Geostatistical Studies ...........................................................................................5

Geohydrologic Setting...................................................................................................................................6

Geostatistical Methods.................................................................................................................................6

Data Selection ......................................................................................................................................6

Data Preparation...................................................................................................................................7

Modified Coordinates..................................................................................................................7

Duplicate Coordinates ................................................................................................................7

Ambiguous Stratigraphic Designations ...................................................................................7

Database Organization ...............................................................................................................7

Analysis Approach................................................................................................................................8

Kriging Approach..................................................................................................................................8

Software Selection...............................................................................................................................8

Geostatistical Analysis and Modeling........................................................................................................9

Exploratory Data Analysis ..................................................................................................................9

General Considerations: Spatial Clustering and Partial Penetration..................................9

Sediment Abundance by Composite Unit ................................................................................9

Sediment Abundance by Model Layer...................................................................................13

Conclusions Derived from Exploratory Analysis ..................................................................16

Autocorrelation Analysis ..................................................................................................................17

Variography of Composite Units..............................................................................................17

Indicator Variography of Model Layers .................................................................................20

Kriging Approach.......................................................................................................................22

Kriging Results.....................................................................................................................................23

Multiple Indicator Kriging.........................................................................................................23

Geohydrologic Implications .....................................................................................................24

Application to Scaling of Aquifer Hydraulic Conductivity ..................................................25

Summary and Conclusions.........................................................................................................................28

Acknowledgments.......................................................................................................................................29

References Cited..........................................................................................................................................29
iv

Figures

Figure 1. Map showing study area, major topographic features, ground-water flow

model extent, and locations of borehole data used for this analysis in and

near the Idaho National Laboratory, Idaho ............................................. 2

Figure 2. Map showing principal volcanic, structural, and geologic features in the

study area, showing the inferred location of the Big Lost Trough in relation to

major volcanic constructional zones...................................................... 4

Figure 3. Diagram showing cross section A-A’ showing composite stratigraphic units 1

through 14, ground-water model layers A through F, and the relation between

model layers and actual stratigraphy ................................................... 5

Figure 4. Maps showing geographic distribution of sediment-rich and sediment-poor

boreholes that fully penetrate composite stratigraphic units 1 through 7 and 8

through 14, as defined by Anderson and Liszewski (1997).............................. 10

Figure 5. Graphs showing comparative distributions of sediment abundances showing

the effect of spatial clustering on summary statistics ................................. 12

Figure 6. Box plots distribution of thickness of total accumulated sediment in

boreholes that fully penetrate individual composite units 1 through 7 and in

boreholes that penetrate the oldest units 8 to 14 combined ........................... 13

Figure 7. Histograms showing sediment abundance in each of composite units 1

through 7 and in combined unit 8 to 14, as a percentage of composite unit

thickness .................................................................................... 14

Figure 8. Histograms showing distributions of sediment percentage in various

groupings of composite units showing evidence of statistical similarity

among the youngest (other than composite unit 1) and the oldest .................. 14

Figure 9. Histograms showing distributions of total sediment thickness among various

groupings of composite units ............................................................ 15

Figure 10. Maps showing availability of boreholes that fully penetrate each model layer

and that were used to model the two-dimensional spatial abundance of

sediment in each model layer ............................................................ 16

Figure 11. Histograms showing distributions of sediment percentages in boreholes that

fully penetrate layers A, B, and C, showing similar distributions, and in the

unsaturated zone (layer U), which has a substantially different distribution

because of the influence of composite unit 1 .......................................... 17

Figure 12. Histograms showing distributions of sediment percentage from fully

penetrating wells in model layers A, B, and C in different areas north and

south of the stratigraphic discontinuity and among clustered and

nonclustered areas shown in figure 10 ................................................... 18

Figure 13. Relative variograms of sediment thickness within composite units 1 through

6, based on the inverted correlogram statistic .......................................... 19

Figure 14. Graphs showing cumulative frequency distributions (CFDs) of sediment

percentage in boreholes that fully penetrate layers A, B, and C showing CFD

probabilities calculated from equation 1 ................................................ 20
Figure 15. Graphs showing examples of multiple indicator semivariograms at selected

thresholds for layer A, showing the nature of the autocorrelation structure and 

the progressive change in nuggets, sills, and ranges as thresholds increase ...... 21

Figure 16. Graph showing results of multiple indicator variogram analysis for layer A,

showing the progressive change in autocorrelation structure across a range

of indicator thresholds ..................................................................... 22

Figure 17. Maps showing median percentage of sediment derived from multiple

indicator kriging of layers A, B, and C. Percentages are equivalent to

thickness, in feet, within each 100-foot thick layer .................................... 23

Figure 18. Map showing ordinary kriging variance (σ OK

2 ) in layers A, B, and C,

reflecting the relative uncertainty of kriging estimates based on the

availability and arrangement of borehole data .......................................... 24

Figure 19. Maps showing distribution of sediment-rich and sediment-poor areas in the

model domain .............................................................................. 25

Figure 20. Graph and map showing comparison of local cumulative frequency

distributions (CFDs) derived from multiple indicator kriging (mIK) at two

locations in Layer A, representing locations with low (location 1) and high

(location 2) estimation variance ......................................................... 26

Figure 21. Maps showing zones of hydraulic conductivity estimated from the kriged

sediment content in layer A ............................................................... 27

Figure 22. Map showing layer A’s bulk hydraulic conductivity (Kbulk) classified according

to the relative confidence in the kriged estimates as defined by the ordinary

kriging variance.............................................................................. 28

Tables

Table 1. Summary of composite stratigraphic units, thicknesses, and time spans as

defined by Anderson and Liszewski (1997) and used in this study..................... 4

Table 2. Results of nonparametric tests of similarity of median sediment content

between various groupings of younger and older composite units, where

sediment content is expressed in percent................................................ 15

Table 3. Results of nonparametric tests of similarity of median sediment content

between groups of younger and older composite units, where sediment

content is expressed as sediment thickness............................................. 15

Figures—Continued
vi

Conversion Factors, Datums, Abbreviations, and Acronyms

Conversion Factors

Multiply By To obtain

cubic foot (ft3

) 0.02832 cubic meter 

cubic mile (mi3

) 4.168 cubic kilometer 

foot (ft) 0.3048 meter

foot per day (ft/d) 0.3048 meter per day

foot squared per day (ft2

/d) 0.09290 meter squared per day 

gallon (gal) 0.003785 cubic meter

gallon per minute per foot [(gal/min)/ft)] 0.2070 liter per second per meter

mile (mi) 1.609 kilometer

square foot (ft2

) 0.09290 square meter

square mile (mi2

) 2.590 square kilometer 

Temperature in degrees Celsius (°C) may be converted to degrees Fahrenheit (°F) as follows:

°F=(1.8×°C)+32

Temperature in degrees Fahrenheit (°F) may be converted to degrees Celsius (°C) as follows:

°C=(°F-32)/1.8

Datums

Vertical coordinate information is referenced to the National Geodetic Vertical Datum of 1929 

(NGVD 29).

Horizontal coordinate information is referenced to the North American Datum of 1927 (NAD 27).

Altitude, as used in this report, refers to distance above the vertical datum.

*Transmissivity: The standard unit for transmissivity is cubic foot per day per square foot times 

foot of aquifer thickness [(ft3/d)/ft2]ft. In this report, the mathematically reduced form, foot 

squared per day (ft2/d), is used for convenience.

Abbreviations and Acronyms

Acronyms Definition

C.U. Composite stratigraphic unit

CFD Cumulative frequency distribution

ESRP Eastern Snake River Plain

GIS Geographic Information System

GPS Ground Positioning System

GSLIB Geostatistical Software Library

INL Idaho National Laboratory

K Hydraulic conductivity

Ka Thousand years

mIK Multiple indicator kriging

RWMC Radioactive Waste Management Complex

TAN Test Area North

USGS U.S. Geological Survey
Abstract

The spatial distribution of sediment in the eastern Snake 

River Plain aquifer was evaluated and modeled to improve 

the parameterization of hydraulic conductivity (K) for a 

subregional-scale ground-water flow model being developed 

by the U.S. Geological Survey. The aquifer is hosted within a 

layered series of permeable basalts within which intercalated 

beds of fine-grained sediment constitute local confining 

units. These sediments have K values as much as six orders 

of magnitude lower than the most permeable basalt, and 

previous flow-model calibrations have shown that hydraulic 

conductivity is sensitive to the proportion of intercalated 

sediment.

Stratigraphic data in the form of sediment thicknesses 

from 333 boreholes in and around the Idaho National 

Laboratory were evaluated as grouped subsets of lithologic 

units (composite units) corresponding to their relative timestratigraphic position. The results indicate that median 

sediment abundances of the stratigraphic units below the water 

table are statistically invariant (stationary) in a spatial sense 

and provide evidence of stationarity across geologic time, as 

well. Based on these results, the borehole data were kriged as 

two-dimensional spatial data sets representing the sediment 

content of the layers that discretize the ground-water flow 

model in the uppermost 300 feet of the aquifer.

Multiple indicator kriging (mIK) was used to model 

the geographic distribution of median sediment abundance 

within each layer by defining the local cumulative frequency 

distribution (CFD) of sediment via indicator variograms 

defined at multiple thresholds. The mIK approach is superior 

to ordinary kriging because it provides a statistically best 

estimate of sediment abundance (the local median) drawn 

from the distribution of local borehole data, independent of 

any assumption of normality. A methodology is proposed 

for delineating and constraining the assignment of hydraulic 

conductivity zones for parameter estimation, based on the 

locally estimated CFDs and relative kriging uncertainty. A 

kriging-based methodology improves the spatial resolution 

of hydraulic property zones that can be considered during 

parameter estimation and should improve calibration 

performance and sensitivity by more accurately reflecting the 

nuances of sediment distribution within the aquifer.

Introduction

Idaho’s largest fresh water resource, the Eastern Snake 

River Plain (ESRP) aquifer, occurs in the ESRP mafic volcanic 

province. This aquifer is composed primarily of pahoehoe 

basalt with intercalated sediments and minor occurrences of 

andesite and rhyolite. It is characterized by transmissivities 

and hydraulic conductivities ranging over six orders of 

magnitude and linear ground-water flow velocities up to 

10 ft/d (Ackerman, 1991; Anderson and others, 1999). 

A large number of well borings, logging records, and 

cores have been collected in the vicinity of the Idaho National 

Laboratory (INL) for defining the geohydrologic regime 

in this part of the ESRP aquifer (fig. 1). Improvements in 

characterization of the subsurface stratigraphy have led to 

increasingly sophisticated conceptual and digital models of 

the aquifer and of ground-water flow and transport at various 

spatial scales. A subregional-scale flow model (about 25 ×

75 mi; see fig. 1) currently is being tested and refined at the 

U.S. Geological Survey-INL (USGS-INL) project office. 

Calibration of the model has shown that assigned hydraulic 

conductivity values are sensitive to the relative abundance of 

sediment within the aquifer and that a more realistic model of 

the spatial variability of sediment will be required to refine the 

flow simulator. 

Geostatistical Modeling of Sediment Abundance in a 

Heterogeneous Basalt Aquifer at the Idaho National 

Laboratory, Idaho

By John A. Welhan1, Renee L. Farabaugh2, Melissa J. Merrick2, and Steven R. Anderson3

1

 Idaho Geological Survey, Department of Geosciences, Idaho State 

University, Pocatello, Idaho 83209.

2

 Department of Geosciences, Idaho State University, Pocatello, Idaho 

83209.

3 Formerly with U.S. Geological Survey, Idaho National Laboratory Project 

Office, P.O. Box 2230, Idaho Falls, Idaho 93403.
Figure 1. Study area, major topographic features, ground-water flow model extent, and locations of borehole data used for 

this analysis in and near the Idaho National Laboratory, Idaho.

 Geostatistical Modeling of Sediment Abundance in a Heterogeneous Basalt Aquifer, INL, Idaho

Ground-water flow model

boundary

Base from U.S. Geological Survey digital data, 1:24,000 and 1:100,000

Albers Equal-Area Conic projection Standard parallels 42°50'N and 

44°10'N, central meridian 113°00'W

Datum is North American Datum of 1927

RWMC

TAN

Idaho National Laboratory boundary

Borehole location

Selected facilities at the Idaho National Laboratory

 Radioactive waste management complex

 Test area north

EXPLANATION

IDAHO

BOISE

EASTERN

SNAKE RIVER

IDAHO PLAIN

NATIONAL

LABORATORY

Twin

Falls

Pocatello

Idaho

Falls

0

0 10

10

20 KILOMETERS

20 MILES

Mud

Lake

Big

Lost

River

Little

Lost

River

Birch

Creek

Mackay

Reservoir

Camas

Creek

SNAKE

RIVER

Big

Southern

Butte

East

Butte Middle

Butte

BITTERROOT

RANGE

LEMHI RANGE

LOST RIVER RANGE

PIONEER MOUNTAINS

EASTERN SNAKE RIVER PLAIN

WHITE KNOB

MOUNTAINS

26

93

20

26

33

22

28

20

93

113°30’ 113°00’ 112°30’

44°00’

30’

15’

45’

43°00’

Atomic

City

Howe

Arco

Mud

Lake

Terreton

TAN

RWMC
Purpose and Scope 

This study, conducted by the USGS in cooperation with 

the U.S. Department of Energy, was undertaken to develop 

a geostatistical model of sediment abundance in the aquifer 

that can be used to constrain future ground-water flow model 

calibrations. 

The spatial model of sediment abundance is based on 

borehole data compiled by Anderson and others (1996) in a 

region of the ESRP corresponding to the geographic extent 

of the ground-water model (fig. 1). The sediment model was 

created in a Geographic Information Systems (GIS)-compliant 

format that can interface with the ground-water flow model 

on a cell-by-cell basis, in specific depth-discretized layers 

corresponding to the digital flow model. The relative 

uncertainty of the sediment abundance estimates also had 

to be computed on a cell-by-cell basis because the degree 

of borehole coverage over the model domain is relatively 

sparse and highly variable. Geostatistical methods were 

employed to analyze the borehole data to select and justify 

an appropriate estimation method (in this case, multiple 

indicator kriging) and to quantify the relative uncertainty of 

the resulting estimates. Finally, the sediment model was used 

to demonstrate an approach by which sediment abundances 

could be used to develop a refined calibration process for the 

ground-water model. 

Previous Investigations

Hundreds of boreholes have been drilled in and near the 

INL during the past 50 years for purposes of characterizing 

the ESRP aquifer and monitoring chemical and radioactive 

wastes in the subsurface. The lithologies encountered in most 

of these boreholes have been characterized using geophysical 

logging techniques; less than 10 percent of all boreholes have 

been cored and selectively analyzed for properties such as 

geochronology, paleomagnetism, and bulk geochemistry in 

an ongoing attempt to characterize subsurface architecture. 

Systematic interpretation of these and other data sets has led 

to the development of a stratigraphic framework (for example, 

Anderson and others, 1996) and conceptual models of the 

aquifer (Lindholm and Vaccaro, 1988; Ackerman and others, 

2006), without which the quantitative analysis undertaken in 

this study would not be possible.

Interpreted Stratigraphy

Walker (1974) conducted one of the earliest stratigraphic 

studies of the INL’s subsurface using a combination of 

geologic and geophysical logs. His work was followed by 

another at the Radioactive Waste Management Complex 

(RWMC) in which sedimentary interbeds and individual 

basalt-flow groups (defined as packages of flows derived 

from a single shield volcano) were differentiated using 

natural-gamma logs (Barraclough and others, 1976). The first 

quantitative stratigraphic study using continuous cores was 

conducted at the RWMC using Potassium-Argon (K-Ar) age 

dating, paleomagnetic polarity and inclination measurements, 

and petrographic characterization (Kuntz and others, 1980). 

Similar studies in other areas have used 40Argon/39Argon 

and K-Ar ages (Champion and others, 1988; Lanphere and 

others, 1993, 1994; Champion and Lanphere, 1997; Champion 

and Herman, 2003). Detailed, small-scale stratigraphic 

interpretations based on basalt chemistry have been described 

by Reed and others (1997), Scarberry (2003), and Grimm 

Chadwick (2004) for selected coreholes. 

The work of Anderson and coworkers (Anderson and 

Lewis, 1989; Anderson, 1991; Anderson and Bowers, 1995; 

Anderson and Bartholomay, 1995) created the first and most 

carefully documented and vetted compilation of subsurface 

lithologic variability and stratigraphic correlations across 

the INL. Starting at selected INL facilities, their work with 

natural-gamma log records that had been calibrated to 

available core properties, led them to eventually describe and 

compile a subsurface stratigraphic database of 333 boreholes 

across the INL and vicinity (Anderson and others, 1996; 

Anderson and Liszewski, 1997). This database constitutes the 

basis of this study; the geographic locations of its boreholes 

are shown in figure 1.

Composited Stratigraphy

The database compiled by Anderson and others (1996) 

indicates that boreholes that penetrate the unsaturated zone of 

the aquifer have penetrated at least 178 identified basalt-flow 

groups, 103 sedimentary interbeds, 6 andesite-flow groups, 

and 4 rhyolite domes. Anderson and Liszewski (1997) first 

proposed a classification based on groups of units of similar 

age that was intended to show the principal stratigraphic 

features in the subsurface. They grouped individual basalt and 

sedimentary units into 14 composite stratigraphic units, each 

made up of 5 to 90 individual units. 

Table 1 summarizes these composite stratigraphic units. 

Composite unit 1, the youngest, is made up of 78 basalt-flow 

groups and 12 sedimentary interbeds. Composite unit 14, the 

oldest, is made up of only four identified basalt-flow groups 

and one sedimentary interbed. The decreasing number of 

individual units in successively older composite units may be 

due to larger and less-frequent volcanic eruptions as well as 

to the decreasing availability of borehole data with depth on 

which to base stratigraphic correlations. 

Introduction 
Composite units 1 through 7 represent time spans of 

approximately 65 to 248 thousand years each and, collectively, 

800 thousand years. Composite units 1, 2, and 3 generally 

comprise most of the unsaturated zone and only units 4 and 

older influence saturated ground-water flow in the aquifer. 

Composite units 8 through 14 span an age range of about 

800 thousand to 1.8 million years and are grouped together 

as one aggregated composite unit. As shown by Anderson 

and others (1996) and by subsequent work, these older units 

occur beneath most of the INL at considerable depths below 

the water table. In the northeastern part of the INL, however, 

where a stratigraphic discontinuity exposes the older units at 

the surface, they make up most or all of the unsaturated zone. 

The approximate locus of the stratigraphic discontinuity 

is shown in figure 2, and its subsurface character is shown 

in figure 3, together with the nomenclature used to identify 

the layers that discretize the aquifer in the subregional 

ground-water flow model. The water table in figure 3 reflects 

hydrologic conditions in 1996, at which time composite 

unit 5 was generally below the water table south of the 

discontinuity, whereas north of it, only composite units 12 

and older coincided with saturated conditions (Anderson and 

Liszewski, 1997). Despite recent drought conditions, a similar 

geohydrologic situation to that shown in figure 3 prevailed in 

2004.

Table 1. Summary of composite stratigraphic units, thicknesses, 

and time spans as defined by Anderson and Liszewski (1997) and 

used in this study.

[Abbreviations: ft, foot; Ka, thousands of years; n.a., not applicable]

Composite 

stratigraphic 

unit name

Number of 

wells fully 

penetrating 

the unit

Average 

thickness of 

the unit 

(ft)

Approximate 

time span of 

the unit 

(Ka)

1 300 62 5 – 250

2 170 108 250 – 350

3 140 163 350 – 440

4 143 137 440 – 515

5 112 126 515 – 580

6 68 107 580 – 650

7 13 283 650 – 800

8 - 14 n.a. n.a. 800 –1,800

Figure 2. Principal volcanic, structural, and geologic 

features in the study area, showing the inferred location 

of the Big Lost Trough in relation to major volcanic 

constructional zones. 

ID19-0017_fig02

A

A

A’

A’

Base from U.S. Geological Survey digital data, 1:24,000 and 1:100,000

Albers Equal-Area Conic projection Standard parallels 42°50'N and 

44°10'N, central meridian 113°00'W

Datum is North American Datum of 1927

EXPLANATION

Axial volcanic high and rift zones—From Geslin and others, 1999 

Approximate location of stratigraphic discontinuity—From Anderson 

and Liszewski, 1997 

Approximate extent of Big Lost Trough—From Geslin and others, 1999 

Trace of section—See figure 3

Idaho National Laboratory boundary

0

0 10

10

20 KILOMETERS

20 MILES

Mud

Lake

Big Lost

River

Little

Lost

River

Birch

Creek

Big

Southern

Butte

East

Butte

Middle

Butte

BITTERROOT

RANGE

LEMHI RANGE

LOST RIVER RANGE

EASTERN

SNAKE RIVER

PLAIN

113°15’ 113°00’ 112°30’

44°00’

43°30’

45’

45’

Atomic

City

Howe

Arco

Mud

Lake

Terreton

 Geostatistical Modeling of Sediment Abundance in a Heterogeneous Basalt Aquifer, INL, Idaho
Subsurface Geostatistical Studies

Inferring spatial continuity is one of the principal 

goals of geostatistics. Spatial autocorrelation is a statistical 

property that characterizes the statistical similarity of 

measurements separated by various distances. Various types 

of autocorrelation statistics known as variogram statistics 

can be used to describe the degree of similarity in space as 

well as to constrain geostatistical predictive models to make 

optimal estimates at unsampled locations (as well as quantify 

estimation uncertainty). Previous geostatistical work at the 

INL has demonstrated that subsurface lithologic properties 

in the ESRP aquifer are characterized by moderate to strong 

spatial autocorrelation, with greater statistical similarity over 

short distances that diminishes with distance. This fact can be 

exploited through various geostatistical modeling techniques 

(for example, kriging) to make best unbiased estimates 

of subsurface properties at unsampled locations. Previous 

geostatistical analysis of borehole data at the INL has been 

conducted with one of two general aims: (1) characterizing 

the nature of subsurface hydrologic variability (Welhan and 

Reed, 1997); and (2) modeling spatial variability to predict 

subsurface lithologic or hydrologic properties (Welhan and 

others, 2002a; Leecaster, 2002, 2004). 

In previous geostatistical modeling at the INL, different 

approaches have been used to achieve different goals. Welhan 

and others (2002a) used a stochastic indicator simulator 

to model the spatial configuration of sediments and basalt 

interflow zones in the aquifer beneath the Test Area North, 

Leecaster (2002, 2004) applied ordinary kriging to predict 

sedimentary interbed thicknesses and subsurface hydraulic 

properties of the aquifer in an area beneath the RWMC, and 

Gego and others (2002) used sequential Gaussian simulation 

to cast particle transport and breakthrough times in a stochastic 

context and demonstrate the application of probabilistic 

(Monte Carlo) ground-water modeling methods. 

Figure 3. Cross section A-A’ showing composite stratigraphic units 1 through 14, ground-water model layers A through 

F, and the relation between model layers and actual stratigraphy. Modified from Anderson and Liszewski (1997).

ID19-0017_fig03

9

10

9

8

7

6

4

3

2

1

5

14

13

12

10

9 8

11

?

?

?

4000

4400

3600

3200

4800

5000

4000

4400

3600

3200

4800

5000

A A'

VERTICAL SCALE GREATLY EXAGGERATED

FEET FEET

Vertical datum is National Geodetic Vertical Datum of 1929

Stratigraphic unit-age range per layer 

 with layer number—See table 1

 Less than 800,000 years

 800,000-1,800,000 years

Undifferentiated basalts and sediments

Ground-water flow model layer

Water-table level, 1996

Well

3 Layer D

Trace of section shown in figure 2

EXPLANATION Approximate location of stratigraphic discontinuity

Layer F

Layer A

Layer B

Layer C

Layer D

Layer E

Layer U

Introduction 
Geohydrologic Setting

The ESRP is an elongate, downwarped basin containing 

approximately 3,500 ft of Quaternary and late Tertiary basalt 

lava flows and sedimentary interbeds overlying late Tertiary 

rhyolitic volcanic rocks. The basin is characterized by plainsstyle volcanism (Greeley, 1982) comprising thin (3–100 ft), 

basaltic, pahoehoe lava flows and coalesced shield volcanoes. 

Estimates of basalt accumulation rates and eruptive recurrence 

intervals in 23 boreholes across the INL were published 

by Champion and others (2002). An analysis of their data 

reveals that accumulation rates are normally distributed, with 

pre-800 Ka (thousand years) basalts having a higher mean 

accumulation rate than post-800 Ka basalts (125 ft/Ka versus 

79 ft/Ka).4 Recurrence intervals are lognormally distributed 

with a ±1 σ interval about the median of 66 [33–167] Ka, and 

pre-800 Ka basalts have a shorter median recurrence interval 

than post-800 Ka basalts (54 Ka versus 100 Ka).5

Basalt lava flows are highly fractured at their margins, 

less so in their interiors, and the rubbly contacts between 

individual flows (interflow zones) can have greatly enhanced 

porosity and permeability that focuses ground-water 

movement (Lindholm and Vaccaro, 1988). The morphology 

and stacking of lava flows imparts a complex, anastomosed 

hydrostratigraphy of preferred flow paths at the local scale and 

a subhorizontal stratigraphic character at larger scales (Welhan 

and others, 2002b). 

Lengthy periods of quiescence during volcanic 

construction were marked by accumulation of fine-grained 

eolian, playa, and lacustrine sediment, and lesser amounts 

of coarser, channelized fluvial deposits within the basalt 

pile. Beneath the Idaho National Laboratory, sedimentary 

interbeds constitute about 15 percent by volume of the Snake 

River Plain aquifer and its unsaturated zone (Anderson and 

Liszewski, 1997). Because most interbeds are fine-grained, 

the presence of sediment in the subsurface tends to reduce the 

aquifer’s bulk hydraulic conductivity. The amount of sediment 

penetrated by individual boreholes ranges from less than 

5 percent of the stratigraphic column to more than 50 percent 

(up to 100 percent in shallow wells) near the terminus of the 

Big Lost and Little Lost Rivers in an inferred depositional 

center known as the Big Lost Trough (Anderson and others, 

1996; Geslin and others, 1999).

The Big Lost Trough is situated between the Axial 

Volcanic Zone to the east, the mountains and associated 

valleys of the Basin and Range to the north, and the northwesttrending volcanic rift zones in the northern and southern parts 

of the INL (fig. 2). The Axial Volcanic Zone is a topographic 

high of linearly clustered volcanic vents that have shed 

numerous basalt flows towards the northwest and southeast. 

Over time, these lava flows interfingered with, and trapped, 

sediments carried onto the Snake River Plain by the Big Lost 

and Little Lost Rivers. Similar clusters of volcanic vents 

associated with the northwest-trending rift zones may have 

further confined the geographic extent of this depositional 

center.

Although volcanic construction likely exerts the primary 

control on local geomorphology and the development of 

depositional centers, it is unclear whether and to what extent 

regional tectonic subsidence influences larger-scale patterns 

of sedimentary and volcanic accumulation. Paleomagnetic 

and radiometric evidence suggests that the overall rate of 

basalt accumulation in the vicinity of the INL decreased 

about 200 Ka before present (Anderson and others, 1997; 

Champion and others, 2002). This conjecture explains the 

increased proportion of sediment in the stratigraphic section 

over the latter half of the Quaternary, but leaves unanswered 

the question of whether sediment supply rates and (or) sources 

also changed.

Of particular significance to this study is the stratigraphic 

discontinuity between older and younger stratigraphic units 

in the northernmost part of the INL (figs. 2 and 3). South of 

the discontinuity, rocks and sediments younger than about 

800 Ka overlie older units, dipping northeastward into the 

discontinuity; north of the discontinuity, these younger units 

are almost completely absent. Because different patterns and 

frequency of sediment accumulation may characterize the 

older and younger composite units, the statistical nature of 

their sediment distribution was investigated to determine if 

units had to be segregated and modeled separately. 

Geostatistical Methods

Geostatistics is a specialized branch of statistical analysis 

concerned with the spatial relationships among data situated 

in two- or three-dimensional coordinate space. The process of 

geostatistical analysis and modeling includes: (1) assembling 

the data for analysis (data preparation); (2) calculating basic 

statistics, evaluating statistical stationarity, and performing 

data transformations to reduce the effect of outliers 

(exploratory analysis); (3) characterizing the spatial nature 

of the data (autocorrelation analysis); and (4) applying these 

statistical insights to create a model of spatial variability using 

a kriging estimator.

Data Selection 

The database of basalt and interbed stratigraphy compiled 

by Anderson and others (1996) was used as the basis for 

this work. However, several limitations in their data set had 

to be addressed prior to analysis. Foremost was the fact that 

4

 Significant at a 95 percent confidence level; a Student’s-t test of similarity 

of the mean accumulation rates returned a p-value of 0.017 (for pre- and post800 thousand year basalt sample sizes of 5 and 18, respectively.

5

 Significant at a 99.9 percent confidence level; a Mann-Whitney rank-sum 

test of similarity of the median recurrence intervals returned a p-value of 

0.0006 (for pre- and post-800 thousand year basalt sample sizes of 13 and 80, 

respectively).

 Geostatistical Modeling of Sediment Abundance in a Heterogeneous Basalt Aquifer, INL, Idaho
their stratigraphic interpretations represented best available 

knowledge at the time of compilation. Subsequent drilling, 

coring, and paleomagnetic and radiometric measurements have 

added to our understanding of subsurface stratigraphy beneath 

the INL. In many locations, new data have corroborated their 

interpretations; in other locations, new data have forced a 

re-evaluation of local stratigraphic relationships. 

Such limitations do not diminish the value of the 1996 

compilation for two reasons. First, the 1996 database defines 

the spatial occurrences of various lithologies, regardless of 

their interpreted stratigraphic context; therefore, the Anderson 

database represents a valid snapshot of subsurface sedimentary 

abundances at the time of data compilation. Second, as new 

subsurface data accrue, an appropriate temporal reference 

point will have to be selected on which future data syntheses 

can build. Although recent stratigraphic compilations were 

and currently are being made by various groups at the 

INL, these are not as extensively documented or vetted as 

the interpretations published by Anderson and coworkers. 

However, such compilations will provide a useful means of 

evaluating the accuracy and quality of earlier compilations 

after selecting an appropriate reference against which new 

information can be compared. 

For these reasons, the database of Anderson and 

others (1996) serves as a reasonable baseline for analyzing 

the subsurface distribution of sediment in the study area. 

Geostatistical analysis and modeling performed during the 

current work did not rely on the interpreted stratigraphy of 

Anderson and others (1996) except as it was used by Anderson 

and Liszewski (1997) to assign a composite unit stratigraphy. 

Data Preparation

The stratigraphic data set compiled by Anderson and 

others (1996) was in the form of a series of ASCII text files 

containing well locations, boring identifiers, lithologies, 

and interpreted stratigraphic unit designations. These were 

converted to spreadsheet form and then to .dbf tables, and 

ultimately incorporated in a relational database. 

Several factors influenced the decision regarding which 

data fields to incorporate in the database tables. Relational 

queries were used to extract information for specific analysis 

tasks, and the information was reformatted to be compatible 

with different graphical display and visualization software, 

as well as with Geographic Information System (GIS) 

and geostatistical software used in subsequent analyses. 

Accordingly, database fields were assigned, containing 

modified identifiers with which data can be located by 

lithologic and stratigraphic codes and, in the future, to 

trace the pedigree of interpreted stratigraphic assignments. 

Appropriate coordinates were added, reflecting the USGS 

subregional ground-water flow-model domain, and various 

modifications were carried out, including removing duplicate 

coordinates, assigning missing-value flags, and modifying 

anomalous stratigraphic designations which, in their original 

form, were incompatible with the software used for graphical 

display and visualization. These modifications are described 

below.

Modified Coordinates

Anderson and others (1996) reported their well locations 

in geographic coordinates of latitude and longitude, as decimal 

degrees. To be compatible with the USGS subregional flow 

model, Albers (NAD 27) coordinates for all wells were 

provided by L. Davis (U.S. Geological Survey, written 

commun., 2004) and used to update the Anderson database. 

Altitude coordinates, originally referenced to the National 

Geodetic Vertical Datum of 1929 (NGVD 29), were not 

updated from those provided by Anderson and others (1996).

Duplicate Coordinates

A total of 12 pairs, triplets, and quadruplets of well 

locations in the original Anderson data set shared identical 

geographic coordinates, which usually reflects wells that were 

drilled close together. These locations were updated with 

more accurate Global Positioning System (GPS) coordinates 

provided by L. Davis (U.S. Geological Survey, written 

commun., 2004). 

Ambiguous Stratigraphic Designations

The stratigraphic assignments made in the original 

database contain almost two dozen instances of a sedimentary 

unit interpreted to be a lens within a basalt unit; for 

example, the “B(1)s” sedimentary unit appears between two 

occurrences of the basalt unit “B(1)b”. Such a designation 

is physically impossible unless the basalt unit is actually 

an aggregate of two or more separate basalt flows. These 

designations had to be corrected to avoid incompatibility 

with the visualization software used. Omission of such lenses 

was considered inappropriate, as was assigning different 

stratigraphic designations to the basalt above and below the 

lens (which would affect the entire database). The solution 

was to aggregate the thickness of the sedimentary lens with its 

stratigraphically nearest interbed. This approach introduced 

the least disruption into the modified database while 

preserving the lens’s contribution to total sediment abundance. 

Database Organization

To facilitate organization and data-query capabilities, a 

number of minor modifications were made to the original data 

set published by Anderson and others (1996). Some borehole 

names were modified, usually by deleting spaces or other 

ambiguous alphanumeric designators, to ensure standardized 

Geostatistical Methods 
well names for linking tables. A field was added to designate 

the composite unit into which a particular stratigraphic unit 

falls (table 1). Finally, a field was incorporated to allow for 

future changes to the interpreted stratigraphic assignments to 

track the “pedigree” of future stratigraphic interpretations and 

database modifications.

Analysis Approach

Almost all of the lithologic information contained in the 

database compiled by Anderson and others (1996) reflects 

a dichotomous designation of sediment or basalt. Thus, the 

proportion of sediment encountered by a borehole can be 

expressed either as a cumulative thickness or as a percentage 

over a designated depth interval. Although both measures were 

evaluated during exploratory data analysis, the use of sediment 

percentages introduces statistical bias where sediment 

constitutes all of a composite unit’s thickness (most frequently, 

in locations where a composite unit is relatively thin or 

incompletely penetrated by boreholes) and skews the statistics 

of sediment abundance relative to basalt. Where practical, 

cumulative thicknesses as well as percentage of sediment were 

statistically analyzed; note that for model layers A, B, and C, 

each of which is 100-ft thick, measures of sediment thickness 

and percentage are numerically equivalent.

Exploratory analysis involved the computation of basic 

statistics for each composite unit’s sediment content, an 

evaluation of the effect of geographically clustered wells on 

the computed statistics, and an examination of the hypothesis 

of statistical stationarity (both in a geographic spatial sense 

and in a time-stratigraphic spatial sense). As one test of 

statistical stationarity, nonparametric hypothesis tests were 

used to evaluate the degree of similarity between sample 

medians. A Mann-Whitney test, equivalent to a Student’s t-test 

for normally distributed populations, was used to compare 

pairs of sample means; the Kruskal-Wallis test (equivalent 

to a one-way analysis of variance for normally distributed 

populations) was used to compare differences among three or 

more sample medians. All hypothesis tests were evaluated at 

a 95 percent confidence level (where a p-value greater than 

0.05 indicates no statistically significant difference among the 

samples).

Autocorrelation analysis focused on the computation of 

various spatial statistics to characterize the scale over which 

the thicknesses of composite units and interbeds persist. 

These statistical measures provided additional support for 

the hypothesis of stationarity and served as the basis for 

developing models of sediment thickness via kriging.

Kriging Approach

Kriging is a nearest-neighbor predictor that makes 

statistically “best” estimates at unsampled locations on the 

basis of available information in local neighborhoods. It 

differs from other linear prediction methods in its reliance 

on a model of spatial autocorrelation (inferred from the 

data) to assign weights to nearest data neighbors. Kriging 

estimates are optimal in the sense that their estimation error 

variances are minimized. In ordinary kriging, the error 

variance provides an estimate of prediction uncertainty 

based on the assumption that estimation errors are normally 

distributed. In this sense, ordinary kriging has been described 

as a “parametric” estimator (Isaaks and Srivastava, 1989). 

It is often unrealistic, however, to assume that estimation 

uncertainty is normally distributed or that the ordinary kriging 

estimate, the mean, is an appropriate estimator of central 

tendency at all kriged locations. To avoid this limitation, 

multiple indicator kriging (mIK) uses nonlinear (indicator) 

transforms of data around multiple thresholds to estimate 

the cumulative frequency distribution (CFD) within a local 

neighborhood. The local CFD is unconstrained by assumptions 

of normality, and provides a true measure of estimation 

uncertainty due to local data variability as well as a more 

appropriate estimate of central tendency, the median. The 

drawback is that autocorrelation statistics must be calculated 

and an autocorrelation model inferred at each threshold used 

to discretize the CFD, making mIK more time consuming to 

implement. However, the power, flexibility, and generality 

of mIK make it the superior kriging technique for modeling 

environmental data (Goovaerts, 1997). All kriging results were 

discretized and converted to industry-standard shapefile format 

on a quarter-mile grid corresponding to the cell centers of the 

USGS subregional flow model.

Software Selection

Data organization, query building, and data extraction 

were performed with Microsoft AccessTM database software; 

Microsoft ExcelTM was used for data manipulation, sorting, 

and indicator transformation. Basic statistical analyses were 

performed with MinitabTM and StatMostTM statistical software, 

and ArcMapTM Geographic Information System software was 

used for spatial data visualization and graphics generation. 

Variogram analysis was performed with the VARIOWINTM 

software package (Pannatier, 1996), and declustering analysis 

and multiple indicator kriging were performed with the 

Stanford Geostatistical Software LibraryTM (GSLIB) (Deutsch 

and Journel, 1998). Results of geostatistical modeling were 

imported as GSLIB grids into ArcMap and converted to 

shapefiles and raster grids for visualization and raster-based 

computations. 

 Geostatistical Modeling of Sediment Abundance in a Heterogeneous Basalt Aquifer, INL, Idaho
Geostatistical Analysis and Modeling

The goal of geostatistical analysis was to understand the 

distribution of sediment in a stratigraphic, three-dimensional 

spatial context in order to make defensible decisions during 

geostatistical modeling. In particular, statistical similarity 

among geologic units was evaluated to determine whether the 

sedimentary data were statistically stationary and could be 

modeled with kriging.

Exploratory Data Analysis 

The total thickness of sediment encountered in boreholes 

around the INL is summarized in figure 4. Most sediment-rich 

boreholes occur within or near the boundary of the putative 

sedimentary basin known as the Big Lost Trough (dashed 

polygon). Note the relative paucity of boreholes that penetrate 

composite units older than unit 6. The data have been 

segregated into groups representing wells that fully penetrate 

each of the composite stratigraphic units listed in table 1. 

More than 95 percent of all wells that encounter at least 

25 ft of sediment in each composite unit lie within or on the 

margins of the Big Lost Trough and are generally within 10 to 

15 mi of the present-day location of the Big Lost River.

General Considerations: Spatial Clustering and 

Partial Penetration

In statistically characterizing a spatial data set, the effect 

of spatial clustering must be considered. Spatially clustered 

data provide information redundancy in densely sampled areas 

and can bias calculated statistics toward values that reflect 

the most densely sampled areas, particularly if redundant 

sampling occurs in areas where the measured attribute is 

valued unusually high or low. In such cases, the spatial data 

must be appropriately weighted (declustered) to minimize this 

bias. In this study, the impact of such information redundancy 

is negligible. As shown in figure 5, the sample statistics 

changed only slightly when appropriate declustering weights 

were assigned, because boreholes are not preferentially 

clustered in areas of unusually high or low sediment content. 

Similar results were obtained when examining subsets of the 

data, broken down either by composite unit or by model layer.

The analysis of clustered data can be critical, however, 

in evaluating such characteristics as spatial autocorrelation 

and statistical stationarity of a population. Comparisons 

between areas of clustered and nonclustered data can provide 

evidence of spatial stationarity and help identify small-scale 

autocorrelation behavior. Furthermore, kriging automatically 

accounts for information redundancy, thus eliminating the 

need for declustering corrections. 

All subsequent analyses in this study were performed 

without declustering. In all cases in which raw and declustered 

data were compared, the effect of declustering was to diminish 

the apparent differences between the sample populations being 

compared; in this sense, the raw data provided a worst-case 

test for comparing statistical similarity among composite units.

Of far greater concern than clustered data was the effect 

of partial penetration. In a borehole which only partially 

penetrates a composite unit, the reported thickness of that 

unit or its sediment content is censored (quantified by an 

inequality, such as “greater than” an apparent thickness; Helsel 

and Hirsch, 2002). Bias arises when the censored information 

is compared with data from other boreholes that fully 

penetrate the same unit. Unlike other types of environmental 

information (where censoring occurs consistently in only 

one sense, such as when imposed by a detection limit, for 

example), borehole partial penetration bias is bi-directional 

and random. That is, if a composite unit happens to be 

sediment rich, then the proportion of sediment sampled by a 

partially penetrating borehole would represent an upper limit; 

alternatively, if the borehole penetrates a sediment-poor zone, 

the measurement would be biased low. In the case of the INL 

data, partial penetration should not introduce bias because 

such censoring is expected to be bi-directional and randomly 

expressed. 

Wherever possible, only those sediment thicknesses 

defined by fully penetrating boreholes were considered. 

Fortunately, sediment typically represents less than 15 percent 

of a borehole’s length, so that partial penetration bias is 

infrequent and has negligible effect. Except in the most 

infrequently sampled composite units, the attendant reduction 

in sample size is about 5 percent or less. It was necessary to 

use all available sedimentary thickness data because of limited 

sample size only in the aggregated composite unit 8 to 14, 

where just five deep boreholes fully penetrate. 

Only the ground-water model layers (A, B, and C 

[fig. 3]), corresponding to the uppermost 300 ft of the 

saturated zone, were considered in this study. An insufficient 

number of fully penetrating boreholes were available in layers 

D and deeper to justify a statistical analysis. For comparative 

purposes, the unsaturated zone (layer U) was also statistically 

analyzed but not kriged.

Sediment Abundance by Composite Unit

The statistical variability of sediment thickness in 

boreholes that fully penetrate each composite unit is shown in 

figure 6. The exception is the aggregate package, composite 

units 8 to 14, whose oldest units are fully penetrated by only 

five boreholes and which was analyzed without an adjustment 

for partial penetration. 

Geostatistical Analysis and Modeling 
ID19-0017_fig04a

Approximate extent of Big Lost Trough

Approximate location of stratigraphic discontinuity

Total thickness of sediment per borehole 

0 - 10 feet

10 - 25 feet

25 - 50 feet

50 - 100 feet

More than100 feet

EXPLANATION

113°00’ 112°30’

44°00’

43°30’

113°00’ 112°30’

44°00’

43°30’

113°00’ 112°30’

44°00’

43°30’

113°00’ 112°30’

44°00’

43°30’

0

0 10

10

20 KILOMETERS

20 MILES

Composite

stratigraphic

unit 1

5-250 thousand years

Composite

stratigraphic

unit 2

250-350 thousand years

Composite

stratigraphic

unit 3

350-440 thousand years

Composite

stratigraphic

unit 4

440-515 thousand years

Base from U.S. Geological Survey digital data, 1:24,000 and 1:100,000

Albers Equal-Area Conic projection Standard parallels 42°50'N and 

44°10'N, central meridian 113°00'W

Datum is North American Datum of 1927

Figure 4. Geographic distribution of sediment-rich and sediment-poor boreholes that fully penetrate composite stratigraphic units 

1 through 7 and 8 through 14, as defined by Anderson and Liszewski (1997). 

10 Geostatistical Modeling of Sediment Abundance in a Heterogeneous Basalt Aquifer, INL, Idaho
ID19-0017_fig04b

Approximate extent of Big Lost Trough

Approximate location of stratigraphic discontinuity

Total thickness of sediment per borehole 

0 - 10 feet

10 - 25 feet

25 - 50 feet

50 - 100 feet

More than100 feet

EXPLANATION

113°00’ 112°30’

44°00’

43°30’

113°00’ 112°30’

44°00’

43°30’

113°00’ 112°30’

44°00’

43°30’

113°00’ 112°30’

44°00’

43°30’

0

0 10

10

20 KILOMETERS

20 MILES

Composite

stratigraphic

unit 5

515-580 thousand years

Composite

stratigraphic

unit 6

580-650 thousand years

Composite

stratigraphic

unit 7

650-800 thousand years

Composite

stratigraphic

units 8-14

800-1800 thousand years

Base from U.S. Geological Survey digital data, 1:24,000 and 1:100,000

Albers Equal-Area Conic projection Standard parallels 42°50'N and 

44°10'N, central meridian 113°00'W

Datum is North American Datum of 1927

Figure 4.—Continued. 

Geostatistical Analysis and Modeling 11
ID19-0017_fig05

0 20 40 60 80 100

10

20

30

40

50

60

 0

10

20

30

40

50

60

 0

FREQUENCY

A.

B.

PERCENT SEDIMENT

Raw data

Mean 18.9

s 14.6

3rd quartile 22.7

2nd quartile 14.6

1st quartile 9.6 

Declustered data

Mean 16.4

s 16.9

3rd quartile 20.0

2nd quartile 11.3

1st quartile 5.7

Figure 5. Comparative distributions of sediment 

abundances showing the effect of spatial clustering on 

summary statistics. (A) sediment abundance data derived 

from spatially clustered boreholes may bias the statistics 

that summarize the central tendency and (or) spread of a 

distribution; (B) effect of redundant information is reduced 

by assigning appropriate declustering weights, leading to 

more representative summary statistics. Sample standard 

deviation is denoted by s.

The box plots in figure 6 clearly show that sediment 

comprises a significantly larger proportion of the thickness 

of composite unit 1 than of any of the older units. A similar 

pattern holds when considering sediment as a percentage of 

composite unit thickness (fig. 7). The distribution of sediment 

in composite unit 1 is strikingly different from that in all other 

composite units.6 Of the older composite units, some of whose 

sediment thickness distributions are statistically similar, the 

distribution in composite unit 5 is the most skewed and has the 

lowest median sediment content. However, it also represents 

the shortest accumulation period of any of the seven youngest 

composite units (65 Ka, table 1), raising the possibility that its 

anomalous character reflects an under-sampled stratigraphic 

interval.

To test this hypothesis and to further evaluate the case 

for statistical similarity over time, different groupings of 

composite units were combined into longer time spans of 

similar duration. Selected results for grouped units are shown 

in figure 8, for sediment percentages and figure 9, for sediment 

thicknesses. Nonparametric tests of similarity were conducted 

to identify whether median sediment abundances are more 

similar than not; statistical comparisons for various groupings 

that have similar distributions of sediment percentage and of 

sediment thickness, respectively, are summarized in tables 2

and 3. Overall, the results indicate that the relative sediment 

contents of composite units older than unit 1 are statistically 

more similar than they are different. An important caveat is 

that some of the sediment data representing composite units 8 

through 14 originates from partially penetrating boreholes and 

therefore may be biased. The inference of statistical similarity 

between these and younger units, therefore, remains tentative. 

These results underscore the difficulty in assessing 

statistical similarity in geologic units over time and space. 

Substantial variation in sediment thickness or in percentage 

of sediment in total thickness among units arises because of 

limited geographic representation or because the geologic time 

intervals being considered are not representative of long-term 

average sediment accumulation rates. The results of tables 2

and 3 indicate that time intervals of about 100–200 Ka or 

larger are required if meaningful comparisons are to be made 

in this sense. 

6

 Tests of similarity with other composite unit medians returned a MannWhitney p-value of 0.00 in all instances.

12  Geostatistical Modeling of Sediment Abundance in a Heterogeneous Basalt Aquifer, INL, Idaho
Figure 6. Distribution of thickness of total accumulated sediment in boreholes that fully penetrate individual composite 

units 1 through 7 and in boreholes that penetrate the oldest units 8 to 14 combined. Approximate time spans are from table 1.

ID19-0017_fig06

Max

90th

75th

50th

Min

10th

25th

ACCUMULATED SEDIMENT THICKNESS, FEET

n = 140

3

30

60

90

120

150

180

 0

2

n = 170

4

n = 143

5

n = 112

1

n = 300

COMPOSITE UNITS

0 200 400 600 800 1000 1200 1400 1600 1800

AGE, IN THOUSANDS OF YEARS

8-14

n = 81

n = Number of samples

7

n = 13

6

n = 68

Other than the relatively high sediment content noted in 

composite unit 1, sediment distributions appear similar when 

examined over sufficiently long time spans to average out 

short-term differences in sediment accumulation rates and 

sampling-induced noise. The data reveal no apparent trend 

over time in either accumulated sediment thickness or in its 

proportion relative to basalt. Composite unit 1’s high sediment 

content may reflect the period of relative volcanic quiescence 

over the past 200 Ka (Anderson and others, 1997; Champion 

and others, 2002), the interval during which most of composite 

unit 1 accumulated. If so, then sedimentation rates during this 

period may not have necessarily increased relative to earlier 

periods during which basalt accumulation rates were higher. 

Sediment Abundance by Model Layer

From the preceding analysis, it appears that except for 

the time span represented by composite unit 1, sediment 

abundance relative to basalt has not changed substantially over 

time. For the purposes of spatial modeling, such a conclusion 

is vital to any approach that aggregates borehole data into still 

larger groupings (such as ground-water model layers) that 

crosscut aquifer stratigraphy (fig. 3).

Locations of fully penetrating boreholes in each of the 

uppermost three ground-water model layers that represent 

the saturated zone are shown in figure 10. All subsequent 

statistical analyses were carried out on the basis of sediment 

percentages within these three, 100-ft thick layers (where 

percentages are numerically equivalent to thickness in feet). 

Geostatistical Analysis and Modeling 13
ID19-0017_fig08

Composite Unit 1 Composite Units 2-5

Composite Units 5-7

Composite Units 6-14

Composite Units 8-14

n = 300 n = 565

n = 193

n = 162

n = 81

Time span: 330 thousand years

Time span: about 1.2 million years

Time span: 250 thousand years

Time span: 285 thousand years

Time span: about 1 million years

n = Number of samples

FREQUENCY

 20

 40

 60

80

100

 0

 120

 140

 160

180

200

 20

 40

 60

80

100

 0

 120

 140

 160

180

200

 20

 40

 60

80

100

 0 0 20 40 60 80 100

SEDIMENT PERCENT

0 20 40 60 80 100

Figure 8. Distributions of sediment percentage in 

various groupings of composite units showing evidence 

of statistical similarity among the youngest (other than 

composite unit 1) and the oldest. 

Figure 7. Sediment abundance in each of composite units 

1 through 7 and in combined unit 8 to 14, as a percentage 

of composite unit thickness. 

ID19-0017_fig07 FREQUENCY

 20

 40

 60

80

100

 0

 20

 40

 60

80

100

 0

 20

 40

 60

80

100

 0

 20

 40

 60

80

100

 0

0 20 40 60 80 100

SEDIMENT PERCENT

0 20 40 60 80 100

Composite Unit 1 Composite Unit 2

Composite Unit 3 Composite Unit 4

Composite Unit 5 Composite Unit 6

Composite Unit 7 Composite Units 8-14

n = 300

156

n = 170

n =140 n = 143

n = 112 n = 68

n = 13 n = 81

Time span: 5-250 thousand years Time span: 250-350 thousand years

Time span: 350-440 thousand years Time span: 440-515 thousand years

Time span: 515-580 thousand years Time span: 580-650 thousand years

Time span: 650-800 thousand years Time span: 800-1,800 thousand years

n = Number of samples

14  Geostatistical Modeling of Sediment Abundance in a Heterogeneous Basalt Aquifer, INL, Idaho
7

 A Kruskal-Wallis nonparametric test of similarity among layers A, B, and 

C returned a p-value of 0.032 indicating that their sediment distributions are 

marginally different; layers A and C are the most different (Mann-Whitney 

p-value = 0.01), whereas layers A and B are statistically similar, as are B and 

C (Mann-Whitney p-values = 0.28 and 0.11, respectively).

Figure 9. Distributions of total sediment thickness among 

various groupings of composite units. For rocks older than 

unit 1, the similarity in histogram shape among various 

age-range groupings and the results of nonparametric 

tests of similarity (table 3) indicate that the amount of 

sediment relative to basalt has remained fairly similar over 

time. 

Table 2. Results of nonparametric tests of similarity of median 

sediment content between various groupings of younger and older 

composite units, where sediment content is expressed in percent.

[Abbreviations: Ka, thousands of years; M-W, Mann-Whitney test of 

similarity; p-value, probability value; vs., versus]

Groupings of 

composite units 

Time spans 

compared 

(Ka)

Test of similarity 

p-value

3–5 vs. 6– 14 230 vs. 2,020 0.504 (M-W)

2 –5 vs. 6– 14 330 vs. 2,020 .847 (M-W)

5 –7 vs. 8– 14 285 vs. 1,800 1.000 (M-W)

3 –7 vs. 8– 14 450 vs. 1,800 .300 (M-W)

ID19-0017_fig09

Time span: 250 thousand years Time span: 330 thousand years

Composite Unit 1 Composite Units 2-5

Composite Units 6-7 Composite Units 8-14

n = 300 n = 565

n = 81 n = 81

Time span: 220 thousand years Time span: about 1 million years

n = Number of samples

FREQUENCY

 20

 40

 60

80

100

 0

 120

 140

 160

180

200

 20

 40

 60

80

100

 0 0 20 40 60 80 120100 0 20 40 60 80 120100

SEDIMENT THICKNESS, IN FEET

Table 3. Results of nonparametric tests of similarity of median 

sediment content between groups of younger and older 

composite units, where sediment content is expressed as 

sediment thickness.

[Abbreviations: Ka, thousands of years; M-W, Mann-Whitney test of 

similarity; K-W, Kruskal-Wallis test of similarity; p-value, probability value; 

vs., versus]

Groupings of 

composite units 

Time spans 

compared 

(Ka)

Test of similarity 

p-value

2 –7 vs. 8–14 550 vs. 1,800 0.792 (M-W)

2 – 5 vs. 6 –7 vs. 8–14 330 vs. 220 vs. 1,800 .221 (K-W)

2 vs. 3 –7 vs. 8–14 100 vs. 350 vs. 1,800 .898 (K-W)

Similarities and differences among the sediment contents 

of these layers, as well as in the layer that represents the 

unsaturated zone, are summarized in figure 11. At the 

95-percent confidence level, nonparametric tests of similarity 

among pairs of model layers show that sediment distributions 

in layers A, B, and C range from marginally distinguishable 

to statistically indistinguishable.7

 Considering this together 

with the similarity of median abundances of sediment among 

composite units, the results indicate that the amount of 

sediment in the uppermost 300 ft of the aquifer is statistically 

stationary in the vertical sense. 

To evaluate statistical stationarity in the geographic sense, 

nonparametric tests of similarity were applied to sediment 

data from three different geographic groupings of boreholes 

in each of the model layers A, B, and C (fig. 10): (1) sediment 

data from nonclustered boreholes south of the stratigraphic 

discontinuity, (2) data from clustered boreholes south of the 

discontinuity, and (3) data from clustered boreholes north of 

the discontinuity. Histograms of these data groupings, together 

with the associated p-values for tests of similarity among 

the geographic groupings within each layer are presented in 

figure 12. Note that insufficient data precluded a comparison 

for the southern clustered area in layer C. Within each layer, 

the three data groupings are statistically indistinguishable 

at the 95-percent confidence level, indicating that median 

sediment abundance does not vary in a geographic sense. 

Geostatistical Analysis and Modeling 15
Figure 10. Availability of boreholes that fully penetrate each model layer and that were used to model the two-dimensional spatial 

abundance of sediment in each model layer. Boxes differentiate north- and south-clustered areas of wells used to evaluate spatial 

stationarity (see fig. 12).

ID19-0017_fig10

Clustered area of wells

Approximate extent of Big Lost Trough—

From Geslin and others, 1999 

Sediment content in borehole, as a 

percentage

0%

0.01-10%

10.01-20%

20.01-100%

EXPLANATION

113°00’ 112°30’

113°30’

44°00’

43°30’

113°00’ 112°30’

113°30’

44°00’

43°30’

113°00’ 112°30’

113°30’

44°00’

43°30’

LAYER A

LAYER C

LAYER B

0

0 10

10

20 KILOMETERS

20 MILES

Base from U.S. Geological Survey digital data, 1:24,000 and 1:100,000

Albers Equal-Area Conic projection Standard parallels 42°50'N and 

44°10'N, central meridian 113°00'W

Datum is North American Datum of 1927

Conclusions Derived from Exploratory Analysis

Sediment abundances within the composite units appear 

to be statistically similar over time when sufficiently long time 

spans (groupings of composite units) are compared. When the 

data are grouped by model layer, or by geographic area within 

model layers, median sediment abundances in the uppermost 

300 ft of the saturated zone appear to be statistically similar, 

further corroborating the inferences made on the basis of 

the composite unit analysis. Dissimilarity is apparent only 

in the very youngest composite unit (representing most of 

the shallowest part of the unsaturated zone) and makes the 

unsaturated zone, as a whole, statistically unique in terms of 

its sediment content. On this basis, it is concluded that the 

median sediment contents of rocks below the water table are 

statistically similar when grouped by either composite unit, 

by model layer, or by geographic area within a model layer, 

indicating that the sediment content of the aquifer may be 

statistically stationary in a three-dimensional spatial sense. 

16  Geostatistical Modeling of Sediment Abundance in a Heterogeneous Basalt Aquifer, INL, Idaho
ID19-0017_fig11 FREQUENCY

 20

 40

 60

80

100

 0

 20

 40

 60

80

100

 0

0 20 40 60 80 100

SEDIMENT PERCENT

0 20 40 60 80 100

Layer U (Vadose zone) Layer A

n = 222 n = 136

Layer C

n = 34

Layer B

n = 72

n = Number of samples

Figure 11. Distributions of sediment percentages in 

boreholes that fully penetrate layers A, B, and C, showing 

similar distributions, and in the unsaturated zone (layer U), 

which has a substantially different distribution because of 

the influence of composite unit 1 (see figures 6 and 7).

Autocorrelation Analysis 

The analysis of autocorrelation behavior involves the 

computation of sample statistics that characterize the degree 

of similarity of data values in space (information redundancy), 

followed by the inference of a model of autocorrelation 

structure from those sample statistics. In a nonrandom process 

such as sediment accumulation, the degree of similarity is 

greatest when pairs of closely spaced measurements are 

compared; measurements taken farther apart tend to be less 

similar. Autocorrelation statistics such as the autocovariance 

or the semivariogram provide a quantitative measure of 

the rate at which this similarity decreases with separation 

distance. These statistical measures, collectively referred to 

as variogram statistics, are the foundation on which all twopoint geostatistical inference and modeling is based. They 

are represented by plots of a variance-like statistical quantity 

(for example, the semivariogram) as a function of separation 

distance. The geometric structure of such plots is described 

in terms of a model of spatial autocorrelation, a curvilinear 

function characterized by a nugget (uncorrelated variance 

at zero separation), a range (separation distance over which 

autocorrelation disappears), and a sill (uncorrelated sample 

variance among pairs of data values separated by distances 

greater than the range). The presence of a horizontal sill that 

corresponds to the sample variance is the strongest indication 

that the data are second-order stationary, as required for most 

kriging applications.

Variography of Composite Units

The results of autocorrelation analysis of sediment 

thickness within individual composite units are summarized in 

figure 13, for composite units 1 through 6 (the spatial density 

of boreholes within composite units 7 and 8-14 was too low to 

warrant a similar analysis). A relative variogram estimator (the 

inverted correlogram) was used to simplify intercomparison 

of autocorrelation structure among the composite units. 

Variogram statistics were calculated at two different spatial 

lags to better assess the structure of the variogram’s transition 

zone and the magnitude of its nugget. Because of its small 

sample size (n = 68), the variogram for composite unit 6 is 

poorly defined, but it appears to fit the same autocorrelation 

model as other composite units. Similar autocorrelation 

structure is evident for all composite units shown, with a welldefined sill and clear transition zone, a relative nugget of about 

40 percent, and similar autocorrelation ranges. 

It is noteworthy that all composite units share a similar 

autocorrelation structural model. Their well-defined sills, 

together with the uniformity of their autocorrelation ranges 

(from 20,000 to 25,000 ft) corroborates the inference of 

statistical stationarity across time-stratigraphic boundaries. In 

particular, the similarity of composite unit 1’s autocorrelation 

structure to other composite units substantiates the conjecture 

that its relative sediment content differs not because of 

fundamental changes in sedimentation rate or spatial patterns, 

but because of a decrease in lava accumulation rate.

The statistical similarity of medians among different 

geographic groupings of boreholes within each model layer 

(fig. 12) and the similarity in autocorrelation structure among 

composite units (fig. 13) indicate that sediment abundance 

is a geographically stationary population. Although this 

conclusion must remain tentative in view of the paucity of 

spatial information at greater borehole depths and in certain 

parts of the study area, a preponderance of evidence supports 

the hypothesis of statistical stationarity within at least the 

uppermost 300 ft of the ground-water flow model domain.

Geostatistical Analysis and Modeling 17
ID19-0017_fig12

FREQUENCY

 10

 20

 30

40

50

 0

 10

 20

 30

40

50

 0

10

20

30

40

50

 0

0 20 40 60 80 100

SEDIMENT PERCENT

0 20 40 60 80 100

0 20 40 60 80 100

South Nonclustered Area North Clustered Area South Clustered Area

n = 65 n = 29 n = 42

n = 34 n = 16 n = 22

n = 27 n = 7

Mann-Whitney p-value = 0.975

Kruskal-Wallis p-value = 0.090

Mann-Whitney p-value = 0.509

Model layer A

Model layer B

Model layer C

n = Number of samples

Figure 12. Distributions of sediment percentage from fully penetrating wells in model layers A, B, 

and C in different areas north and south of the stratigraphic discontinuity and among clustered and 

nonclustered areas shown in figure 10. Probability values (p-values) refer to nonparametric tests of 

similarity among geographic groupings of wells.

18  Geostatistical Modeling of Sediment Abundance in a Heterogeneous Basalt Aquifer, INL, Idaho
ID19-0017_fig13 INVERTED CORRELOGRAM STATISTIC, ρ(h)

0.2

0.4

0.6

0.8

1.0

 0

1.4

1.6

1.2

0.2

0.4

0.6

0.8

1.0

 0

1.4

1.6

1.2

0.2

0.4

0.6

0.8

1.0

 0

1.4

1.6

1.2

0 20 40 60 80 100 120 0 20 40 60 80 100 120

PAIR SEPARATION DISTANCE, IN THOUSANDS OF FEET

Composite Unit 1 Composite Unit 2

Composite Unit 3 Composite Unit 4

Composite Unit 5 Composite Unit 6

s2 = 1153

s2 = 326

s2 = 184

s2 = 215

s2 = 918

s2 = 243

Range

Nugget

Exponential variogram with range=

 20,000-25,000 feet

Exponential variogram with range=

 15,000-30,000 feet and sills varying 

 by ± 5%

Sample variance, s2

Variogram statistics calculated with

 Small lag size

 Large lag size 

Figure 13. Relative variograms of sediment thickness within composite units 1 through 6, based on the inverted 

correlogram statistic. Variograms are plotted at two lag sizes to elucidate the variogram range and nugget. All fitted models 

have a relative nugget of 40 percent of the sample variance. 

Geostatistical Analysis and Modeling 19
Indicator Variography of Model Layers

For the purposes of mIK, the cumulative frequency 

distribution (CFD) of sediment percentage is discretized and 

approximated by a series of indicator transformations of the 

original data. By way of example, the borehole data for layer 

A (expressed as sediment percentages) was transformed 

around twelve indicator thresholds at 1 percent increments 

from 1 to 12 percent sediment (few boreholes in this layer 

have higher sediment contents, making transformation 

at higher thresholds impractical). At each threshold, the 

following indicator transformation was applied: for example, 

at a 7 percent sediment threshold, each data value was 

transformed using the following relationship:

I

i

 = 1 if Zi < 7 percent

I

i

 = 0 if Zi > 7 percent,

where

Zi are the original data values (sediment percentages), and the

I

i

are their indicator-transformed values.

The average of the transformed values provides an 

estimate of the global CFD’s p-quantile (probability value) at 

that threshold:

p = (1/n)Σ I

i

 . (1)

The thresholds at which the indicator transformations 

were carried out for each layer and the resulting indicatorderived probabilities are shown in figure 14. The dashed 

lines portray the global CFDs that would be plotted using 

the untransformed data. In effect, multiple indicator 

transformation allows a distribution’s shape to be estimated 

from a series of nonlinear transformations of the original data. 

A different set of thresholds was applied to each model layer 

because the shapes of their global CFDs and the amount of 

data in the upper tails of their distributions differed greatly.

Multiple indicator kriging uses the autocorrelation 

structure of the indicator variograms at each threshold to 

calculate optimal kriging weights, from which the local CFD’s 

conditional probabilities are estimated at each of several 

thresholds. The principal strength of mIK is that different 

spatial autocorrelation structure can be accommodated across 

the range of the data being kriged. Figure 15 shows examples 

of layer A’s indicator variograms and how the autocorrelation 

ID19-0017_fig14

PROPORTION OF SEDIMENT IN LAYER, AS PERCENT

0 20 40 60 80 100

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

PROBABILITY

Layer B

Layer C

Layer A

Global cumulative frequency distribution

Indicator-derived probability

Threshold used for indicator transformation

Figure 14. Cumulative frequency distributions (CFDs) 

of sediment percentage in boreholes that fully penetrate 

layers A, B, and C showing CFD probabilities calculated 

from equation 1. 

structure changes from one threshold to the next, reflecting 

the degree to which high and low values are spatially 

autocorrelated. All indicator variograms display clearly 

defined sills, providing further support for the hypothesis 

that the sediment abundance data reflect a globally stationary 

variable (Isaaks and Srivastava, 1989). 

20 Geostatistical Modeling of Sediment Abundance in a Heterogeneous Basalt Aquifer, INL, Idaho
ID19-0017_fig15

γ (h)

 0 0

 0 5 10 15 0 5 10 15

 0

0.3

0.4

0.2

0.1

 0

0.3

0.4

0.2

0.1

 0

0.3

0.4

0.2

0.1

50 100 150 50 100 150

PAIR SEPARATION DISTANCE, IN THOUSANDS OF FEET

9 percent sediment threshold

1 percent sediment threshold 3 percent sediment threshold

5 percent sediment threshold 7 percent sediment threshold

11 percent sediment threshold

Exponential variogram

Sample variance, s2

Calculated variogram statistics 

Figure 15. Examples of multiple indicator semivariograms at selected thresholds for layer A, showing the nature of the 

autocorrelation structure and the progressive change in nuggets, sills, and ranges as thresholds increase. Note that the 

horizontal axes for the 1- and 3-percent thresholds are 10 times smaller than for larger thresholds.

Geostatistical Analysis and Modeling  21
Indicator variogram model parameters for the 12 

thresholds subsequently used to make mIK estimates for 

layer A are summarized in figure 16. Because the choice 

of CFD thresholds is arbitrary, care must be taken to 

ensure that the variogram model parameters inferred at one 

threshold are compatible with parameter values inferred 

at neighboring thresholds (Goovaerts, 1997). A relatively 

smooth variation of nugget, sill, and range parameters, such 

as shown in figure 16, ensures that order-relation corrections 

in the resulting kriging estimates are minimized and that the 

mIK-estimated local CFDs will accurately reflect local data 

dispersion (Deutsch and Journel, 1998). Note that the indicator 

variogram ranges inferred for layer A (5,000 to 50,000 ft) 

bracket the 20,000-25,000 value deduced for composite units 

4 and 5 which comprise the bulk of layer A, as expected for 

a second-order stationary variable (Goovaerts, 1997). The 

autocorrelation ranges tend to be higher at higher thresholds, 

that is, where thicker accumulations of sediment extend over 

larger geographic areas. These nuggets, sills and ranges were 

subsequently used in multiple indicator kriging to define local 

CFDs at each kriging location.

Kriging Approach

Kriging relies on data within a nearest-neighbor search 

area to make estimates at unsampled locations on the basis 

of the autocorrelation structure of the data. Ideally, where 

sampling density is high, locally-defined variograms may be 

calculated solely on the basis of the neighboring data (Whelan 

and others, 2001), but in most environmental and subsurface 

situations, data density is almost never sufficient to permit 

such an approach and a global variogram must be inferred 

from the entire data set. Evidence of statistical stationarity 

is then used to justify the relevance of the global variogram 

to all local neighborhoods where kriging estimates are to be 

made. In this study, evidence for stationarity was based on the 

similarity of medians within different geographic portions of 

each model layer as well as on the autocorrelation structure 

of their indicator variograms and of the composite units’ 

variograms. 

Multiple indicator kriging was performed with the 

GSLIBTM software package (Deutsch and Journel, 1998). 

Indicator variogram parameters were defined at sediment 

thresholds of 1, 3, 5, 7, 9, and 11 percent for all layers, as 

well as at thresholds of 18 and 25 percent for layer B and 20, 

40, and 60 percent for layer C. Linear tail approximations 

were used to approximate the CFD above the highest 

indicator threshold considered (fig. 14). Even though the 

maximum range of the indicator variograms was about 

60,000 ft (fig. 16), a kriging search radius of 150,000 ft was 

used to return estimates in even the most data-poor regions 

of each model layer. The disadvantage of this approach is 

that unsubstantiated estimates can arise in data-poor areas 

(reflecting simple averages of distant data values), but it is 

the only way to generate estimates over the entire domain. If 

preferred, such artifacts could be filtered and removed using a 

relative variance criterion (such as that generated by ordinary 

kriging; Isaaks and Srivastava, 1989) or assigned a reduced 

confidence level.

0.10

0.20

0.15

0.05

0.00

0 2 4 6 8 10 12 14

15

0

60

45

30

CUMULATIVE FREQUENCY DISTRIBUTION THRESHOLD, PERCENT SEDIMENT

SEMIVARIOGRAM NUGGET, SILL

SEMIVARIOGRAM RANGE, IN THOUSANDS OF FEET

Layer A

Nugget

Partial sill

Range

Figure 16. Results of multiple indicator variogram analysis for layer A, showing the progressive change in 

autocorrelation structure across a range of indicator thresholds. 

22  Geostatistical Modeling of Sediment Abundance in a Heterogeneous Basalt Aquifer, INL, Idaho
Kriging Results

Multiple Indicator Kriging

The median sediment contents of each layer (as read from 

the local CFDs estimated by mIK) are shown in figure 17, 

in relation to the surface extent of the Big Lost Trough 

depositional basin (dashed polygon; after Geslin and others, 

1999). Estimation variances also can be determined from 

the local CFDs, but they do not directly reflect the relative 

uncertainty that arises from the proximity of neighboring data 

Sediment content in borehole, as 

a percentage of layer thickness 

0%

0.01-10%

10.01-20%

20.01-100%

EXPLANATION

50%

0

Median

sediment

content

80%

0

Median

sediment

content

60%

0

Median

sediment

content

113°00’ 112°30’

113°30’

44°00’

43°30’

113°00’ 112°30’

113°30’

44°00’

43°30’

113°00’ 112°30’

113°30’

44°00’

43°30’

LAYER A

LAYER C

LAYER B

0

0 10

10

20 KILOMETERS

20 MILES

Extent of

Big Lost Trough

Extent of

Big Lost Trough

Extent of

Big Lost Trough

INL boundary

Big

Lost River

INL boundary

Big

Lost River

INL boundary

Big

Lost River

Active model 

boundary

for this layer

Active model 

boundary

for this layer

Active model 

boundary

for this layer

Base from U.S. Geological Survey digital data, 1:24,000 and 1:100,000

Albers Equal-Area Conic projection Standard parallels 42°50'N and 

44°10'N, central meridian 113°00'W

Datum is North American Datum of 1927

values. Because the latter is a far more important criterion 

when kriged sediment content is used to constrain hydraulic 

conductivity estimates, kriging uncertainty was quantified 

using the ordinary kriging variance rather than the local 

CFD variance returned by mIK. A global variogram of the 

untransformed data in each layer (with a nugget of 0 and 

normalized to a sill of 1.0) was used to generate maps of 

ordinary kriging variance, σ OK

2 (for example, figure 18). The 

ordinary kriging variance provides a quantitative basis for 

classifying the modeled sediment abundance according to the 

uncertainty in data-poor areas relative to data-rich areas.

Figure 17. Median percentage of sediment derived from multiple indicator kriging of layers A, B, and C. Percentages are 

equivalent to thickness, in feet, within each 100-foot thick layer. 

Geostatistical Analysis and Modeling  23
Geohydrologic Implications

Figure 17 shows the highest sediment contents in layers 

A and B generally occur within the boundary of the Big Lost 

Trough as delineated by Geslin and others (1999). For layer 

C, substantial thicknesses of sediment are detected outside 

the boundary, notably to the east of the depositional basin, 

indicating that the Big Lost Trough’s spatial extent may 

have changed over time, possibly in response to volcanic 

construction. 

The predicted distribution of sediment within individual 

layers bears little resemblance to the spatial extents of 

sediment-rich and sediment-poor areas currently implemented 

in the USGS subregional flow model. A map (fig. 19) of the 

ID19-0017_fig18

Ordinary kriging variance

(σ2

OKrelative)

Sediment content in borehole, as 

a percentage of layer thickness 

0%

0.01-10%

10.01-20%

20.01-100%

EXPLANATION

1

0

113°00’ 112°30’

113°30’

44°00’

43°30’

113°00’ 112°30’

113°30’

44°00’

43°30’

113°00’ 112°30’

113°30’

44°00’

43°30’

LAYER A

LAYER C

LAYER B

0

0 10

10

20 KILOMETERS

20 MILES

INL boundary

Big

Lost River

INL boundary

Big

Lost River

INL boundary

Big

Lost River

Active model 

boundary

for this layer

Active model 

boundary

for this layer

Active model 

boundary

for this layer

Base from U.S. Geological Survey digital data, 1:24,000 and 1:100,000

Albers Equal-Area Conic projection Standard parallels 42°50'N and 

44°10'N, central meridian 113°00'W

Datum is North American Datum of 1927

Figure 18. Ordinary kriging variance (σ OK

2 ) in layers A, B, and C, reflecting the relative uncertainty of kriging estimates based 

on the availability and arrangement of borehole data.

subregional flow model’s sediment-rich and sediment-poor 

areas, (Ackerman and others, 2006), shows the boundary 

between sediment-rich and sediment-poor areas inferred 

from borehole data collected from the uppermost 300 ft of 

the aquifer, regardless of depth. Their classification is based 

on an 11-percent sediment classification threshold (fig. 5A). 

Figure 19B shows the combined mIK results for the model 

layers recalculated as the median sediment content of the 

uppermost 300 ft of the aquifer. The mIK model provides a 

refined picture of the spatial distribution of sediment in the 

upper part of the aquifer, by representing it on a layer-by-layer 

basis and reflecting greater confidence in areas with the most 

borehole control. 

24  Geostatistical Modeling of Sediment Abundance in a Heterogeneous Basalt Aquifer, INL, Idaho
Application to Scaling of Aquifer Hydraulic 

Conductivity 

The USGS modeling team has shown that optimized 

hydraulic conductivity values assigned during calibration 

of the subregional ground-water flow model tend to scale 

proportionately with zones of high- and low-sediment 

abundance shown in figure 19A (D.J. Ackerman, U.S. 

Geological Survey, oral commun., 2005). That is, the 

calibration process tends to assign lower values of hydraulic 

conductivity where sediment comprises more than 11 percent 

of the aquifer’s uppermost 300 ft. Similar relations between 

aquifer hydraulic properties and sediment abundance have 

been noted for calibrated vertical:horizontal anisotropy 

ratios and specific yields (D.J. Ackerman, U.S. Geological 

Survey, oral commun., 2005). However, the coarse sediment 

zonation currently implemented in the subregional flow model 

introduces undesirable modeling artifacts, such as streamline 

refraction across a zonal boundary where hydraulic properties 

change abruptly (G.W. Rattray, U.S. Geological Survey, 

written commun., 2005).

ID19-0017_fig19

Sediment content of the 

uppermost 300 feet of the 

aquifer

More than 11%

Less than 11%

EXPLANATION

EXPLANATION

50%

2

Median sediment content 

of the uppermost 300 

feet of the aquifer

Borehole

113°00’ 112°30’

113°30’

44°00’

43°30’

113°00’ 112°30’

113°30’

44°00’

43°30’

A. B.

Sediment-rich

Sediment-poor

0

0 10

10

20 KILOMETERS

20 MILES

Extent of

Big Lost Trough

INL boundary

Big

Lost River

INL boundary

Big

Lost River

Active model 

boundary

for this layer

Active model 

boundary

for this layer

Boundary

between sediment-rich 

and sediment-poor areas

Base from U.S. Geological Survey digital data, 1:24,000 and 1:100,000

Albers Equal-Area Conic projection Standard parallels 42°50'N and 

44°10'N, central meridian 113°00'W

Datum is North American Datum of 1927

Figure 19. Distribution of sediment-rich and sediment-poor areas in the model domain. (A) as currently implemented in the U.S. 

Geological Survey subregional ground-water flow model, based on an 11-percent sediment threshold, regardless of depth; and 

(B) median sediment content of the uppermost 300 feet of the aquifer, as determined by multiple-indicator kriging. The extent of 

the Big Lost Trough (dashed polygon) and the boundary between sediment-rich and sediment-poor areas in (A) are shown for 

comparison.

The enhanced spatial resolution of the mIK models 

allows us to consider more sophisticated approaches to the 

scaling of hydraulic properties that depend on sediment 

content. Two types of methodologies are proposed: (1) those 

based solely on a threshold sediment content, such as the 

median, and (2) those based on a combination of a threshold 

and the uncertainty associated with the threshold estimate. 

In what follows, we consider hydraulic conductivity as an 

example of how scaling can be implemented, although similar 

lines of reasoning could be developed for porosity, storativity 

or any property whose values differ sufficiently between pure 

sediment and pure basalt.

Scaling involves expressing hydraulic conductivity 

as a function of the aquifer’s sediment content, such as its 

median value (or some other, operationally defined level 

like the 11-percent threshold currently implemented in the 

USGS model). In a horizontally layered aquifer composed of 

fractured, high-permeability basalt and fine-grained, lowpermeability sediment, the vertical and horizontal components 

of hydraulic conductivity (or, more precisely, the components 

Geostatistical Analysis and Modeling  25
of hydraulic conductivity parallel and normal to layering) can 

be expressed as a function of the relative thicknesses of the 

two lithologies (Freeze and Cherry, 1979):

	 KH = [Z50 • KS

 + (100–Z50) • KB]/100 (2)

and

	 KV = 100/[Z50/KS + (100–Z50) • KB]/100 , (3)

where

KS is the hydraulic conductivity of sediment;

KB is the hydraulic conductivity of basalt;

Z50 is the median (50th percentile) sediment content; 

and

KH and

KV

are the resulting components of hydraulic 

conductivity.

A threshold other than the median could also serve as 

the scaling variable. For example, at the 11-percent threshold 

currently implemented in the USGS subregional flow model, 

a relative bulk hydraulic conductivity can be estimated 

from the probability, P11, of not exceeding the 11-percent 

sediment threshold (the threshold currently implemented in the 

subregional model): 

Kbulk = [(1–P11) • KS + P11 • KB] . (4)

Examples of two mIK-derived local CFDs in a sedimentpoor and a sediment-rich location of layer A are shown 

in figure 20, having median sediment contents of Z50-1 = 

4 percent and Z50-2 = 36 percent, respectively. For values 

of KS = 0.2 ft/d and KS = 11,000 ft/d (Garabedian, 1989), 

equation 2 results in KH-1 = 10,500 ft/d and KH-2 = 7,000 ft/d; 

corresponding KV values derived from equation 3 are 5.0 

and 0.6 ft/d. The resulting anisotropy ratios (KH/KY) range 

from 2,100:1 to 13,000:1, similar to those obtained in the 

flow model calibration process (4,000:1 to 15,000:1; D.J. 

Ackerman, oral commun., 2005).

Using equation 4, with values of P11-1 = 0.92 and P11-2 = 

0.32, and the same KS and KB values, the calculated values of 

Kbulk at sites 1 and 2 are 10,100 and 3,500 ft/d, respectively. 

Figure 21 illustrates how such a sediment-based scaling 

approach could be used to delineate hydraulic conductivity 

zones for zone-based parameter estimation. Although only 

three zones are considered in this illustration, any number 

could be defined. The mIK median sediment contents were 

classified into three categories (0–5-, 5–20-, 20–50-percent 

sediment) and hydraulic conductivity values representing the 

midpoints of those ranges were calculated from equations 2 

and 3. The resulting values are displayed as color-classified 

maps of KH and KH/KV in figure 21A and 21B, respectively. 

Similarly, local CFD p-quantiles representing the mid-points 

ID19-0017_fig20

Big

Lost River

113°00’ 112°30’

113°30’

44°00’

43°30’

LAYER A

0

0 10

10

20 KILOMETERS

20 MILES

Z50-1 Z50-2

P11-1

P11-2

0 40 80

PERCENT SEDIMENT

 0

0.2

0.4

0.6

0.8

 1

CUMULATIVE FREQUENCY DISTRIBUTION

20 60 100

Location 1

Location 2

Base from U.S. Geological Survey digital data, 1:24,000 and 1:100,000

Albers Equal-Area Conic projection Standard parallels 42°50'N and 

44°10'N, central meridian 113°00'W

Datum is North American Datum of 1927

Median sediment content

Sediment content in borehole, 

as a percentage of layer 

thickness

0%

0.01-10%

10.01-20%

20.01-100%

EXPLANATION

50%

0%

Active model 

boundary

for this layer

INL boundary

Location 2

Location 1

Figure 20. Comparison of local cumulative frequency 

distributions (CFDs) derived from multiple indicator kriging (mIK) 

at two locations in Layer A, representing locations with low 

(location 1) and high (location 2) estimation variance. Median 

sediment contents (Z50) computed from the local CFDs are used 

to estimate horizontal and vertical hydraulic conductivity with 

equations 2 and 3. Local CFD probabilities for an 11 percent 

sediment threshold were used to estimate bulk hydraulic 

conductivities with equation 4. 

26  Geostatistical Modeling of Sediment Abundance in a Heterogeneous Basalt Aquifer, INL, Idaho
ID19-0017_fig21

113°00’ 112°30’

113°30’

44°00’

43°30’

113°00’ 112°30’

113°30’

44°00’

43°30’

113°00’ 112°30’

113°30’

44°00’

43°30’

A. Horizontal 

hydraulic

conductivity

C. Bulk 

hydraulic

conductivity

B. Horizontal hydraulic 

conductivity/

Vertical hydraulic 

conductivity ratio

0

0 10

10

20 KILOMETERS

20 MILES

INL boundary

Big

Lost River

INL boundary

Big

Lost River

INL boundary

Big

Lost River

Active model 

boundary

for this layer

Active model 

boundary

for this layer

Active model 

boundary

for this layer

Base from U.S. Geological Survey digital data, 1:24,000 and 1:100,000

Albers Equal-Area Conic projection Standard parallels 42°50'N and 

44°10'N, central meridian 113°00'W

Datum is North American Datum of 1927

Kbulk, in feet per day

9,900

6,600

3,300

EXPLANATION

KH, in feet per day

10,700

9,600

7,200

EXPLANATION

KH/Kv, ratio

12,500

6,000

1,300

EXPLANATION

of the P11 ranges (0.2-0.4), (0.4-0.8), and (0.8-1.0) were 

used in equation 4 to calculate values of Kbulk (fig. 21C). 

The resulting zone maps are applied in the same way that 

figure 19A is used in the parameter estimation process.

During the parameter optimization process, it may be 

advantageous to allow K values to vary more widely in areas 

where the zone boundaries are least certain. To achieve this, 

the zones defined in figure 21 can be further subdivided 

according to the kriging variance (fig. 18) into zones of low 

and high uncertainty. An example is shown in figure 22, 

where layer A’s domain was classified around an ordinary 

kriging variance threshold of 0.6 on a relative scale of 0 to 1. 

The hydraulic conductivity zones shown in figure 21C were 

subsequently reclassified into five new categories representing 

areas of greater and lesser confidence in the scaled K values. 

Such an approach would allow K assignments during the 

calibration process to range more widely in zones of lower 

confidence. 

Figure 21. Zones of hydraulic conductivity estimated from the kriged sediment content in layer A. (A) horizontal hydraulic 

conductivity (KH) calculated as a weighted arithmetic mean of sediment and basalt end-members, using equation 2 with Kbasalt

= 11,000 and Ksediment = 0.2 feet per day; (B) the anisotropy ratio derived from KH and the vertical hydraulic conductivity (KV

) 

calculated using equation 3; and (C) bulk hydraulic conductivity (Kbulk) calculated as a weighted arithmetic mean, using equation 4 

and the probability of not exceeding an 11-percent sediment threshold.

Geostatistical Analysis and Modeling  27
Summary and Conclusions

The relative sediment content of the eastern Snake 

River Plain aquifer is a key geologic variable that affects the 

assignment of hydraulic conductivity values during calibration 

of ground-water flow models. Geostatistical analysis was used 

to estimate sediment content in the uppermost 300 ft of the 

U.S. Geological Survey (USGS) subregional ground-water 

flow model in order to improve future model calibrations.

The statistical characteristics of the sediment content of 

the aquifer, as measured in 333 boreholes located across the 

Idaho National Laboratory, were analyzed using conventional 

(nonspatial) and spatial statistical methods. Borehole 

data were segregated into various groups by stratigraphic 

age, depth, and geographic location in order to determine 

the degree of statistical stationarity that characterizes the 

model domain. The results indicate that the stratigraphic 

distribution of sediment within the aquifer, when evaluated 

over sufficiently large intervals of time (200 thousand years 

or greater), appears to be statistically stationary in a vertical 

sense as well as in a geographic sense. This finding makes it 

possible to model the spatial distribution of sediment in the 

uppermost three layers of the USGS subregional ground-water 

flow model using two-dimensional kriging within individual 

layers. 

ID19-0017_fig22

113°00’ 112°30’

113°30’

44°00’

43°30’

INL boundary

Big Lost River

Active model 

boundary

for this layer

LAYER A

0

0 10

10

20 KILOMETERS

20 MILES

Kbulk, in feet per day

9,900

6,600

3,300

}

}

Sediment content in borehole, 

as a percentage 

0%

0.01-10%

10.01-20%

20.01-100%

EXPLANATION

Base from U.S. Geological Survey digital data, 1:24,000 and 1:100,000

Albers Equal-Area Conic projection Standard parallels 42°50'N and 

44°10'N, central meridian 113°00'W

Datum is North American Datum of 1927

Figure 22. Layer A’s bulk hydraulic conductivity (Kbulk) classified according to the relative confidence in the kriged estimates 

as defined by the ordinary kriging variance (fig. 18, Layer A). Within each hydraulic conductivity zone, lighter colors indicate 

less uncertainty (better borehole control in data-rich areas) and darker colors indicate greater uncertainty in data-poor areas.

Multiple indicator variogram analysis and multiple 

indicator kriging (mIK) were used to model the spatial 

distribution of sediment. The use of mIK has three principal 

advantages over ordinary kriging: (1) local uncertainty of 

the kriging estimate is quantified more accurately than is 

possible in ordinary kriging; (2) differing degrees of spatial 

autocorrelation can be used to describe sediment-rich and 

sediment-poor areas; and (3) local frequency distributions 

that characterize a layer’s sediment contents can be used to 

delineate hydraulic-property zones that depend on sediment 

content. 

The kriging approach honors all available borehole 

sediment data, provides better spatial resolution for delineating 

zones of relative sediment abundance, and provides a basis 

for making parameter assignments that reflect the uncertainty 

of the kriged sediment estimates. By delineating hydraulic 

conductivity zones with mIK and considering kriging 

uncertainty in the model calibration process, the ultimate 

benefit of such an approach will be to more tightly constrain 

parameter assignments in data-rich areas. This approach 

could be incorporated in future development of the USGS 

subregional ground-water flow model and development of 

methods for integrating the uncertainty associated with zonal 

boundaries into the parameter estimation process. 

28  Geostatistical Modeling of Sediment Abundance in a Heterogeneous Basalt Aquifer, INL, Idaho
Acknowledgments

We thank Joseph Rousseau of the U.S. Geological Survey 

(USGS) for his encouragement and support of this work, for 

his helpful comments as an informal reviewer, and for the 

opportunity to publish these results. We also thank Linda 

Davis (USGS) for assistance with the borehole database and 

GIS-related work, and Daniel Ackerman and Gordon Rattray 

of the USGS for sharing insights and data on the subregional 

flow model. We are especially grateful to Molly Leecaster 

(CH2M-Washington Group), John Nimmo (USGS), and Kip 

Bossong (USGS) for their thoughtful reviews that materially 

improved the manuscript and to Mary Hodges (USGS) for her 

comments as an informal reviewer. Finally, the senior author 

wishes to recognize co-author Steven Anderson for his work 

on elucidating the stratigraphic framework of the Eastern 

Snake River Plain, without which this study would not have 

been possible.

References Cited

Ackerman, D.J., 1991, Transmissivity of the Snake River 

Plain aquifer at the Idaho National Engineering Laboratory, 

Idaho: U.S. Geological Survey Water-Resources 

Investigations Report 91-4058 (DOE/ID-22097), 35 p.

Ackerman, D.J., Rattray, G.W., Rousseau, J.P., Davis, L.C., 

and Orr, B.R., 2006, A conceptual model of ground-water 

flow in the eastern Snake River Plain aquifer at the Idaho 

National Laboratory and vicinity with implications for 

contaminant transport: U.S. Geological Survey Scientific 

Investigations Report 2006-5122 (DOE-ID-22198), 62 p.

Anderson, S.R., 1991, Stratigraphy of the unsaturated zone 

and uppermost part of the Snake River Plain aquifer at 

the Idaho Chemical Processing Plant and Test Reactors 

Area, Idaho National Engineering Laboratory, Idaho: U.S. 

Geological Survey Water-Resources Investigations Report 

91-4010 (DOE/ID-22095), 71 p.

Anderson, S.R., Ackerman, D.J., Liszewski, M.J., and 

Freiburger, R.M., 1996, Stratigraphic data for wells at and 

near the Idaho National Engineering Laboratory, Idaho: 

U.S. Geological Survey Open-File Report 96-248 (DOE/ID22127), 27 p. and 1 diskette.

Anderson, S.R., and Bartholomay, R.C., 1995, Use of naturalgamma logs and cores for determining stratigraphic 

relations of basalt and sediment at the Radioactive Waste 

Management Complex, Idaho National Engineering 

Laboratory, Idaho: Journal of the Idaho Academy of 

Science, v. 31, no. 1, p. 1-10.

Anderson, S.R., and Bowers, Beverly, 1995, Stratigraphy 

of the unsaturated zone and uppermost part of the Snake 

River Plain aquifer at Test Area North, Idaho National 

Engineering Laboratory, Idaho: U.S. Geological Survey 

Water-Resources Investigations Report 95-4130 (DOE/ID22122), 47 p.

Anderson, S.R., Kuntz, M.A., and Davis, L.C., 1999, Geologic 

controls of hydraulic conductivity in the Snake River Plain 

aquifer at and near the Idaho National Engineering and 

Environmental Laboratory, Idaho: U.S. Geological Survey 

Water-Resources Investigations Report 99-4033 (DOE/ID22155), 38 p.

Anderson, S.R., and Lewis, B.D., 1989, Stratigraphy of the 

unsaturated zone at the Radioactive Waste Management 

Complex, Idaho National Engineering Laboratory, Idaho: 

U.S. Geological Survey Water-Resources Investigations 

Report 89-4065 (DOE/ID-22080), 54 p.

Anderson, S.R., and Liszewski, M.J., 1997, Stratigraphy of 

the unsaturated zone and the Snake River Plain aquifer at 

and near the Idaho National Engineering Laboratory, Idaho: 

U.S. Geological Survey Water-Resources Investigations 

Report 97-4183 (DOE/ID-22142), 65 p.

Anderson, S.R., Liszewski, M.J., and Cecil, L.D., 1997, 

Geologic ages and accumulation rates of basalt-flow groups 

and sedimentary interbeds in selected wells at the Idaho 

National Engineering Laboratory, Idaho: U.S. Geological 

Survey Water-Resources Investigations Report 97-4010 

(DOE/ID-22134), 39 p.

Barraclough, J.T., Robertson, J.B., Janzer, V.J., and Saindon, 

L.G., 1976, Hydrology of the solid waste burial ground, as 

related to the potential migration of radionuclides, Idaho 

National Engineering Laboratory: U.S. Geological Survey 

Open-File Report 76-471, (IDO 22056), 183 p.

Champion D.E., and Herman, T.C., 2003, Paleomagnetism 

of basaltic lava flows in coreholes ICPP-213, ICPP-214, 

ICPP-215, and USGS 128 near the Vadose Research Park, 

Idaho Nuclear Technology and Engineering Center, Idaho 

National Engineering and Environmental Laboratory, Idaho: 

U.S. Geological Survey Open-File Report 03-483, 15 p.

Champion, D.E., and Lanphere, M.A., 1997, Age and 

paleomagnetism of basaltic lava flows in corehole ANLOBS-AQ-014 at Argonne National Laboratory-West, Idaho 

National Engineering and Environmental Laboratory: U.S. 

Geological Survey Open-File Report 97-700, 34 p.

Champion, D.E., Lanphere, M.A., and Kuntz, M.A., 1988, 

Evidence for a new geomagnetic reversal from lava flows in 

Idaho: Discussion of short polarity reversals in the Brunhes 

and late Matuyama polarity chrons: Journal of Geophysical 

Research, v. 93, p. 11667-11680.

References Cited  29
Champion, D.E., Lanphere, M.A., Anderson, S.R., and 

Kuntz, M.A., 2002, Accumulation and subsidence of late 

Pleistocene basaltic lava flows of the eastern Snake River 

Plain, Idaho, in Link, P.K., and Mink, L.L., eds., Geology, 

Hydrogeology, and Environmental Remediation: Idaho 

National Engineering and Environmental Laboratory, 

Eastern Snake River Plain, Idaho: Boulder, Colo., 

Geological Society of America, Special Paper 353, p. 175-

192.

Deutsch, C.V., and Journel, A.G., 1998, GSLIB: Geostatistical 

Software Library and User’s Guide (2d ed.): New York, 

Oxford University Press, 340 p. 

Freeze, R.A., and Cherry, J.A., 1979, Groundwater: 

Englewood Cliffs, N.J., Prentice-Hall, 604 p.

Garabedian, S.P., 1989, Hydrology and digital simulation 

of the regional aquifer system, eastern Snake River Plain, 

Idaho: U.S. Geological Survey Open-File Report 87-237, 

151 p.

Gego, E.L., Johnson, G.S., Hankin, M.R., and Welhan, 

J.A., 2002, Modeling ground water flow and transport 

in the Snake River Plain aquifer: a stochastic approach, 

in Link, P.K., Mink, R., and Ralston, D., eds., Geology, 

hydrogeology and environmental remediation, Idaho 

National Engineering and Environmental Laboratory, 

eastern Snake River Plain, Idaho: Geological Society of 

America, Special Paper 353, p. 249-261.

Geslin, J.K., Link, P.K., and Fanning, C.M., 1999, Highprecision provenance determination using detrital-zircon 

ages and petrography of Quaternary sands on the eastern 

Snake River Plain, Idaho: Geology, v. 27, no. 4, p. 295-298. 

Goovaerts, P., 1997, Geostatistics for natural resources 

evaluation: New York, Oxford University Press, 483 p. 

Greeley, Ronald, 1982, The Snake River Plain, Idaho: 

Representative of a new category of volcanism: Journal of 

Geophysical Research, v. 87, p. 2705-2712.

Grimm Chadwick, Claire, 2004, Petrogenesis of an evolved 

olivine tholeiite and chemical stratigraphy of cores 

USGS 127, 128, and 129, Idaho National Engineering 

and Environmental Laboratory: Pocatello, Idaho State 

University, M.S. thesis, 100 p., 4 app.

Helsel, D.R., and Hirsch, R.M., 2002, Statistical methods 

in water resources: Techniques of Water-Resources 

Investigations of the United States Geological Survey, book 

4, chap. A3, 510 p. 

Isaaks, E.H., and Srivastava, R.M., 1989, An introduction to 

applied geostatistics: New York, Oxford University Press, 

561 p.

Kuntz, M.A., Dalrymple, G.B., Champion, D.E., and Doherty, 

D.J., 1980, Petrography and paleomagnetism of volcanic 

rocks at the Radioactive Waste Management Complex, 

Idaho National Engineering Laboratory, Idaho, with an 

evaluation of potential volcanic hazards: U.S. Geological 

Survey Open-File Report 80-388, 63 p.

Lanphere, M.A., Champion, D.E., and Kuntz, M.A., 1993, 

Petrography, age, and paleomagnetism of basalt lava flows 

in coreholes well 80, NRF 89-04, NRF 89-05, and ICPP 

123, Idaho National Engineering Laboratory, Idaho: U.S. 

Geological Survey Open-File Report 93-327, 40 p.

Lanphere, M.A., Kuntz, M.A., and Champion, D.E., 1994, 

Petrography, age, and paleomagnetism of basalt lava flows 

in coreholes at Test Area North (TAN), Idaho National 

Engineering Laboratory, Idaho: U.S. Geological Survey 

Open-File Report 94-686, 49 p.

30 Geostatistical Modeling of Sediment Abundance in a Heterogeneous Basalt Aquifer, INL, Idaho
Leecaster, M.K., 2002, Geostatistic modeling of subsurface 

characteristics in the Radioactive Waste Management 

Complex region, Operable Unit 7-13/14: Idaho Falls, Idaho, 

U.S. Department of Energy, INEEL/EXT-02-00029, 26 p.

Leecaster, M.K., 2004, Fiscal year 2004 geostatistical 

modeling of lithologic characteristics in the Radioactive 

Waste Management Complex for OU 7-13/14: Idaho Falls, 

Idaho, U.S. Department of Energy, ICP/EXT-04-00494, 

27 p.

Lindholm, G.F., and Vaccaro, J.J., 1988, Region 2, Columbia 

Lava Plateau, in Hydrogeology, Back, W., Rosenshein, 

J.S., and Seabar, P.R., eds.: Geology of North America, v. 

O-2, Decade of North American Geology, Boulder, Colo., 

Geological Society of America, p. 37-50.

Pannatier, Y., 1996, VARIOWIN: Software for Spatial Data 

Analysis in 2D: New York, Springer-Verlag, 91 p.

Reed, M.F., Bartholomay, R.C., and Hughes, S.S., 1997, 

Geochemistry and stratigraphic correlation of basalt lavas 

beneath the Idaho Chemical Processing Plant, Idaho 

National Engineering Laboratory: Environmental Geology, 

v. 30, p. 108-118.

Scarberry, K.C., 2003, Volcanology, geochemistry, and 

stratigraphy of the F basalt flow group, eastern Snake River 

Plain, Idaho: Pocatello, Idaho State University, M.S. thesis, 

139 p.

Walker, G.P.L., 1974, Compound and simple lava flows and 

flood basalts. Bulletin of Volcanology 35, p. 579-590.

Welhan, J.A., and Reed, M.F., 1997, Geostatistical analysis 

of regional hydraulic conductivity variations in the Snake 

River Plain aquifer, eastern Idaho: Geological Society of 

America, Bulletin 109, p. 855-868.

Welhan, J.A., Clemo, T. and Gego, E., 2002a, Stochastic 

simulation of aquifer heterogeneity in a layered basalt 

aquifer system, eastern Snake River Plain, Idaho, in Link, 

P.K., Mink, R. and Ralston, D., eds., Geology, hydrogeology 

and environmental remediation, Idaho National Engineering 

and Environmental Laboratory, eastern Snake River Plain, 

Idaho: Geological Society of America Special Paper 353, 

p. 225-247.

Welhan, J.A., Johannesen, C.M., Reeves, K., Clemo, T., 

Glover, J.A., and Bosworth, K., 2002b, Morphology of 

inflated pahoehoe lavas and spatial architecture of their 

porous and permeable zones, eastern Snake River Plain, 

Idaho, in Link, P.K., Mink, R. and Ralston, D., eds., 

Geology, hydrogeology and environmental remediation, 

Idaho National Engineering and Environmental Laboratory, 

eastern Snake River Plain, Idaho: Geological Society of 

America Special Paper 353, p. 135-150

Whelan, B.M., McBratney, A.B., and Minasny, B., 2001, 

VESPER: Spatial prediction software for precision 

agriculture, in Grenier, G. and Blackmore, S., eds., 

Proceedings of the European Conference on Precision 

Agriculture: Montpelier, France, Agro-Montpelier, 

p. 139-144.

References Cited  31
This page is intentionally left blank.

32  Geostatistical Modeling of Sediment Abundance in a Heterogeneous Basalt Aquifer, INL, Idaho
Manuscript approved for publication, December 9, 2006

Prepared by the USGS Publishing Network, 

Publishing Service Center, Tacoma, Washington

Bill Gibbs

Bob Crist

Debra Grillo

Sharon Wahlstrom

Bobbie Jo Richey

For more information concerning the research in this report, contact the 

Director, Idaho Water Science Center 

230 Collins Road

Boise, Idaho 83702

http://id.water.usgs.gov
SIR 2006-5316 Welhan and others

Geostatistical Modeling of Sediment Abundance in a Heterogeneous 

Basalt Aquifer at the Idaho National Laboratory, Idaho