U.S. Department of the Interior

U.S. Geological Survey

Scientific Investigations Report 2014–5127

Version 1.1, January 2015

Groundwater Resources Program

Numerical Simulation of Groundwater Flow in the 

Columbia Plateau Regional Aquifer System, Idaho, 

Oregon, and Washington
Cover: Lower left: Mean annual groundwater pumpage for current conditions (2000–2007).

Center foreground: Model-generated upper Grande Ronde unit composite drawdowns 

from predevelopment to current conditions (2000–2007).

Lower right: Mean annual groundwater recharge for current conditions (2000–2007).

Background: Photograph of Frenchman’s Coulee in the vicinity of Quincy, Washington. 

(Photograph taken by Peter Prehn, private citizen, 2008. Used with permission.)
Numerical Simulation of Groundwater Flow in the 

Columbia Plateau Regional Aquifer System, Idaho, 

Oregon, and Washington

By D. Matthew Ely, Erick R. Burns, David S. Morgan, and John J. Vaccaro

Groundwater Resources Program

Scientific Investigations Report 2014–5127

Version 1.1, January 2015 

U.S. Department of the Interior

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

SALLY JEWELL, Secretary

U.S. Geological Survey

Suzette M. Kimball, Acting Director 

U.S. Geological Survey, Reston, Virginia

First release: 2014

Revised January 2015 (ver. 1.1)

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

resources, natural hazards, and the environment, visit http://www.usgs.gov or call 1–888–ASK–USGS

For an overview of USGS information products, including maps, imagery, and publications, 

visit http://www.usgs.gov/pubprod

To order this and other USGS information products, visit http://store.usgs.gov

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

U.S. Government.

Although this information product, for the most part, is in the public domain, it also may contain copyrighted materials 

as noted in the text. Permission to reproduce copyrighted items must be secured from the copyright owner.

Suggested citation:

Ely, D.M., Burns, E.R., Morgan, D.S., and Vaccaro, J.J., 2014, Numerical simulation of groundwater flow in the 

Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington (ver. 1.1, January 2015): U.S. Geological 

Survey Scientific Investigations Report 2014–5127, 90 p., http://dx.doi.org/10.3133/sir20145127.

ISSN 2328-0328 (online)
iii

Contents

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

Introduction.....................................................................................................................................................2

Purpose and Scope ..............................................................................................................................2

Previous Investigations........................................................................................................................2

Description of Study Area ............................................................................................................................4

Location and Setting ............................................................................................................................4

History of Water Resource Development.........................................................................................4

Hydrogeologic Setting .........................................................................................................................6

Hydrologic Budget Components ........................................................................................................6

Description of Groundwater-Flow System ..............................................................................................10

Hydrogeologic Units...........................................................................................................................10

Overburden Unit.........................................................................................................................11

Saddle Mountains Unit .............................................................................................................11

Mabton Interbed Unit................................................................................................................11

Wanapum Unit............................................................................................................................11

Vantage Interbed Unit...............................................................................................................11

Grande Ronde Unit ....................................................................................................................11

Older Bedrock Unit ....................................................................................................................11

Hydraulic Characteristics of Hydrogeologic Units........................................................................12

Groundwater Flow and Occurrence in Basalt Units.....................................................................12

Numerical Simulation of Groundwater Flow...........................................................................................15

Modeling Strategy ..............................................................................................................................15

Upstream Weighting Package.................................................................................................15

Newton Solver............................................................................................................................15

Layering and Spatial Discretization.................................................................................................16

Temporal Discretization.....................................................................................................................23

Hydraulic Properties...........................................................................................................................23

Horizontal Hydraulic Conductivity...........................................................................................23

Vertical Hydraulic Conductivity ...............................................................................................24

Storage Properties ....................................................................................................................24

Geologic Structures ..................................................................................................................24

Hydrologic Boundaries ......................................................................................................................26

Hydraulic Conditions Along the Periphery of the Model.....................................................26

Groundwater Pumpage.............................................................................................................26

Groundwater Recharge ............................................................................................................30

Streams and Surface-Water Features ...................................................................................33

Model Calibration and Sensitivity ....................................................................................................36

Observations Used in Model Calibration ..............................................................................36

Water-Level Altitudes, Water-Level Altitude Changes, and Associated Errors.....36

Streamflow Observations and Associated Errors.......................................................38
iv

Contents—Continued

Numerical Simulation of Groundwater Flow—Continued 

Model Calibration and Sensitivity—Continued

Model Evaluation .......................................................................................................................40

Parameterization and Regularization ............................................................................40

Optimal Parameter Estimates .........................................................................................45

Statistical Measures of Model Fit..................................................................................50

Steady-State Model Fit and Model Error......................................................................50

Comparison of Simulated and Measured Hydraulic Heads .............................51

Comparison of Simulated and Measured Streamflow ......................................51

Transient Model Fit and Model Error.............................................................................56

Comparison of Simulated and Measured Hydraulic Heads .............................56

Comparison of Simulated and Measured Streamflow ......................................62

Simulated Potentiometric Surfaces......................................................................69

Model Uncertainty and Limitations.........................................................................................69

Model Application...............................................................................................................................78

Groundwater Budget.................................................................................................................78

Simulated Effects of Pumping, Commingling Wells, and Irrigation Recharge.................80

Major Findings From Numerical Simulation of Groundwater Flow .....................................................84

Summary........................................................................................................................................................85

References Cited..........................................................................................................................................86

Figures

1. Map showing Columbia Plateau Regional Aquifer System study area and structural

regions, Idaho, Oregon, and Washington.................................................................................3

2. Map showing location of major river basins and study subareas, Oregon and

Washington....................................................................................................................................5

3. Map showing surficial distribution of generalized hydrogeologic units and

locations of hydrogeologic sections, Columbia Plateau Regional Aquifer System,

Idaho, Oregon, and Washington.................................................................................................7

4. Generalized hydrogeologic sections of the Columbia Plateau Regional Aquifer

System, Idaho, Oregon, and Washington .................................................................................8

5. Diagram of features within a typical Columbia River Basalt Group flow..........................13

6. Map showing estimated groundwater levels and directions of lateral groundwater

movement for the Grande Ronde unit, Columbia Plateau Regional Aquifer System,

Idaho, Oregon, and Washington...............................................................................................14

7. Map showing location and extent of the groundwater-flow model of the Columbia

Plateau Regional Aquifer System, Idaho, Oregon, and Washington .................................17

8. Maps showing extent and thickness of the Overburden unit, combined Saddle

Mountains basalt and Mabton interbed unit,combined Wanapum basalt and

Vantage interbed unit, and Grande Ronde basalt unit, Columbia Plateau Regional

Aquifer System, Idaho, Oregon, and Washington.................................................................19
v

9. Map showing location and hydraulic characteristic group of the groundwater

model structure (horizontal-flow barriers), Columbia Plateau Regional Aquifer

System, Idaho, Oregon, and Washington ...............................................................................25

10. Map showing location of the groundwater-flow model river and drain cells, and the

no-flow boundary, Columbia Plateau Regional Aquifer System, Idaho, Oregon, and

Washington..................................................................................................................................27

11. Graph showing estimated groundwater pumpage, by category, Columbia Plateau

Regional Aquifer System, Idaho, Oregon, and Washington, 1985–2007............................28

12. Graph showing estimates of groundwater pumpage from Cline and Collins (1992)

and SOil WATer (SOWAT) balance (Kahle and others (2011), and simulated

groundwater-pumpage, Columbia Plateau Regional Aquifer System, Idaho, Oregon,

and Washington..........................................................................................................................29

13. Map showing mean annual groundwater pumpage for current conditions

(2000–2007), Columbia Plateau Regional Aquifer System, Idaho, Oregon, and

Washington..................................................................................................................................31

14. Maps showing mean annual groundwater recharge for predevelopment and

current conditions (2000–2007), Columbia Plateau Regional Aquifer System, Idaho,

Oregon, and Washington...........................................................................................................34

15. Graph showing distance of the well bottom from the estimated top of Grande

Ronde Basalt compared with water levels measured in wells for the Palouse Slope

eastern well groups, Washington............................................................................................37

16. Graph showing hydraulic-head observation weights used in the steady-state

model, Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington...39

17. Map showing locations of base-flow observation sites used to calibrate the

steady-state and transient models of the Columbia Plateau Regional Aquifer

System, Idaho, Oregon, and Washington ...............................................................................41

18. Map showing initial model parameter zones and geologic structure for the

Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington................44

19. Maps showing simplified final model parameter zones for the Overburden unit and

Columbia River Basalt Group, Columbia Plateau Regional Aquifer System, Idaho,

Oregon, and Washington...........................................................................................................46

20. Map showing unweighted hydraulic-head residuals (simulated minus measured)

from the steady-state model, Columbia Plateau Regional Aquifer System, Idaho,

Oregon, and Washington...........................................................................................................52

21. Graph showing comparison of unweighted and weighted hydraulic-head residuals

with weight, Columbia Plateau Regional Aquifer System, Idaho, Oregon, and

Washington..................................................................................................................................53

22. Map showing comparison of the unweighted residuals with simulated current

conditions (2000–2007) average drawdown for the Saddle Mountains unit,

Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington................54

23. Graph showing measured hydraulic heads as a function of simulated heads,

Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington................57

24. Map showing differences between simulated and measured water levels

(residuals), Columbia Plateau Regional Aquifer System, Idaho, Oregon, and

Washington..................................................................................................................................58

25. Map showing well groups with similar hydraulic response within the Umatilla

subarea, Oregon .........................................................................................................................59

Figures—Continued
vi

26. Graphs showing simulated water levels for the North-South; and East-West well

groups near the Oregon Water Resources Department administrative areas in the

Umatilla subarea, Oregon..........................................................................................................60

27. Map showing well groups with similar hydrologic responses and generalized

groundwater-flow paths under 2000–10 conditions near the Washington State

Department of Ecology administrative areas in parts of the Palouse Slope/eastern

Yakima Fold Belt and the Columbia Basin Ground Water Management Area,

Washington..................................................................................................................................63

28. Graphs showing simulated water levels for the western flow path; eastern flow

path; middle flow path; and southern flow path well groups near the Washington

State Department of Ecology administrative areas in parts of the Palouse

Slope/eastern Yakima Fold Belt and the Columbia Basin Ground Water

Management Area, Washington..............................................................................................64

29. Map showing change in simulated base flow from predevelopment (pre-1900) to

current conditions (2000–2007) for the Columbia Plateau Regional Aquifer System,

Idaho, Oregon, and Washington...............................................................................................68

30. Maps showing simulated composite hydraulic head for Overburden unit, Saddle

Mountains unit, Wanapum unit, and upper Grande Ronde unit for the Columbia

Plateau Regional Aquifer System, Idaho, Oregon, and Washington .................................70

31. Maps showing simulated composite drawdowns from predevelopment to current

conditions (2000–2007) for Overburden unit, Saddle Mountains unit, Wanapum unit,

and (upper Grande Ronde unit for the Columbia Plateau Regional Aquifer System,

Idaho, Oregon, and Washington...............................................................................................74

32. Graph showing simulated annual water-budget flux for predevelopment to 2007 for

the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington.........79

33. Graph showing simulated annual water-budget flux for predevelopment, wet (1997),

average (2000), and dry (2001) years, and current conditions (2000–2007), Columbia

Plateau Regional Aquifer System, Idaho, Oregon, and Washington .................................79

34. Maps showing simulated effects of commingling wells only (no pumping); pumping

wells only (no commingling); and groundwater recharge from irrigation in the

Wanapum basalt unit, Columbia Plateau Regional Aquifer System, Idaho, Oregon,

and Washington..........................................................................................................................81

Figures—Continued

Tables

1. Correlation chart showing relation between generalized stratigraphy and

hydrogeologic units, Columbia Plateau Aquifer System, Idaho, Oregon, and

Washington....................................................................................................................................9

2. Correlation between geologic model units and groundwater-flow model layers of

the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington.........18

3. Mean annual base-flow estimates, and upper and lower bounds used to calibrate

the steady-state and transient models of the Columbia Plateau Regional Aquifer

System, Idaho, Oregon, and Washington ...............................................................................42

4. Final calibrated model parameters for the Columbia Plateau Regional Aquifer

System, Idaho, Oregon, and Washington ...............................................................................48

5. Simulated steady-state base flow and errors, Columbia Plateau Regional Aquifer

System, Idaho, Oregon, and Washington ...............................................................................55
vii

Conversion Factors, Datums, and Abbreviations and 

Acronyms

Conversion Factors

Inch/Pound to SI

Multiply By To obtain

Length

inch (in.) 2.54 centimeter (cm)

foot (ft) 0.3048 meter (m)

mile (mi) 1.609 kilometer (km)

Area

acre 4,047 square meter (m2)

square mile (mi2) 2.590 square kilometer (km2) 

Volume

cubic mile (mi3) 4.168 cubic kilometer (km3) 

acre-foot (acre-ft) 1,233 cubic meter (m3)

Flow rate

acre-foot per year (acre-ft/yr) 1,233 cubic meter per year (m3/yr)

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

cubic foot per second (ft3/s) 0.02832 cubic meter per second (m3/s)

inch per year (in/yr) 25.4 millimeter per year (mm/yr)

Hydraulic conductivity

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

Transmissivity*

foot squared per day (ft2/d) 0.09290 meter squared per day (m2/d) 

SI to Inch/Pound

Multiply By To obtain

Length

kilometer (km) 0.6214 mile (mi)

*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.

Datums

Vertical coordinate information is referenced to North American Vertical Datum of 1988 

(NAVD 88).

Horizontal coordinate information is referenced to the North American Datum of 1983 (NAD 83).

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

Conversion Factors, Datums, and Abbreviations and 

Acronyms

Abbreviations and Acronyms

BCF Block-Centered Flow package in MODFLOW

CBIP Columbia Basin Irrigation Project

CONDFACT conductance factor

CPRAS Columbia Plateau Regional Aquifer System

CRBG Columbia River Basalt Group

DRN Drain package in MODFLOW

ELEVATION drain elevation

ET evapotranspiration

GWMA Ground Water Management Area

HFB Horizontal-Flow Barrier

HUF Hydrogeologic-Unit Flow package in MODFLOW

HYSEP hydrograph separation technique

Kh horizontal hydraulic conductivity

Kv vertical hydraulic conductivity

LPF Layer Property Flow package in MODFLOW

MNW2 Multi-Node Well package in MODFLOW

NHD National Hydrography Dataset

MODFLOW USGS modular three-dimensional finite-difference groundwater-flow model

NWT Newton Solver package in MODFLOW

PEST parameter estimation software package

PRISM Parameter-elevation Regressions on Independent Slopes Model

RASA Regional Aquifer-System Analysis program

RBOT river bottom elevation

RIV River package in MODFLOW

RMS root-mean-square error

SOWAT SOil WATer balance

Ss specific storage

STAGE river stage

Sy specific yield

UPW Upstream Weighting package in MODFLOW-2005

USGS U.S. Geological Survey
Numerical Simulation of Groundwater Flow in the 

Columbia Plateau Regional Aquifer System, Idaho, 

Oregon, and Washington

By D. Matthew Ely, Erick R. Burns, David S. Morgan, and John J. Vaccaro

Abstract

A three-dimensional numerical model of groundwater 

flow was constructed for the Columbia Plateau Regional 

Aquifer System (CPRAS), Idaho, Oregon, and Washington, 

to evaluate and test the conceptual model of the system and 

to evaluate groundwater availability. The model described in 

this report can be used as a tool by water-resource managers 

and other stakeholders to quantitatively evaluate proposed 

alternative management strategies and assess the long-term 

availability of groundwater. The numerical simulation 

of groundwater flow in the CPRAS was completed with 

support from the Groundwater Resources Program of the 

U.S. Geological Survey Office of Groundwater.

The model was constructed using the U.S. Geological 

Survey modular three-dimensional finite-difference 

groundwater-flow model, MODFLOW-NWT. The model uses 

3-kilometer (9,842.5 feet) grid cells that subdivide the model 

domain by 126 rows and 131 columns. Vertically, the model 

domain was subdivided into six geologic model units. From 

youngest to oldest, the units are the Overburden, the Saddle 

Mountains Basalt, the Mabton Interbed, the Wanapum Basalt, 

the Vantage Interbed, and the Grande Ronde Basalt.

Natural recharge was estimated using gridded historical 

estimates of annual precipitation for the period 1895–2007. 

Pre-development recharge was estimated to be the average 

natural recharge for this period. Irrigation recharge and 

irrigation pumping were estimated using a remote-sensing 

based soil-water balance model for the period 1985–2007. 

Pre-1985 irrigation recharge and pumping were estimated 

using previously published compilation maps and the history 

of large-scale irrigation projects. Pumping estimates for 

municipal, industrial, rural, residential, and all other uses were 

estimated using reported values and census data. Pumping was 

assumed to be negligible prior to 1920.

Two models were constructed to simulate groundwater 

flow in the CPRAS: a steady-state predevelopment model 

representing conditions before large-scale pumping 

and irrigation altered the system, and a transient model 

representing the period 1900–2007. Automated parameterestimation techniques (steady-state predevelopment model) 

and traditional trial-and-error (transient model) methods were 

used for calibration. To calibrate the steady-state and transient 

models, 10,525 and 46,460 water level measurements, 

respectively, and 50 base-flow estimates were used.

The steady-state model simulated the shape, slope, and 

trends of a potentiometric surface that was generally consistent 

with mapped water levels. For the transient model, the mean 

and median difference between simulated and measured 

hydraulic heads is -10 and 4 ft, respectively, with a standard 

deviation of 164 ft over a 5,648 ft range of measured heads. 

The residuals for the simulation period show that 52 percent of 

the simulated heads exceeded measured heads with a median 

residual value of 43 ft, and 48 percent were less than measured 

heads with a median residual value of -76 ft.

The CPRAS model was constructed to derive 

components of the groundwater budget and help understand 

the interactions of stresses, such as recharge, groundwater 

pumping, and commingling wells on the groundwater and 

surface-water system. Through these applications, the model 

can be used to identify trends in groundwater storage and 

use, and quantify groundwater availability. The annual 

groundwater budgets showed several patterns of change over 

the simulation period. Groundwater pumping was negligible 

until the 1950s and began to increase significantly during the 

1970s and 1980s. Recharge was highly variable due to the 

interannual variability of precipitation, but began to increase 

in the late 1940s due to the increase in surface-water irrigation 

projects. Groundwater contributions to streamflow (base flow) 

followed recharge closely. However, in areas of significant 

groundwater-level decline, base flow is reduced.

Groundwater pumping had the greatest effect on 

water levels, followed by irrigation enhanced recharge. 

Commingling was a larger factor in structurally complex 

upland areas where hydraulic-head gradients are 

naturally high.
2 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Groundwater pumping has increased substantially over 

the past 40–50 years; this increase resulted in declining water 

levels at depth and decreased base flows over much of the 

study area. The effects of pumping are mitigated somewhat 

by the increase of surface-water irrigation, especially in the 

shallow Overburden unit, and commingling wells in some 

areas. During dry to average years, groundwater pumping 

causes a net loss of groundwater in storage and current 

condition (2000–2007) groundwater pumping exceeds 

recharge in all but the wettest of years.

Introduction

The Columbia Plateau Regional Aquifer System 

(CPRAS) covers approximately 44,000 mi2 of northeastern 

Oregon, southeastern Washington, and western Idaho (fig. 1). 

The area supports a $6 billion-per-year agricultural industry, 

leading the Nation in production of apples and nine other 

commodities (State of Washington Office of Financial 

Management, 2014; U.S. Department of Agriculture, 2007). 

Groundwater availability in the area is a critical water-resource 

management issue due to the large water demand for 

agriculture, economic development, and ecological needs.

The primary aquifers of the CPRAS are basalts of 

the Columbia River Basalt Group (CRBG) and in places, 

overlying basin-fill sediments (Overburden). Waterresources issues that have implications for current (2013) 

and future groundwater availability in the region include 

(1) widespread water-level declines associated with 

withdrawals of groundwater for irrigation and other uses, 

(2) decreases in base flow to rivers and the associated effects 

on river temperature and water quality, and (3) current and 

potential effects of global climate change on recharge, base 

flow, and ultimately, groundwater availability.

The U.S. Geological Survey (USGS) Groundwater 

Resources Program began a study of the CPRAS in 2007 

with the broad goals of (1) characterizing the hydrologic 

status of the system, (2) identifying trends in groundwater 

storage and use, and (3) quantifying groundwater availability. 

The study approach included updating and refining the 

regional hydrogeologic framework, documenting changes 

in the status of the system, quantifying the hydrologic 

budget, and developing a groundwater-flow model for the 

system. The simulation model, presented here, was used to 

evaluate and test the conceptual model of the system and to 

evaluate groundwater availability. Groundwater availability 

is not only a function of the quantity and quality of water 

in an aquifer system, but also the physical structures, laws, 

regulations, and socioeconomic factors that control its demand 

and use (Reilly and others, 2008). This report discusses the 

physical characteristics that are important as indicators of 

groundwater availability.

Purpose and Scope

This report describes the construction, calibration, and 

application of a numerical model of groundwater flow in the 

Columbia Plateau Regional Aquifer system. The hydrogeologic 

framework for constructing the model was documented in a 

previous report (Kahle and others, 2011) as part of this study. 

The purposes for constructing the model were to test the 

conceptual model and to provide an improved understanding of 

the groundwater-flow system and groundwater availability. The 

model development is presented and described, and includes 

information on the spatial and temporal discretization of the 

aquifer system, boundary conditions, stresses, and hydraulic 

properties of the hydrogeologic units constituting the aquifer 

system. Predevelopment and 1900–2007 conditions were 

simulated to provide a better understanding of current demands 

on groundwater in the study area.

Previous Investigations

Numerous geologic, hydrologic, and hydrogeologic studies 

of or within the CPRAS have been done. The earliest works of 

Smith (1901), Calkins (1905), Waring (1913), Schwennsen and 

Meinzer (1918), Piper (1932), and Taylor (1948) formed the 

foundation of our current understanding of the water resources 

of the Columbia Plateau. A substantial body of work on the 

geology and hydrology of the CPRAS was produced as part of 

the USGS Regional Aquifer-System Analysis (RASA) Program 

of the 1980s and 1990s. A description of the hydrogeologic 

framework, characteristics of the hydrogeologic units, water 

budget components, geochemistry of the aquifer system, and 

regional groundwater-flow system are provided in Whiteman 

and others (1994).

As part of the USGS groundwater availability study, 

aspects of the CPRAS, such as the hydrogeology and waterbudget components, have been investigated and provided 

the basis for much of the updated information incorporated 

in the groundwater-flow model. The model presented in this 

report is based on the existing USGS groundwater-flow model 

(CP-RASA) developed by Hansen and others (1994) for 

the RASA program. The boundaries and much of the initial 

conceptual understanding of the current model is based on the 

CP-RASA model. Groundwater-flow models exist for subareas 

of the CPRAS, including the Columbia Basin Irrigation Project 

(CBIP) area (Tanaka and others, 1974), Umatilla and Horse 

Heaven Hills (Davies-Smith and others, 1988), PullmanMoscow (Barker, 1979; Lum and others, 1990), Horse Heaven 

Hills (Packard and others, 1996), Hanford Site (Bergeron and 

others, 1986; Wurstner and others, 1995; Vermeul and others, 

2001 and 2003), Yakima River Basin (Ely and others, 2011), 

Mosier, Oregon (Burns, 2012a), and the Columbia Basin 

Ground Water Management Area (GWMA) Adams, Franklin, 

Grant, and Lincoln Counties, Washington (Porcello and 

others, 2010).
Introduction 3

tac13-0856_fig01

Grand Coulee Dam

B L U E

M O U N T A I N S

CASCADE RANGE

OKANOGAN HIGHLANDS

CASCADE RANGE

WALLOWA

MOUNTAINS

R O C K Y M O U N T A I N S 

R O C K Y M O U N T A I N S

HORSE HEAVEN HILLS

BEEZLEY HILLS

Snohomish

Skagit

King

Pierce

Lewis

Chelan

Okanogan

Kittitas

Yakima

Skamania

Klickitat

Hood

River

Clackamas

Marion

Linn

Ferry

Stevens

Pend

Orielle

Benton

Grant

Lincoln

Spokane

Adams

Franklin

Walla Walla

Whitman

Columbia

Garfield

Asotin

Sherman

Jefferson

Crook

Gilliam

Morrow

Umatilla

Wheeler

Malheur

Baker

Wallowa

Bonner

Kootenai

Benewah

Latah

Lincoln

Sanders

Clearwater

Nez Perce

Washington Valley

Adams

Idaho

Douglas

Wasco

Grant

Union

Shoshone

Lewis

Snake

River

Crab Creek

Columbia

River

COLUMBIA

Naches River

Yakima River

Umatilla River

Columbia

River

RIVER

N orth Fork John Day River

Snake River

Imnaha River

Banks

Lake Coeur

D'Alene

Lake

Grande Ronde River

Salmon River

John D ay River

Moses Coulee

Walla Walla River

Deschutes River

Klickitat River

White Salmon R

Clearwater

River

Hangman Cr

F.D Ro osevelt La ke

Hood R

Wenatch e e R.

Palouse River

Satus

T

Cr

oppenish Creek

Spokane

River

Potholes 

Reservoir

Moses

Lake

Rocky

Ford Cr.

WASHINGTON

OREGON IDAHO

Location of

study area

Ellensburg

Umatilla

Walla Walla

WASHINGTON

OREGON

IDAHO

Yakima

Heppner

MiltonFreewater

Pasco

Lewiston

Moscow

Spokane

The

Dalles

LaGrande

122° 120° 118° 116°

46°

48°

47°

45°

0 20 40 60 80 Miles

0 20 40 60 80 Kilometers

Base modified from U.S. Geological Survey digital data, various scales. 

EXPLANATION

Structural regions of Columbia 

Plateau Regional Aquifer System

Yakima Fold Belt

Palouse Slope

Blue Mountains

Clearwater Embayment

Study area boundary

Groundwater-flow model boundary

Figure 1. Columbia Plateau Regional Aquifer System study area and structural regions, Idaho, Oregon, and Washington. Structural 

regions modified from Reidel and others (2002).
4 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Description of Study Area

A complete description of the Columbia Plateau is 

available in Kahle and others (2011), parts of which are 

presented here.

Location and Setting

The Columbia Plateau is a structural and topographic 

basin within the drainage of the Columbia River (fig. 1). It 

is bounded on the west by the Cascade Range, on the east 

by the Rocky Mountains, and on the north by the Okanogan 

Highlands. Its southern boundary corresponds to the mapped 

extent of the CRBG. The Columbia Plateau is underlain by 

massive basalt flows having an estimated composite thickness 

of at least 14,000 ft at one of the lowest points of the plateau 

near Pasco, Washington (Drost and others, 1990; Reidel 

and others, 2002). Sedimentary deposits overlie the basalt 

over large areas of the plateau, and may exceed 2,000 ft in 

thickness in places.

The Columbia Plateau was divided into four informal 

structural regions—the Yakima Fold Belt, Palouse Slope, 

Blue Mountains, and the Clearwater Embayment (fig. 1). The 

Yakima Fold Belt includes most of the western half of the 

plateau and is characterized by a series of east-west trending 

anticlinal ridges and synclinal basins. The Palouse Slope, in 

the northeast quarter of the plateau, is much less deformed 

and has a gently southwestward dipping slope. The other 

structural regions within the CPRAS are the Blue Mountains, 

a composite anticlinal structure that forms the southern extent 

of the Columbia River Basin, and the Clearwater Embayment, 

which marks the eastward extent of the CPRAS along the 

foothills of the Rocky Mountains and includes a series of folds 

extending into Idaho. The geologic and hydraulic properties 

of these structural regions may influence groundwater 

flow owing to flow barriers which may impede flow or to 

compartmentalization (the creation of distinct zones within an 

aquifer with limited interconnectivity) of the hydrogeologic 

units. The presence and importance of flow barriers and 

compartments in the CPRAS have been recognized and 

discussed in numerous studies (for example, Newcomb, 1969; 

Porcello and others, 2010).

Much of the Columbia Plateau is semiarid, the mean 

annual precipitation for 1895–2007 is about 17 in. (about 

40 million acre-ft) and ranges from about 7 in. in the center 

of the study area to more than 60 in. in the northwesternmost extent of the study area (PRISM Climate Group, 2004; 

calculated from annual values). The types and amounts of 

natural vegetation on the Columbia Plateau vary according to 

precipitation and land-surface altitudes. The vegetation ranges 

from sagebrush and grasslands at lower altitudes to grasslands 

and forest at mid-altitudes to barren rock and conifer forests 

that are representative of the mountainous topography at the 

upper altitudes. Dry land agriculture mostly includes winter 

and spring wheat and lentils. Irrigated agriculture includes 

apples, hops, and other crops.

Overviews of the geology and hydrology of the CPRAS 

presented in this report summarize detailed descriptions 

in reports by (1) Kahle and others (2009), who discuss the 

geologic framework used in this report; (2) Burns and others 

(2011), who describe the three-dimensional characteristics of 

the geology of the CPRAS; and (3) Kahle and others (2011), 

who discuss the hydrogeologic framework and the hydrologic 

budget components of the CPRAS.

History of Water Resource Development

The cultural and economic development of the Columbia 

Plateau has depended heavily on the availability of water and 

the ability to store and redistribute water from the Columbia 

River and its major tributaries. Irrigation began as early as the 

1840s and 1850s at missions in Walla Walla, Lewiston, and the 

Yakima Valley. In the early 1900s, the lumber and agricultural 

industries grew steadily and small-grain production on dryland 

farms and dairy and poultry farming were especially profitable 

because these enterprises did not need large quantities 

of water.

The dry summer climate in the region forced the early 

settlers to develop water supplies for irrigation wherever 

possible. By the 1930s, economic growth was relatively slow 

because of the 1929 depression and severe droughts that 

occurred after 1919. Those who survived these hardships lived 

along surface-water bodies or in areas where groundwater was 

available at shallow depths. A notable exception to this pattern 

was the Yakima River Basin, where reservoirs, diversion 

dams, and canals were constructed during 1892–1933. By 

1902, about 120,000 acres were under mostly surface-water 

irrigation in the Yakima River basin (Parker and Storey, 1916; 

Bureau of Reclamation, 1999). Irrigation water projects by 

the Bureau of Reclamation (Reclamation) between 1910 

and 1933 allowed the irrigated acreage to grow to more than 

500,000 acres.

The start of Reclamation’s CBIP (fig. 2) in 1933 and 

construction of the Hanford Site (between the Yakima and 

Columbia Rivers in Benton County, Washington [fig. 1]) in 

the 1940s brought a large influx of workers and associated 

service industries. By 1946, irrigation water was supplied to 

about 850,000 acres within the study area (Simons, 1953). 
Description of Study Area 5

tac13-0856_fig 02

Base map modified from USGS and other digital data, various 

scales. Horizontal datum: North American Datum of 1983.

B L U E

M O U N T A I N S

121° 120° 119° 118° 117° 116°

48°

47°

46°

45°

WASHINGTON

OREGON IDAHO

Location of

study area

EXPLANATION

Columbia Plateau Regional 

 Aquifer System

Columbia Basin Irrigation Project

Ground Water Management Area

Odessa Subarea

Umatilla Subarea

Walla Walla River Basin

Yakima River Basin

0 50 100 MILES

0 50 100 KILOMETERS

COLUMBIA RIVER

Figure 2. Location of major river basins and study subareas, Oregon and Washington.
6 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Most of these acres were irrigated with water supplied by 

Reclamation projects created under the Federal Reclamation 

Act of 1902. Water from Franklin D. Roosevelt Lake, the 

reservoir formed by Grand Coulee Dam, became available 

in 1952, and by 1972 more than 0.5 million acres were being 

irrigated by this project in Washington, producing more than 

60 crops. Currently (2013), about 2 million acres of croplands 

are irrigated with surface water and groundwater; much of 

the surface water is supplied by Reclamation projects. The 

surface-water withdrawals account for about 75 percent of the 

total irrigation water and groundwater withdrawals account for 

the remaining 25 percent (Kahle and others, 2011).

With the advent of new technology in about 1950, a 

rapid and intensive expansion of deep-well irrigation practices 

took place in areas not served by surface-water irrigation 

projects. These areas included parts of the Yakima, Pasco, 

Umatilla, The Dalles, and Walla Walla River Basins and the 

Odessa subarea in western Adams County, Washington. By 

1984, about 0.5 million acres were irrigated with groundwater 

(Whiteman and others, 1994. By 2002, the total number of 

irrigated acres generally had stabilized (U.S. Department of 

Agriculture, 2007).

Hydrogeologic Setting

The Columbia Plateau is an intermontane basin between 

the Rocky Mountains and the Cascade Range that is filled with 

mostly Cenozoic basalt and sediment. Most rocks exposed 

in the region are the CRBG, intercalated sedimentary rocks, 

and younger sedimentary rocks and deposits that include 

Pleistocene cataclysmic flood deposits, eolian deposits, and 

terrace gravels of modern rivers. The CRBG consists of a 

series of more than 300 lava flows that erupted during various 

stages of the Miocene Age, 17 million to 6 million years ago. 

Individual flows range in thickness from 10 to more than 

300 ft (Tolan and others, 1989; Drost and others, 1990) with 

a total thickness that might be greater than 14,000 ft in the 

central part of the study area near Pasco, Washington (Reidel 

and others, 2002). Soils derived from flows or sediments 

deposited on the surface of a flow were sometimes preserved, 

creating sedimentary interbeds between flows.

Generalized hydrogeologic units of the Overburden and 

CRBG, listed in order of generally increasing age, include 

Overburden, Saddle Mountains, Mabton Interbed, Wanapum, 

Vantage Interbed, and Grande Ronde (figs. 3 and 4; table 1). 

The Overburden unit consists predominantly of undivided, 

unconsolidated to semi-consolidated sedimentary deposits 

ranging from Pliocene to Holocene in age (Drost and others, 

1990). These include many types of deposits of local and 

(or) regional extent including flood gravels and slack water 

sediments, terrace gravels of modern rivers, and eolian 

deposits that can range in thickness from 0 to 1,300 ft. The 

Saddle Mountains hydrogeologic unit consists mostly of the 

Saddle Mountains Basalt and interbed members and is the 

least extensive and youngest formation of the CRBG. Most of 

the unit is in the west-central part of the study area, with less 

continuous occurrences in the Blue Mountains and eastward 

into Idaho. Thickness of the Saddle Mountains unit can range 

from 0 to about 1,000 ft. The Wanapum hydrogeologic unit, 

composed mostly of basalt and interbed members of the 

Wanapum Basalt, mostly is in the north-central part of the 

study area. Much of the unit lies beneath the Overburden 

and Saddle Mountains units. Thickness of the Wanapum 

unit ranges from 0 to about 1,200 ft. The Grande Ronde 

hydrogeologic unit is the oldest and most extensive of the 

basalt units and constitutes the vast majority of the CRBG. 

This unit underlies most of the study area, except for an area 

along the southern boundary of the CPRAS in Oregon and 

along the eastern edge of the CPRAS in Idaho. The Grande 

Ronde unit contains the Grande Ronde Basalt and associated 

interbed members. Younger Wanapum and Saddle Mountains 

basalts were deposited in synclines and fault bounded valleys 

underlain by Grande Ronde basalts, indicating folding and 

faulting was active during CRBG deposition, continuing into 

the Quaternary. Distribution and thickness of younger basalt 

flows are determined by structurally controlled valleys. During 

the Pleistocene, the surface of the basalt units was modified 

greatly during repeated catastrophic outburst flooding, 

which caused erosion of vast channels as well as removal or 

deposition of overlying sediment. The basement confining 

unit, referred to as Older Bedrock, is composed of various 

rock types older than the CRBG (Kahle and others, 2009). In 

Washington and Idaho, the rocks bordering the CPRAS consist 

mostly of sedimentary and granitic rocks. In Oregon, the 

CPRAS is bordered by sedimentary, volcaniclastic, volcanic, 

plutonic, and metamorphic rocks (Drost and others, 1990). The 

bedrock unit generally has much lower permeabilities than the 

basalts and is considered the base of the regional flow system.

Hydrologic Budget Components

Discussion and estimates of the regional-scale hydrologic 

budget components for the CPRAS is from Kahle and others 

(2011). The extent of the study area used by Kahle and others 

(2011) differs from the model domain presented in this report; 

therefore, a direct comparison of budget component values 

with simulated values is not appropriate. The conceptual 

model for the aquifer system is that water (1) enters the 

system as recharge from precipitation (rainfall and snowmelt) 

and recharge from the delivery and application of surfacewater and groundwater irrigation; and (2) exits the system as 

streamflow, evapotranspiration, and groundwater pumpage. A 

complete description of the budget components can be found 

in Kahle and others (2011).
Description of Study Area 7

tac13-0856_fig 03

Base map modified from USGS and other digital data, various 

scales. Horizontal datum: North American Datum of 1983.

From Kahle and others (2009)

B L U E

M O U N T A I N S

A

A'

B

B'

River

Snake

121° 120° 119° 118° 117° 116°

48°

47°

46°

45°

WASHINGTON

OREGON IDAHO

Location of

study area

EXPLANATION

Columbia Plateau Regional 

 Aquifer System

Well 

Deepened well 

A A' Hydrogeologic section

Well used in construction

 of section 

Generalized hydrogeologic units

Overburden

Saddle Mountains

Wanapum

Grande Ronde

Older Bedrock

0 50 100 MILES

0 50 100 KILOMETERS

Columbia River

Figure 3. Surficial distribution of generalized hydrogeologic units and locations of hydrogeologic sections, Columbia Plateau 

Regional Aquifer System, Idaho, Oregon, and Washington. Hydrogeologic sections are shown in figure 4.
8 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

tac13-0856_fig 04

Columbia River

Snake River

Snake River

Snake River

SECTION

A–A'

SECTION

B–B'

FEET

A

FEET

B'

SOUTH

FEET METERS

B

NORTH

FEET

A'

METERS

WEST EAST

4,000

2,000

3,000

1,000

2,000

4,000

3,000

5,000

1,000

7,000

9,000

8,000

10,000

6,000

0

NORTH AMERICAN

VERTICAL DATUM

OF 1988

6,000

5,000

4,000

2,000

3,000

1,000

2,000

4,000

3,000

5,000

1,000

7,000

9,000

8,000

6,000

0

1,000

1,000

2,000

3,000

6,000 1,800

5,000

4,000

2,000

3,000

1,000

2,000

4,000

3,000

5,000

1,000

7,000

9,000

8,000

10,000

6,000

0

NORTH AMERICAN

VERTICAL DATUM

OF 1988

11,000

12,000

14,000

16,000

15,000

17,000

13,000

4,000

2,000

3,000

1,000

2,000

4,000

3,000

5,000

1,000

7,000

9,000

8,000

10,000

6,000

0

11,000

12,000

14,000

16,000

15,000

17,000

13,000

1,000

1,000

2,000

3,000

4,000

5,000

Contact

Contact, approximately located

Section trace intersection

EXPLANATION

Generalized hydrogeologic units

Overburden

Saddle Mountains Vantage Interbed

Mabton Interbed

Wanapum

Grande Ronde

Older Bedrock

VERTICAL SCALE GREATLY EXAGGERATED (X 40)

0 25 50 75 MILES

0 25 50 75 KILOMETERS

Well

Figure 4. Generalized hydrogeologic sections of the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and 

Washington. Adapted from Kahle and others (2009).
Description of Study Area 9

Table 1. Correlation chart showing relation between generalized stratigraphy and hydrogeologic units, Columbia Plateau Aquifer 

System, Idaho, Oregon, and Washington.

[From Kahle and others, 2009]

Era Period Epoch Sediment stratigraphy Basalt stratigraphy Hydrogeologic

unit

Cenozoic

Quaternary

Holocene Alluvial, colluvial, eolian, 

glacial, glacial outburst 

flood, lacustrine, landslide, 

terrace, and peat deposits; 

ash, debris-avalanche and 

debris‐flow deposits, talus;

Touchet Beds, Palouse 

Formation Quaternary and Pliocene Basalts Overburden

Pleistocene

Tertiary

Pliocene

Alluvial fan deposits; 

Alkali Canyon, Chenoweth, 

Deschutes, Madras and 

Ringold Formations; Dalles 

Group; Thorpe Gravel; 

and unknown continental 

sedimentary deposits

Miocene

Ellensburg, Deschutes, 

Latah, Madras, Payette, and 

Ringold Formations; Dalles 

Group; Snipes Mountain 

deposits; Deer Creek Beds; 

and unknown continental 

sedimentary deposits

Columbia 

River 

Basalt 

Group

Saddle Mountains Basalt flow 

members and interbeds

Saddle Mountains 

unit

Mabton interbed (Mabton Member 

of the Ellensburg Formation) Mabton unit

Wanapum Basalt Flow

members and interbeds Wanapum unit

(Vantage Member of the 

Ellensburg Formation) Vantage unit

Grande Ronde 

Basalt flow 

members and 

interbeds

Grande Ronde 

unit

Prineville

Basalt

Picture 

Gorge 

Basalt

Imnaha Basalt

pre-Columbia River Basalt Group rocks, undivided Older Bedrock
10 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Mean annual recharge from infiltration of precipitation 

was estimated on the basis of annual precipitation and the 

results of previous modeling of recharge (Bauer and Vaccaro, 

1990). The other budget components, excluding streamflow, 

were developed using estimates of actual evapotranspiration 

from a Simplified Surface Energy Balance method (Senay and 

others, 2007) and a monthly SOil WATer balance (SOWAT; 

Kahle and others, 2011) model that determined irrigation 

water demand, flux into and out of the groundwater system 

(recharge or discharge), direct runoff, and soil moisture in 

irrigated areas. The SOWAT model used simple relations 

among climatic, soils, land-cover, and irrigation data to 

compute monthly irrigation requirements and surplus moisture 

available for recharge. Estimates of groundwater pumping for 

irrigation and recharge from irrigation return flow from the 

application of the SOWAT model were refined as initial input 

to the groundwater-flow model.

The estimated mean annual recharge from infiltration 

of precipitation for the CPRAS was 4.6 in/yr (10.8 million 

acre-feet per year [acre-ft/yr]) for 1985–2007. The spatial 

distribution in recharge mirrors that of annual precipitation, 

with the highest mean annual recharge (more than 20 in.) 

occurring in the upper Yakima River Basin, the Blue 

Mountains southeast of Walla Walla, and adjacent to the 

Columbia River where it leaves the study area through the 

gorge in the Cascade Range. Mean annual recharge from 

infiltration of precipitation is less than 1 in. in a large part 

of the study area where precipitation is low, adjacent to the 

Columbia and Yakima Rivers.

Annual irrigation water use in the study area averaged 5.3 

million acre-feet (acre-ft) during 1985–2007, with 1.4 million 

acre-ft (26 percent) supplied from groundwater and 3.9 

million acre-ft (74 percent) supplied from surface water. Mean 

monthly irrigation throughout the study area peaks in July 

at 1.6 million acre-ft (1985–2007 average), of which, 0.45 

and 1.15 million acre-ft is from groundwater and surfacewater sources, respectively. Mean annual recharge from 

irrigation return flow in the study area was 4.2 million acre-ft 

(1985–2007) with 2.1 million acre-ft (50 percent) occurring 

within the predominately surface-water irrigated regions of 

the Yakima and Umatilla River Basins and Columbia Basin 

Irrigation Project. Annual recharge rates range from less 

than 5 in/yr in predominately sprinkler-irrigated areas where 

groundwater is the source to more than 20 in/yr in surfacewater supplied areas where conveyance losses and lessefficient application methods are common.

Annual groundwater use (1985–2007) for purposes 

other than irrigation was estimated for the study area using 

information from multiple sources (Kahle and others, 

2011). Public-supply groundwater use increased from about 

201,500 acre-ft/yr in 1985 (12 percent of total groundwater 

use) to about 262,800 acre-ft/yr in 2007 (10 percent of total 

groundwater use). Domestic self-supplied groundwater use 

increased from about 54,600 acre-ft/yr in 1985 (3 percent 

of total groundwater use) to about 70,100 acre-ft/yr in 2007 

(3 percent of total groundwater use). Industrial groundwater 

use decreased from about 53,400 acre-ft/yr in 1985 (3 percent 

of total groundwater use) to about 43,900 acre-ft/yr in 2007 

(2 percent of total groundwater use). Other groundwater 

use, including water used for mining, thermoelectric needs, 

livestock, and aquaculture combined, increased from 

16,900 acre-ft/yr (1 percent of total groundwater use) in 

1985 to about 43,600 acre-ft/yr in 2007 (2 percent of total 

groundwater use).

Description of Groundwater-Flow 

System

Hydrogeologic Units

The conceptual groundwater model developed for the 

study area divides the aquifer system into seven hydrogeologic 

units (table 1)—the Overburden unit, three aquifer units 

in the permeable basalt rock, two sedimentary interbeds, 

and the basement confining unit (older bedrock) (Vaccaro, 

1999). The three basalt hydrogeologic units are the Saddle 

Mountains, Wanapum, and Grande Ronde Basalts and 

their intercalated sediments. To distinguish these basalt 

hydrogeologic units from basalt formations in this study, they 

are referred to as “units.” For example, the Wanapum Basalt 

and intercalated sediments are referred to as the Wanapum 

unit (Vaccaro, 1999). In the southeastern part of the study 

area, the Imnaha Basalt and any intercalated sediments are 

included with the Grande Ronde unit. Similarly, the Picture 

Gorge and Prineville Basalts in the southern part of the study 

area are included in the Grande Ronde unit. The interbeds 

between the Saddle Mountains and Wanapum basalt units, 

and between the Wanapum and Grande Ronde basalt units, 

are referred to in this study as the Mabton and Vantage 

Interbeds, respectively (Kahle and others, 2009, table 1). The 

interbed units are fairly extensive laterally, but are thin when 

compared with the thickness of the basalt units. The basement 

confining unit, referred to as “Older Bedrock,” consists of 

pre-CRBG sedimentary, volcaniclastic, volcanic, plutonic, 

and metamorphic rocks that generally have much lower 

permeabilities than the overlying basalts and is considered 

the base of the regional flow system (Vaccaro, 1999; Kahle 

and others, 2009). The approximate surficial distribution of 

overburden and the three basalt hydrogeologic units are shown 

in figure 3. The approximate subsurface distribution of the 

CPRAS hydrogeologic units is shown in figure 4. Detailed 

descriptions of the units are available in Drost and others 

(1990), Whiteman and others (1994), Vaccaro (1999), and 

Jones and Vaccaro (2008), Kahle and others (2009), and Burns 

and others (2011).
Description of Groundwater-Flow System 11

Overburden Unit

The Overburden unit consists of undivided, 

unconsolidated to semi-consolidated sedimentary deposits and 

minor basalt and andesite, ranging from Pliocene to Holocene 

in age (Drost and others, 1990). The Overburden unit includes 

several formations of local or regional extent and numerous 

types of deposits including alluvial, colluvial, eolian, glacial, 

glacial outburst flood, lacustrine, landslide, terrace, and peat 

deposits; talus; and other continental sedimentary deposits of 

undetermined origin (table 1). Thickness of the Overburden 

unit ranges from 0 to 1,300 ft, with a median thickness of 

about 47 ft (Kahle and others, 2009).

Saddle Mountains Unit

The Saddle Mountains unit consists mostly of the Saddle 

Mountains Basalt and intercalated sediments. Most of the 

unit is in the west-central part of the study area, with less 

continuous occurrences in the Blue Mountains and eastward 

into Idaho (fig. 3). Kahle and others (2009) estimated an areal 

extent of about 8,000 mi2, with a range in altitude of the unit 

top from about 4,000 to -280 ft relative to North American 

Vertical Datum of 1988 (NAVD 88). Thickness of the Saddle 

Mountains unit ranges from about 0 to 990 ft, with a median 

thickness of 280 ft (Kahle and others, 2009).

Mabton Interbed Unit

The Mabton unit is the sedimentary interbed between the 

overlying Saddle Mountains unit and the underlying Wanapum 

unit. The Mabton unit consists of the Mabton Member of the 

Ellensburg Formation and is mostly in the west-central part 

of the study area. Outcrops of the Mabton unit are scarce in 

the study area and the extent of the Mabton unit is assumed 

to be within the extent of the Saddle Mountains unit. The 

Mabton unit generally consists of clay, shale, claystone, clay 

with basalt, clay with sand, and sandstone. Thickness of the 

Mabton unit ranges from about 0 to 520 ft, with a median 

thickness of about 44 ft (Kahle and others, 2009).

Wanapum Unit

The Wanapum unit, composed mostly of basalt and 

intercalated sediments of the Wanapum Basalt, is in most 

of the north-central part of the study area (fig. 3) and has an 

estimated areal extent of about 25,000 mi2 with the altitude 

of the top ranging from about 3,400 to -1,000 ft relative 

to NAVD 88 (Kahle and others, 2009). Much of the unit 

lies beneath the Overburden and Saddle Mountains units. 

Thickness of the Wanapum unit ranges from about 0 to 

1,200 ft, with a median thickness of about 330 ft (Kahle and 

others, 2009).

Vantage Interbed Unit

The Vantage unit is the sedimentary interbed between 

the overlying Wanapum unit and the underlying Grande 

Ronde unit. Over most of the study area, this unit consists of 

the Vantage Member of the Ellensburg Formation; however, 

this unit includes sediment of the Latah Formation in the 

northeastern part of the study area. Outcrops of this unit are 

scarce in the study area and the extent is assumed to be within 

the extent of the Wanapum unit (fig. 3). The Vantage unit 

consists of clay, shale, sandstone, tuff with claystone, and clay 

with basalt, but also may contain small amounts of sand and 

sand-and-gravel. A few well logs also indicate that the Vantage 

unit is not present in the southeastern part of the Yakima River 

Basin and near the Cold Creek Syncline and Rattlesnake Hills 

Structure (Jones and Vaccaro, 2008). Thickness of the Vantage 

unit ranges from about 0 to 320 ft, with a median thickness of 

about 20 ft (Kahle and others, 2009).

Grande Ronde Unit

The Grande Ronde unit is the oldest and most extensive 

of the basalt units. This unit underlies most of the study area, 

except for an area along the southern boundary of the CPRAS 

in Oregon and along the eastern extent in Idaho (fig. 3). The 

estimated areal extent of the Grande Ronde unit is about 

42,000 mi2 (Kahle and others, 2009). The Grande Ronde 

unit predominantly contains the basalt and interbed members 

associated with the Grande Ronde Basalt. Sedimentary 

interbeds within the unit generally are rare and are only a few 

feet thick where present. The top of the Grande Ronde unit 

ranges from 4,300 to -2,100 ft relative to NAVD 88 based 

on Kahle and others (2009). Thickness of the unit is largely 

unknown, but is estimated to be greater than 14,000 ft near the 

central part of the basin.

Older Bedrock Unit

The Older Bedrock unit that borders and underlies the 

CPRAS is composed of various rock types older than the 

CRBG (Kahle and others, 2009). In Washington and Idaho, 

the rocks bordering the CPRAS consist mostly of sedimentary 

and granitic rocks. In Oregon, the CPRAS is bordered 

by sedimentary, volcaniclastic, volcanic, plutonic, and 

metamorphic rocks (Drost and others, 1990).
12 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Hydraulic Characteristics of 

Hydrogeologic Units

The hydraulic characteristics of the aquifers determine 

how a groundwater-flow system functions and how it will 

respond to stresses such as pumping. These characteristics 

include horizontal and vertical hydraulic conductivity, 

porosity, and the storage coefficient. Because of the 

heterogeneity of the geologic materials that constitute the 

CPRAS, the hydraulic characteristics can vary widely. The 

overburden deposits are diverse in lithology and the large 

variation in grain size, depositional regimes, and age of 

the deposits account for the large range of their hydraulic 

characteristics (Kahle and others, 2011, p. 20). Hydraulic 

characteristics also vary greatly within and between the 

individual basalt flows and the basalt hydrogeologic units. 

Horizontal hydraulic conductivities generally are greatest 

in the interflow zones formed from the combination of 

basalt flow tops with the base of an overlying basalt 

flow (fig. 5), and an intervening sedimentary interbed, 

if present. The flow tops generally are brecciated and 

(or) vesicular, and the flow bases are brecciated and may 

contain pillow complexes if the basalt extruded within or 

flowed into water and exhibit high horizontal hydraulic 

conductivity (Lindholm and Vaccaro, 1988). The interflow 

zones commonly are separated by a less-transmissive flow 

interior (fig. 5) in which fractures typically are vertically 

oriented. The older bedrock generally has lower values of 

porosity and hydraulic conductivity than the overburden and 

CRBG units.

Kahle and others (2011, p. 21) estimated the median 

horizontal hydraulic conductivity values for the overburden, 

basalt units, and bedrock as 161, 70, and 6 ft/d, respectively, 

based on specific capacity data reported in previous studies. 

Vertical hydraulic conductivity of the hydrogeologic units 

in the CPRAS largely is unknown. Estimates of vertical 

hydraulic conductivity range from about 0.009 to 2 ft/d 

for the Overburden unit although values for some parts of 

this unit may be as low as 10-10 to 10-7 ft/d; values for the 

CRBG units range from 4×10-7 to 4 ft/d (Kahle and others, 

2011, p. 57).

The storage characteristics of an aquifer are described 

by the storage coefficient, which is defined as the volume 

of water that an aquifer takes into or releases from storage 

per unit surface area, per unit change in head. The storage 

coefficient is expressed in units of cubic feet per cubic feet, 

a dimensionless quantity. Previous estimates of the storage 

coefficient in the CPRAS typically range from about 0.1 

to 0.2 for the unconfined parts of the Overburden unit and 

from about 10-6 to 0.01 for the CRBG basalt units (Kahle 

and others, 2011, p. 26–27).

Groundwater Flow and Occurrence in Basalt Units

Groundwater moves through the regional aquifer system 

through preferential pathways (fig. 5) developed during lava 

deposition. Generally, each CRBG lava flow consists of a 

dense flow interior and irregular flow tops and flow bottoms 

with a variety of textures (Reidel and others, 2002). Although 

flow interiors have joints and fractures, they typically do not 

transmit water easily. Flow tops and bottoms commonly are 

vesicular or brecciated, and may or may not be permeable. Local 

permeability of flow tops and bottoms can be highly variable over 

short distances as a result of depositional processes, but tends 

to be high over long distances, resulting in highly transmissive 

aquifers at the regional scale. The CRBG thus comprises a stack 

of laterally extensive lava flows with relatively thin permeable 

productive zones at flow tops and flow bottoms separated by 

relatively thick dense flow interiors of low permeability. Thin 

permeable aquifers are estimated to occupy about 10 percent of 

the total flow thickness. Flow interiors have low permeability and 

low storage characteristics, so that they form effective confining 

units between permeable flow tops. As a result, the aquifer 

system is highly anisotropic, with effective horizontal hydraulic 

conductivity controlled by the fraction of the thickness occupied 

by thin aquifers, and the effective vertical hydraulic conductivity 

controlled by the dense flow interiors.

Except for groundwater flow in the deeply buried parts of the 

CPRAS, large-scale structural features compartmentalize the flow 

system in places (Kinnison and Sceva, 1963; Hansen and others, 

1994; Bauer and Hansen, 2000; Vaccaro and others, 2009). The 

compartmentalization limits the length of the flow paths, resulting 

in relatively short paths for such a large aquifer system. Structural 

control is exerted primarily by the major ridges in the Yakima 

Fold Belt (fig. 1; Hansen and others, 1994; Reidel and others, 

2002; Jones and others, 2006).

Groundwater flow in the Grande Ronde unit is 

compartmentalized, but not to the same extent as in the Saddle 

Mountains or Wanapum units. The large spatial extent of the 

Grande Ronde unit results in a large flow system with more 

interconnections than in the other two CRBG units. Where 

the unit outcrops, the water-level contours mimic landsurface topography and become a more subdued expression 

of topography as the unit becomes buried. In the more deeply 

buried parts of the Grande Ronde unit, the contours are smoother 

than contours for the other CRBG units. Similarly, water-level 

contours near geologic structures in the eastern part of the area 

are more subdued and smoother. The flow system in the Grande 

Ronde unit is controlled by the regional discharge locations along 

the Columbia River and major tributaries; that is, the hydraulic 

gradient in the Grande Ronde unit is toward the major streams 

(fig. 6). There may be a separate regional flow system in the 

deeper part of the unit, but data are insufficient to verify the 

presence of such a system.
Description of Groundwater-Flow System 13

tac13-0856_fig 05

Upper Colonnade

Platy Fracturing

Vesicle Zone

Vesicle Sheets

and Cylinders

Fanning Columns

FLOW TOP ZONE

Vesicular to rubbly and (or)

brecciated basalt

Pillow-palagonite complex,

vesicular base, rubble, or breccia

ENTABLATURE

FLOW INTERIOR

COLONNADE

FLOW BOTTOM ZONE

Dominant directions of

groundwater movement

EXPLANATION

Figure 5. Features within a typical Columbia River Basalt Group flow. Modified from Vaccaro (1986) and Reidel and others (2002). 
14 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

tac13-0856_fig 06

Grand Coulee Dam

B L U E

M O U N T A I N S

CASCADE RANGE

OKANOGAN HIGHLANDS

CASCADE RANGE

WALLOWA

MOUNTAINS

R O C K Y M O U N T A I N S 

R O C K Y M O U N T A I N S

HORSE HEAVEN HILLS

BEEZLEY HILLS

SNOHOMISH

SKAGIT

KING

PIERCE

LEWIS

CHELAN

OKANOGAN

KITTITAS

YAKIMA

SKAMANIA

KLICKITAT

HOOD

RIVER

CLACKAMAS

MARION

LINN

FERRY

STEVENS

PEND

ORIELLE

BENTON

GRANT

LINCOLN

SPOKANE

ADAMS

FRANKLIN

WALLA WALLA

WHITMAN

COLUMBIA

GARFIELD

ASOTIN

SHERMAN

JEFFERSON

CROOK

GILLIAM

MORROW

UMATILLA

WHEELER

MALHEUR

BAKER

WALLOWA

BONNER

KOOTENAI

BENEWAH

LATAH

LINCOLN

SANDERS

CLEARWATER

NEZ PERCE

WASHINGTON VALLEY

ADAMS

IDAHO

DOUGLAS

WASCO

GRANT

UNION

SHOSHONE

LEWIS

Sna k e

Ri v e r

Cr ab Cre e k

Col umbi a Ri v e r

Col umbi a

Nac he s R i v e r

U m a t i l l a Ri v e r

Col um

bia

Ri v e r

Ri v e r No rt h

For k John

Sna k e R i v e r

Imn a h a Rive r

Banks

Lake

Ri v e r

S po ka ne

Coeur

D'Alene

Lake

Moses

Lake

Grande Ro n d e Ri ver

Salmo n Ri v e r

J o h n D a y Ri v e r

M os es Co ul ee

Wa l l a Walla River

D e s c hu t e s

Roc k y

F o rd

Cr.

Kl i c k i t at River

Cl e a rwate r Ri v e r

Hangma n Cre e k

F. D . R o o s e v e l t La k e Palouse R i v er

S a t u s Creek

To p p e ni s h

Yakima River

Ellensburg

Umatilla

Walla Walla

Yakima

Heppner

MiltonFreewater

Pasco

Lewiston

Moscow

Spokane

The

Dalles

LaGrande

WASHINGTON

OREGON IDAHO

Location of

study area Estimated groundwater levels

Columbia Plateau Regional

 Aquifer System

Extent of Grande Ronde unit

Insufficient data

Generalized groundwater-

 flow direction

Lowest Highest

EXPLANATION

122° 120° 118° 116°

46°

48°

47°

45°

Base map modified from U.S. Geological Survey 

digital data, 1:2,000,000, 1972

0 50 100 MILES

0 50 100 KILOMETERS

Figure 6. Estimated groundwater levels and directions of lateral groundwater movement for the Grande Ronde unit, Columbia 

Plateau Regional Aquifer System, Idaho, Oregon, and Washington (from Snyder and Haynes, 2010).
Numerical Simulation of Groundwater Flow 15

Numerical Simulation of 

Groundwater Flow

Development of a calibrated groundwater-flow model 

allows for simulation and analysis of the movement of 

water through the hydrogeologic units that constitute the 

CPRAS. The USGS modular three-dimensional finitedifference groundwater-flow model, MODFLOW-NWT 

(Niswonger and others, 2011), a standalone program that is 

intended for solving problems involving drying and rewetting 

nonlinearities of the unconfined groundwater-flow equation, 

was used to simulate groundwater flow in the overburden 

deposits, interbeds, and basalt units of the aquifer system, 

and the interaction of the groundwater-flow system with 

surface water.

Modeling Strategy

The groundwater-flow model was developed to improve 

understanding of the hydrogeologic system. The model 

was formulated to test processes or aspects of the geology 

that previously have been identified or hypothesized to be 

important for controlling the flow and storage of groundwater 

in parts of the CPRAS (see section, “Previous Investigations”). 

In particular, the model was designed to investigate:

• Quantification of groundwater pumping;

• Flow through wells that are open to multiple aquifers;

• Enhanced groundwater recharge due to anthropogenic 

activities, especially large-scale irrigation projects;

• Horizontal barriers to groundwater flow corresponding 

to geologic structure;

• The large number of barriers to vertical groundwater 

flow corresponding to hundreds of dense lava flow 

interiors; and

• The connection between the groundwater and surfacewater systems.

The goal of the model is to help researchers understand 

the role of each of these mechanisms in producing observed 

long-term (decadal) changes in groundwater (Snyder and 

Haynes, 2010; Burns and others, 2012b) and changes to base 

flow in rivers and streams.

Because there is strong correlation of hydraulic 

measurements over tens of miles (Burns and others, 2012b) 

and steep vertical hydraulic gradients are present in the 

CPRAS, a coarse horizontal grid and fine vertical grid was 

created. Representation of appropriate horizontal hydraulic 

connectivity was achieved by designing the grid to follow 

the general pattern of hydrogeologic deposits by using 

the geologic model of Burns and others (2011), and using 

vertically thin cells to approximate typical basalt flow 

thicknesses of about 100 ft (Burns and others, 2011, 2012a). 

Using thin cells allowed the horizontal permeability of each 

layer to represent the transmissivity of individual CRBG 

aquifers, and the vertical permeability to represent individual 

basalt flow interior confining units. Although individual 

basalt flows and aquifers are not mapped to the level of 

detail necessary to define aquifers, using thin model layers 

ensured that shallow aquifers were appropriately connected 

to local surface-water features, and that deep aquifers were 

appropriately disconnected. The model grid was supplemented 

with a limited number of features representing various internal 

flow boundaries such as discrete barriers to horizontal flow, 

which generally were correlated to geologic structure, and 

rivers, streams, and wells.

Upstream Weighting Package

Because the system is being represented as a large 

number of highly anisotropic layers that may drain at various 

altitudes, numerical instabilities may occur as model cells go 

dry or rewet. Historically, groundwater-flow models of the 

CPRAS have assumed confined conditions for all layers to aid 

model stability; however, Burns and others (2012a) showed 

that in order to represent typical long-term declines in CRBG 

aquifers, it may be necessary to allow the drainage of regions 

of the model domain. To gain the necessary stability for the 

mathematical solution of the groundwater-flow equation, the 

recently developed Upstream Weighting (UPW) Package 

(Niswonger and others, 2011) was used.

The UPW Package is an internal flow package for 

MODFLOW-2005 that is intended to be used with the Newton 

Solver (NWT) Package for problems involving drying and 

rewetting nonlinearities of the unconfined groundwater-flow 

equation. The UPW Package treats nonlinearities of cell drying 

and rewetting by use of a continuous function rather than the 

discrete approach of drying and rewetting that is used in the 

Block-Centered Flow (BCF), Layer Property Flow (LPF), and 

Hydrogeologic-Unit Flow (HUF) Packages (Anderman and 

Hill, 2000; Harbaugh, 2005). This further enables application 

of the Newton solution method for unconfined groundwaterflow problems because conductance derivatives required 

by the Newton solution method are smooth over the full 

range of head for a model cell. A complete description of the 

UPW Package is available in Niswonger and others (2011).

Newton Solver

The Newton solution method is commonly used in the 

Earth sciences to solve nonlinear equations. The method is 

advantageous because many recently developed MODFLOW 

packages apply nonlinear boundary conditions. Additionally, 

the Newton solution method has been particularly beneficial 

in solving flow equations representing unconfined aquifers. 

During the course of this study, the Newton solution 

method provided greater model stability and improved 

model convergence compared to the solvers used in the 

standard MODFLOW-2005 code. A complete description 

of the Newton formulation is available in Niswonger and 

others (2011).
16 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Layering and Spatial Discretization

MODFLOW-NWT uses data sets that describe internal 

characteristics and boundary conditions of the groundwaterflow system, and calculates hydraulic heads and flow at 

discrete points (cell centers) within the model domain. 

The program requires that the groundwater-flow system be 

subdivided, vertically and horizontally, into rectilinear blocks 

called cells. The CPRAS model domain was subdivided 

into a horizontal grid of 126 rows and 131 columns using 

a uniform 3-km (9,842.5-ft) grid (fig. 7) that aligns with 

the 1-km SOWAT model grid (Kahle and others, 2011) 

allowing straightforward extraction of model boundary 

conditions (output from nine SOWAT cells were averaged 

for input to one MODFLOW cell). The large cell size 

and uniform grid spacing were selected as an appropriate 

trade-off between computational limitations, the geometry 

of the flow system, and the areal scale of observed regional 

head and flow patterns, which typically extend for tens to 

hundreds of kilometers (Snyder and Haynes, 2010, Burns and 

others, 2012b).

Vertically, the model domain was subdivided into the 

six geologic model units of Burns and others (2011): Grande 

Ronde Basalt, Vantage Interbed, Wanapum Basalt, Mabton 

Interbed, Saddle Mountains Basalt, and Overburden (table 2). 

Although individual aquifers are not reliably mapped at the 

regional scale, each of these basalt model units was discretized 

vertically into layers that are approximately 100-ft thick 

for most of the model domain, roughly approximating a 

typical CRBG lava flow thickness (Burns and others, 2012a). 

Simulation using 100 ft layers preserved appropriate localscale connectivity of individual basalt aquifers with hydrologic 

features. Historical and anecdotal evidence suggests that 

vertical hydraulic head gradients are orders of magnitude 

larger than horizontal gradients in basalt aquifers (Vaccaro 

and others, 2009; Burns and others, 2012a, 2012b), indicating 

that a relatively fine vertical discretization is appropriate for 

representing system dynamics.

The groundwater-flow model included 100 model layers, 

representing up to 17,000 ft (at the thickest part of the active 

flow system). Although the target thickness of individual 

basalt layers was 100 ft, model cell thickness ranged from 

10 to 200 ft for the Saddle Mountains and Wanapum units 

and the upper 30 layers of the Grande Ronde unit. Where 

necessary to represent the Grande Ronde with more than 

30 layers (total thickness greater than 3,000 ft), the remaining 

thickness was equally divided so that the total number of 

model layers was 100. These thicker cells were limited to very 

deep parts of the model domain, far below the parts of the 

system currently being pumped for groundwater. The thickest 

cells were 365 ft.

The basalt units were assigned to a discrete range 

of layers (based on total basalt unit thickness); no layer 

represented more than one CRBG unit, and the sedimentary 

Overburden and Interbed units were simulated in all remaining 

layers (table 2). Interbeds may be thin, but they were 

simulated across the entire extent of the associated basalt 

unit. Where a CRBG unit is not present, the hydrogeologic 

unit properties assigned to the model layers assigned to the 

missing CRBG unit were based on the overlying surface: 

deposits directly overlain by Wanapum Basalt were defined 

as Vantage Interbed, deposits directly overlain by Saddle 

Mountains Basalt were defined as Mabton Interbed, and 

deposits directly overlain by the land surface were defined as 

Overburden. The thickness of each interbed and Overburden 

cell was determined by dividing the total thickness of the 

respective unit into equal parts. The thinnest model cells occur 

under the unusual circumstances where Saddle Mountains 

Basalt overlies Grande Ronde Basalt with no Wanapum Basalt 

between. Because the overlying basalt is Saddle Mountains, 

the cells are simulated as Mabton Interbed. These cells occur 

infrequently (16 locations throughout the model domain), with 

the thinnest cell being 0.16 ft. The extents and thickness of 

the Overburden and CRBG units in the model are shown in 

figures 8A–D.

When generating the flow simulation model grid for the 

CRBG, the lowest layer was assumed to be the most extensive 

(consistent with the general trend of less voluminous lava 

flows as time progressed over geologic time [Lindholm and 

Vaccaro, 1988]), and all remaining layers were built upon this 

lowest layer until the total thickness of the basalt unit was 

represented. After the total thickness of basalt was assigned to 

model cells, sedimentary deposits were assigned to remaining 

overlying cells as described in the previous paragraph. The 

resulting groundwater-flow model preserved the possibility 

that low-lying sediments are in hydraulic connection 

with discontinuous basalt layers in both valley bottoms 

and at basalt-lava-flow margins, providing a reasonable 

representation of the interaction between CRBG aquifers and 

sedimentary units, which is consistent with the approach by 

Burns and others (2011).

The combination of lateral and vertical discretization 

resulted in 1,650,000 cells within the model grid, of which 

423,394 cells are active. The active cells include an area 

of almost 33,000 mi2 and are composed of 32,000 mi3 of 

aquifer-system material. Total model domain thickness ranged 

from 100 to 17,200 ft. Model layers 1–75 were simulated as 

convertible (confined/unconfined), and layers 76–100 were 

simulated as confined.
Numerical Simulation of Groundwater Flow 17

tac13-0856_fig 07

Major row or column

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

Pasco

EXPLANATION

Active model domain

Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

Inactive

Active model cells

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

John Day

River

John Day

River Deschutes River Deschutes River

Palouse RiverPalouse River

Yakima RiverYakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Active

model

domain

WASHINGTON

OREGON

IDAHO

Figure 7. Location and extent of the groundwater-flow model of the Columbia Plateau Regional Aquifer System, Idaho, 

Oregon, and Washington.
18 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Table 2. Correlation between geologic model units and groundwater-flow model layers of the Columbia Plateau Regional Aquifer 

System, Idaho, Oregon, and Washington.

[Geologic model units: Rock units are Saddle Mountains Basalt, Wanapum Basalt, Grande Ronde Basalt, and Older Bedrock. ft, foot]

Geologic model units 

(Burns and others, 2011)

Groundwater-flow 

model layers Description

Overburden 1–35 Approximately 100 ft model layers between the uppermost rock unit and land surface.

Saddle Mountains Basalt 11–18 Approximately 100 ft layers between the geologic model top of Mabton Interbed and the

geologic model top of Saddle Mountains Basalt. Cells are constructed from the bottom up, 

with the lowest layer being the most laterally extensive.

Mabton Interbed 19–35 Equally spaced layers between the geologic model bottom of Saddle Mountains Basalt and 

the next rock unit below.

Wanapum Basalt 20–34 Approximately 100 ft layers between the geologic model top of Vantage Interbed and the

geologic model top of Wanapum Basalt. Cells are constructed from the bottom up, with the

lowest layer being the most laterally extensive.

Vantage Interbed 35 Layer between the geologic model bottom of Wanapum Basalt and the next 

rock unit below.

Grande Ronde Basalt 36–100 Approximately 100 ft layers between the geologic model top of Older Bedrock and the 

geologic model top of Grande Ronde Basalt. Cells were constructed from the top down, 

using a trend surface for top of Grande Ronde Basalt as a guide surface, allowing the 

representation of river and stream incision exceeding model cell thickness.

Older Bedrock No flow
Numerical Simulation of Groundwater Flow 19

tac13-0856_fig 08a

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

EXPLANATION

Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

EXPLANATION

Active model domain

Overburden unit extent and thickness, in feet

 10.8 – 100.0

100.1 – 200.0

200.1 – 300.0

300.1 – 400.0

400.1 – 500.0

500.1 – 600.0

600.1 – 700.0

 700.1 – 800.0

 800.1 – 900.0

 900.1 – 1,000.0

1,000.1 – 1,100.0

1,100.1 – 1,200.0

1,200.1 – 1,300.0

1,300.1 – 1,359.9

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

A

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 8. Extent and thickness of the (A) Overburden unit, (B) combined Saddle Mountains basalt and Mabton interbed unit,

(C) combined Wanapum basalt and Vantage interbed unit, and (D) Grande Ronde basalt unit, Columbia Plateau Regional Aquifer 

System, Idaho, Oregon, and Washington.
20 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

tac13-0856_fig 08b

13.3 – 100.0

100.1 – 200.0

200.1 – 300.0

300.1 – 400.0

400.1 – 500.0

500.1 – 600.0

600.1 – 700.0

700.1 – 800.0

800.1 – 900.0

900.1 – 1,000.0

1,000.1 – 1,100.0

1,100.1 – 1,200.0

1,200.1 – 1,300.0

1,300.1 – 1,400.0

1,400.1 – 1,500.0

1,500.1 – 1,600.0

1,600.1 – 1,700.0

1,700.1 – 1,745.1

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

EXPLANATION Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

EXPLANATION

Active model domain

Saddle Mountains unit extent and thickness, in feet

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

B

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 8.—Continued
Numerical Simulation of Groundwater Flow 21

tac13-0856_fig 08c

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

EXPLANATION Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

Active model domain

Wanapum unit extent and thickness, in feet

10.2 – 100.0

100.1 – 200.0

200.1 – 300.0

300.1 – 400.0

400.1 – 500.0

500.1 – 600.0

600.1 – 700.0

700.1 – 800.0

800.1 – 900.0

900.1 – 1,000.0

1,000.1 – 1,100.0

1,100.1 – 1,200.0

1,200.1 – 1,300.0

1,300.1 – 1,400.0

1,400.1 – 1,500.0

1,500.1 – 1,600.0

1,600.1 – 1,700.0

1,700.1 – 1,800.0

1,800.1 – 1,900.0

1,900.1 – 2,000.0

2,000.1 – 2,100.0

2,100.1 – 2,200.0

2,200.1 – 2,300.0

2,300.1 – 2,400.0

2,400.1 – 2,500.0

2,500.1 – 2,600.0

2,600.1 – 2,700.0

2,700.1 – 2,800.0

2,800.1 – 2,900.0

2,900.1 – 2,995.0

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

C

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 8.—Continued
22 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

tac13-0856_fig 08d

10.9 – 500.0

500.1 – 1,000.0

1,000.1 – 1,500.0

1,500.1 – 2,000.0

2,000.1 – 2,500.0

2,500.1 – 3,000.0

3,000.1 – 3,500.0

3,500.1 – 4,000.0

4,000.1 – 4,500.0

4,500.1 – 5,000.0

5,000.1 – 5,500.0

5,500.1 – 6,000.0

6,000.1 – 6,500.0

6,500.1 – 7,000.0

7,000.1 – 7,500.0

7,500.1 – 8,000.0

8,000.1 – 8,500.0

8,500.1 – 9,000.0

9,000.1 – 9,500.0

9,500.1 – 1,0000.0

1,0000.1 – 10,500.0

10,500.1 – 11,000.0

11,000.1 – 11,500.0

11,500.1 – 12,000.0

12,000.1 – 12,500.0

12,500.1 – 13,000.0

13,000.1 – 13,500.0

13,500.1 – 14,000.0

14,000.1 – 14,500.0

14,500.1 – 15,000.0

15,000.1 – 15,500.0

15,500.1 – 15,761.8

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

EXPLANATION

Active model domain

Grande Ronde unit extent and thickness, in feet

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

D

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 8.—Continued
Numerical Simulation of Groundwater Flow 23

Temporal Discretization

Two models were constructed to simulate groundwater 

flow in the CPRAS: (1) a steady-state predevelopment model 

representing conditions before large-scale pumping and 

irrigation altered the system, and (2) a transient simulation 

representing the period 1900–2007. Agricultural irrigation, 

from surface-water and groundwater sources, occurred 

throughout the Columbia Plateau from the late 1800s 

through early 1900s. However, large-scale groundwater 

level changes were not recorded until after the 1950s, when 

large-scale irrigation projects and groundwater development 

greatly increased. The steady-state predevelopment model 

was constructed to represent important aspects of the flow 

system without the complexities of large-scale groundwater 

withdrawals and aquifer storage properties.

Time discretization in MODFLOW has two levels of 

division: stress periods and time steps. External stresses, 

such as pumping rates and recharge, change at the beginning 

of each stress period. Stress periods, which are commonly 

monthly, quarterly, or annual in regional models, are divided 

into time steps. The groundwater flow equation is solved 

(meaning heads and flows are calculated) for each time 

step. The first stress period of the MODFLOW simulation 

represented steady-state predevelopment conditions. The 

steady-state predevelopment model provided initial conditions 

for the transient model representing early development 

(1900–1919), followed by annual stress periods that simulated 

the groundwater development that occurred throughout most 

of the 20th century and early 21st century (1920–2007), for 

a total of 90 stress periods. Recharge, pumping, and crossconnection of aquifers through open well boreholes (hereafter 

called commingling wells) were estimated for each stress 

period (see section, “Hydrologic Boundaries”).

The long simulation using an annual stress period 

allowed for a temporal assessment that accounted for large 

changes in pumping and a range in climatic conditions and, 

thus, large ranges in simulated groundwater declines. Because 

the goal of this modeling effort was to understand the decadalscale trends in water levels and base flow conditions, seasonal 

effects of pumping, and recharge were not considered with the 

current model.

Hydraulic Properties

Simulated flow through the model domain is controlled 

by the hydrologic boundaries (where water enters and leaves 

the domain) and the hydraulic properties in the flow system 

(that control the amount of water that may be stored and the 

ease with which water moves through the system). Each of 

the geologic units was represented as locally homogeneous 

and vertically anisotropic with bulk horizontal hydraulic 

conductivity generally being significantly higher than vertical 

conductivity. Discrete barriers to groundwater flow (Burns 

and others, 2012b) were simulated as linear horizontal flow 

barriers with much lower conductance than surrounding 

geologic units. 

Horizontal Hydraulic Conductivity

Horizontal hydraulic conductivity (Kh) controls the ease 

with which water moves horizontally through geologic units. 

For the 100-ft basalt layers, Kh was a bulk conductivity that 

represents the relatively thin aquifer that typically occupies 

approximately 10 percent of the total thickness (Burns and 

others, 2012a). The sedimentary interbeds generally were 

thin compared to the total basalt thickness at each location; 

therefore, a single value of Kh for each interbed in each 

region of the model was appropriate for a regional scale-flowsimulation model. Thick Overburden units were separated 

horizontally by CRBG unit structural highs (fig. 3), and 

these units may exhibit significant amounts of heterogeneity 

horizontally and vertically (Ely and others, 2011). However, 

there are few large-volume sedimentary basins, and from a 

regional perspective, these basins are small compared to the 

regional CRBG aquifer system. Because the primary purpose 

of this groundwater-flow modeling effort is to understand 

the flow in the basalt units, the Overburden unit in a basin or 

region was simulated as locally homogeneous (a single value 

of Kh), although Kh and other hydraulic properties were 

allowed to vary between basins during parameter estimation 

and uncertainty analysis.

Initial values of Kh were assigned based on values 

tabulated from numerous previous studies and analysis of 

specific-capacity/aquifer-test data (Kahle and others, 2011, 

table 3), and on spatial distributions of these hydraulic 

properties in the calibrated CP-RASA model (Hansen and 

others, 1994). The initial estimates of Kh for the basalt 

interflow zones and flow interiors decreased with depth 

based on Weiss (1982) and Hansen and others (1994). The 

assumption is that overburden pressure and secondary 

mineralization have reduced pore space with depth, over time. 

Hansen and others (1994) decreased the Kh values with depth 

as a parabolic expression; resulting in a 40-percent decrease 

of Kh from the surface to a depth of 3,000 ft. The method 

provided reasonable estimates of Kh for the CP-RASA model 

(Hansen and others, 1994); therefore, it was used as starting 

values for the current model. The reduction of Kh with depth 

has been used by others in modeling groundwater flow in deep 

systems (for example, the Death Valley regional groundwaterflow system [Faunt and others, 2010] and the California 

Central Valley [Faunt and others, 2009]). The role of Kh 

decrease with depth for the CPRAS was examined during 

model calibration.
24 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Vertical Hydraulic Conductivity

Vertical hydraulic conductivity (Kv) controls the ease 

with which water moves vertically through geologic units. 

For the 100-ft basalt layers, Kv is a bulk conductivity that 

represents the sequential vertical transmission of water 

through a heterogeneous basalt flow (both interflow and flow 

interior). Because the flow interior is thick and generally 

impermeable, basalt Kv is dominated by this part of the basalt 

flow. Similar to Kh, Kv for the sedimentary interbeds and 

Overburden was assumed locally homogeneous although 

variable between basins. Initial values of Kv were assigned 

based on values tabulated from previous studies and analysis 

of specific-capacity/aquifer-test data (Kahle and others, 

2011, table 4).

Storage Properties

Two storage parameters represent two different storage 

mechanisms: specific storage (Ss) represents the change in 

storage associated with compressibility of water and the 

geologic material in confined aquifers, and specific yield (Sy) 

represents the water that may be gravity drained from the pore 

space within the geologic material in unconfined aquifers. 

In confined aquifers, which are completely filled with water, 

changes in the amount of water stored are due to compression 

or expansion of the water and the aquifer framework. Because 

water and rock are not very compressible, specific storage 

values typically are quite small (approximately 10-4 ft-1). In 

unconfined aquifers, changes in the amount of water stored are 

due to actual draining or filling of the open space within the 

aquifer, so values generally are large (approximately 10-1 ft-1). 

Both confined and unconfined conditions occur in the CPRAS.

Storage properties of the aquifers are highly variable, and 

a general lack of information makes reliable areal estimates 

difficult to obtain. A wide range of values for storage terms 

(Kahle and others, 2011, table 5), which typically reflect Sy 

when values are greater than 0.1, Ss when values are less 

than 5.0 × 10-4 ft-1, and a mixed model-estimated calibration 

storage coefficient when values are between Sy and Ss. When 

considering estimated storage coefficients, an additional 

complication is the volume being analyzed. For example, the 

porosity of a flow interior may be significantly smaller than 

the porosity of an interflow zone, although each zone may be 

represented in a single model cell. For the current analysis 

of the CRBG system, storage coefficients were required that 

represent the average storage of the entire vertical profile of 

the basalt flow (both interflow and interior).

Neutron porosity estimates from boreholes on the 

Hanford Site were used to estimate a bulk specific yield of 

0.1 (dimensionless). Values of Ss for basalt were estimated as 

2.5 × 10-5 ft-1, which were consistent with the reported values 

of Kahle and others (2011), and were used as the initial values in 

the model. Model stability and convergence were not sensitive 

to Ss and Sy values.

Geologic Structures

The CRBG effective Kh is decreased by several orders of 

magnitude in areas of intense folding and faulting, especially 

fault-associated anticlines (Hansen and others, 1994; Packard 

and others, 1996; Reidel and others, 2002, Ely and others, 

2011). Where the distribution of water-level measurements 

are closely spaced, high Kh regions often can be separated by 

narrow transition zones (commonly coincident with a fault 

or fold trace) that are significantly less permeable (Burns and 

others, 2012b). Low permeability along geologic structures 

may be due to the offsetting of interflow zones through faulting, 

which juxtaposes thin basalt aquifers with thick confining units 

and produces low-conductivity fault breccia and gouge material 

at that interface (as described by Stearns [1942] and Newcomb 

[1965, 1969]), and the decreasing pore space through deposition 

of secondary minerals. Similarly, intense folds accommodate 

slip along planes of low strength (interflows), grinding and 

thinning deposits in the slip plane, potentially reducing 

permeability. Not all folds and faults act as flow barriers and 

currently, there is no predictive relation between fold and 

fault type and the extent that the structure forms a barrier to 

groundwater flow.

Displacement of individual basalt flows along faults, 

however, locally can enhance vertical movement of water by 

providing fractured zones across basalt flows that could serve 

as conduits for vertical groundwater flow. Preferential pathways 

may be developed locally, but faults commonly tend to act 

as barriers to flow regionally. Because CRBG flow interiors 

chemically weather to clays when broken rock (fault gouge) is 

exposed to moving water, faults are believed to be “self-sealing” 

on the geologic timescale, which explains why most fault zones 

across the CPRAS are not vertical conduits for groundwater 

flow.

The Horizontal-Flow Barrier (HFB) package of Hsieh 

and Freckleton (1993) was used to simulate low-permeability 

geologic structures. Additionally, HFBs were used at all high 

offset faults identified in the geologic model of Burns and others 

(2011). Flow barriers comprising 34,465 model cells from 

model layers 10–100 were grouped into 25 hydraulic categories 

(fig. 9) and initial hydraulic characteristic values of the 25 

hydraulic categories ranged from 9.95×10-16 to 1.17×10-8 1/d, 

with the smallest value having the largest control. The hydraulic 

characteristic of the horizontal-flow barrier for unconfined 

layers is the barrier hydraulic conductivity divided by the width 

of the horizontal-flow barrier along the flow path. The hydraulic 

characteristic of the HFB for confined layers is the barrier 

transmissivity divided by the barrier width.
Numerical Simulation of Groundwater Flow 25

tac13-0856_fig 09

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

41

10

21

22

23

32

33

34

51

61

62

63

64

65

66

67

70

71

68

69

52

53

54

35

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

EXPLANATION

Horizontal-flow barrier, by group

1

10

21

22

23

32

33

34

35

41

51

52

53

54

61

62

63

64

65

66

67

68

69

70

71

Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 9. Location and hydraulic characteristic group of the groundwater model structure (horizontal-flow barriers), Columbia 

Plateau Regional Aquifer System, Idaho, Oregon, and Washington.
26 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Model fit was considered reasonable with no zones 

of preferential vertical flow associated with faults, which 

supported the assumption that most vertical preferential flow 

associated with faulting has sealed over geologic time with 

respect to the regional groundwater-flow system. However, 

some vertical flow may be occurring at a local scale. For 

regions with a high density of mapped geologic folds and 

faults and steep hydraulic gradients, but where the distribution 

of water-level measurements was not dense enough to identify 

the most important structures for controlling flow, a lower Kh 

was used to represent the effective permeability of the region.

Hydrologic Boundaries

Boundary conditions define the locations and manner 

in which water enters and exits the active model domain. 

The specified boundaries of the model coincide as much as 

possible with natural hydrologic boundaries. Three types of 

model boundaries were used in the CPRAS model: (1) no-flow 

boundaries (groundwater divides and basalt unit extents), 

(2) head-dependent flux boundaries (drains, rivers, and 

multi-node wells), and (3) specified-flux boundaries (recharge 

and pumping). Head-dependent and specified-flux boundary 

conditions were developed for predevelopment conditions and 

for 1900–2007.

Evapotranspiration (ET) was not simulated using the 

groundwater-flow model. Instead, net recharge (precipitation 

plus irrigation minus evapotranspiration minus runoff) was 

computed using the SOWAT model (Kahle and others, 2011), 

and used as input to the groundwater-flow model. This 

approach accounted for most ET, with only a small part of the 

water budget neglected by ignoring groundwater-fed wetlands.

Hydraulic Conditions Along the Periphery of the 

Model

All lateral model boundaries were assumed to be no-flow 

boundaries (fig. 10), with the understanding that some water 

might be entering the model domain, but the amount likely is 

small compared with recharge. During model calibration, this 

assumption was determined to be reasonable, although during 

transient simulation, some wells along the periphery could 

not support the desired pumping, indicating that (1) wells 

were extracting most of their water from beneath the CRBG, 

(2) local stream capture was occurring, but was not simulated 

by the model, or (3) some lateral inflow of water exists.

The model domain was selected by identifying hydraulic 

boundaries with low (ideally zero) groundwater flow across 

the boundary into or out of the model domain. These hydraulic 

divides commonly are coincident with major topographic 

divides associated with mountains and rivers. For example, 

hydraulic evidence shows that both surface water and 

groundwater flows away from the crest of the Blue Mountains 

(fig. 6). To the north, the Columbia River is deeply incised 

into older pre-Miocene rock, and was represented as a no-flow 

boundary because of extremely low permeability of preMiocene rock. However, incised rivers at no-flow boundaries 

still may have river and drain cells simulated to allow water to 

enter or leave the domain. This can allow water to drain from 

the aquifer system into canyons, if appropriate.

Near the depositional margins of the CRBG, the basalt 

was deposited as intra-canyon flows, rather than large sheet 

flows, and lateral inflow of water from adjacent, hydraulically 

connected deposits might occur more easily in these locations 

than in areas of unconnected margins. However, recent 

simulation results for the Yakima River Basin (Ely and 

others, 2011) indicated that most water that is recharged in 

the Cascade Range beyond the extent of the CRBG deposits 

(fig. 1) returns to the stream rather than flowing laterally 

into the CRBG. This finding was assumed to be applicable 

to CRBG deposits in the southwestern boundary and in the 

foothills of western Idaho (fig. 10). The foothills of western 

Idaho are significantly different from the Cascade Range, but 

recharge potential also is much smaller in Idaho; therefore, the 

boundary was selected under the presumption that some water 

might be entering the CPRAS there.

The two remaining major areas of the CRBG boundary 

are where the Deschutes River crosses the model boundary in 

Oregon, and near Lewiston, Idaho (fig. 7). Near Lewiston, a 

fault offset of several thousand feet was interpreted as a likely 

barrier to flow. For the Deschutes River, the continued trend of 

Blue Mountains geologic structures was used to estimate the 

boundary. River streamgage data and Deschutes River Basin 

groundwater-flow model results (Gannett and Lite, 2004) 

support this as a reasonable model boundary.

The base of the CRBG was selected as the lower 

(no-flow) boundary for the model. Previous studies have 

used this boundary (Hansen and others, 1994) and the low 

permeability of the underlying deposits support the location 

of the no-flow boundary, especially near the center of the 

study area.

Groundwater Pumpage

The spatial distribution of pumping was estimated on a 

1-km grid coincident with the SOWAT grid for 1920–2007. 

SOWAT-derived estimates of irrigation pumping for 

1985–2007 (Kahle and others, 2011) were supplemented 

with previously published estimates of large-capacity well 

withdrawals (Cline and Collins, 1992) and new estimates for 

municipal, industrial, rural, and residential uses.
Numerical Simulation of Groundwater Flow 27

tac13-0856_fig 10

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

EXPLANATION Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

River (RIV) cell locations

Drain (DRN) cell locations

No-flow boundary

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Snake River

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 10. Location of the groundwater-flow model river and drain cells, and the no-flow boundary, Columbia Plateau Regional 

Aquifer System, Idaho, Oregon, and Washington.
28 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Amounts of historical groundwater pumping in the 

Columbia Plateau for the period 1945–84 were determined 

by Cline and Collins (1992) as part of the RASA study 

of the 1980s. Methods used to estimate well withdrawals 

are: (1) compiling electrical power consumption records, 

(2) multiplying irrigated acreage by a water-application rate, 

and (3) using other methods, such as data from published 

reports, field visits, population, water-use, water rights, and 

well construction logs, and information from telephone 

surveys. The resulting, gridded estimates could be considered 

the most detailed and accurate estimates of groundwater 

pumping and were used to hind cast and augment pumpage 

estimates for 1920–84.

For the period 1985–2007, estimates of total pumpage 

were made by adding estimates of irrigation pumpage 

extracted from the 1-km grid-spacing SOWAT model results 

of Kahle and others (2011) to estimates of pumpage for all 

other anthropogenic activities (fig. 11). For SOWAT grid 

cells occupied by irrigated agriculture, the SOWAT model 

uses estimates of actual evapotranspiration to calculate when 

irrigation should occur to meet plant needs. When irrigation 

water is required, SOWAT estimates the amount of water 

necessary to fill the soil and uses an irrigation inefficiency 

to estimate the total amount of water applied from each 

irrigation source. Pumpage was smoothed in the same manner 

as pumping-associated recharge (see section, “Groundwater 

Recharge”) so that model inflow and outflow estimates were 

made in a consistent manner. Pumpage also was estimated for 

each year for public-supply wells, private domestic wells, and 

other water uses requiring water-rights permits. If the well 

location was known precisely, pumpage was assigned to the 

associated SOWAT 1-km grid cell. If water use was assigned 

to a mapped area, then pumpage was assigned to the 1-km grid 

proportionally to the fraction of the mapped area intersecting 

each 1-km grid cell. For example, census blocks were used 

to estimate rural residential pumpage for census years based 

on typical per capita water use, and pumpage was distributed 

equally across the entire census block. Between census years, 

pumpage was assumed to linearly vary between the years. 

Pumpage estimates using this method were compared to 

independent estimates for the Umatilla River Basin (Kahle 

and others, 2011, fig. 21B), the Odessa subarea in east-central 

Washington (Cline, 1984), and the area of the Yakima River 

Basin that intersects the CPRAS model domain (Ely and 

others, 2011). All estimates indicated good agreement.

Pumpage for the period 1920–84 was constructed by 

hind casting the pumpage estimates for 1985–2007 using the 

tabulated pumpage estimates of Cline and Collins (1992). The 

Cline and Collins tabulated pumpage estimates for 1984 were 

compared to 1985 estimates made using the SOWAT 1-km 

grid methods to ensure that estimates using these dissimilar 

methods were consistent. To compare the two methods of 

tac13-0856_fig11 Groundwater pumpage, in acre-feet

0

200,000

400,000

600,000

800,000

1,000,000

1,200,000

1,400,000

1,600,000

1980 1985 1990 1995 2000 2005 2010

Water year

Total

Irrigation

Domestic

Industrial Public supply

and other

Estimated 1-kilometer-grid space SOil WATer model

Figure 11. Estimated groundwater pumpage, by 

category, Columbia Plateau Regional Aquifer System, 

Idaho, Oregon, and Washington, 1985–2007.

estimation, the 1-km grid was intersected with the quartertownship (9 mi2) grid of Cline and Collins (1992), and the 

total intersected pumpage was compared with tabulated 

values by Cline and Collins. These methods showed good 

agreement (fig. 12).

Estimating periods of historical pumpage was 

accomplished differently, depending on whether or not the 

1-km grid cell intersected a Cline and Collins (1992) grid 

cell with appreciable pumpage. If Cline and Collins did 

not estimate appreciable pumpage, the 1985 pumpage rate 

at each 1-km cell was tied to an existing nearby deep well 

(typically within a 3- by 3-km grid centered on the 1-km cell), 

and simulated pumping began when the presumed well was 

constructed. For Cline and Collins cells with non-negligible 

amounts of pumpage, the refined 1-km pattern of 1985 

pumping from SOWAT was assumed to be the correct pattern, 

and the magnitude of pumpage using this pattern was scaled 

by the estimates of Cline and Collins using linear interpolation 

for years between estimates, ensuring that both magnitude and 

distribution of pumpage matched 1985 1-km grid estimates 

exactly. Pumpage was assumed to be negligible prior to 

1920, because pumps with significant hydraulic lift were not 

common before then (Cline and Collins, 1992). The earliest 

estimates of pumpage by Cline and Collins were linearly 

decreased to zero in 1920, even though there are a few records 

for wells prior to this time. 
Numerical Simulation of Groundwater Flow 29

tac13-0856_fig12

Water year

1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000 2010

Groundwater pumpage, in acre-feet

0

200,000

400,000

600,000

800,000

1,000,000

1,200,000

1,400,000

EXPLANATION

Idaho (Cline and Collins, 1992)

Oregon (Cline and Collins, 1992)

Washington (Cline and Collins,1992)

Idaho (SOWAT 1-kilometer [km] grid)

Oregon (SOWAT 1-km grid)

Washington (SOWAT 1-km grid)

Idaho (intersection of grids)

Oregon (intersection of grids)

Washington (intersection of grids)

Idaho (simulated)

Oregon (simulated)

Washington (simulated)

Figure 12. Estimates of groundwater pumpage from Cline and Collins (1992) and SOil WATer (SOWAT) balance (Kahle 

and others (2011), and simulated groundwater-pumpage, Columbia Plateau Regional Aquifer System, Idaho, Oregon, and 

Washington.
30 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Groundwater pumping was simulated using the MultiNode Well (MNW2) Package of MODFLOW-NWT (Konikow 

and others, 2009). Initially, pumping at each 1-km grid cell 

was assigned to the deepest nearby well. For 1-km grid cells 

with pumping estimates derived from Cline and Collins (1992) 

estimates, the pumping well is the deepest nearby well during 

the stress period, so pumping moves from shallow wells to 

deeper wells over time. Because wells in the CRBG typically 

are open from the base of the Overburden to the bottom of 

the well, water can be supplied from any open aquifer to meet 

water demand. If well records indicated full penetration of the 

CRBG into the deeper pre-Miocene rocks (not simulated), the 

amount of pumping assigned to these wells was proportional 

to the fraction of the total length of the borehole that intersects 

model units, as estimated using the geologic model of Burns 

and others (2011).

All wells with construction records were simulated in 

the model if the well either potentially commingles aquifers 

(the length of the borehole spans multiple model layers) or 

has pumping assigned to it. Simulation of these wells allows 

the net pumping and net commingling influence of CPRAS 

wells to be represented over time. If pumping was indicated 

by the 1-km grid pumping estimates, but no well record was 

available within the MODFLOW model cell, pumping with no 

commingling was simulated from the uppermost aquifer.

Use of the MNW2 Package allowed representation of 

both pumping and commingling effects. Commingling wells 

were simulated using the SKIN option in MNW2 where 

the hydraulic conductivity of a cylindrical region near the 

borehole is assumed to be less than the hydraulic conductivity 

of the aquifer system (Konikow and others, 2009). The skin 

conductivity (Kskin) allows simulation of the imperfect 

connection between the aquifer system and the borehole that 

results from geologic and well-construction heterogeneity, 

and this parameter allows simulation and adjustment of the 

net effect of commingling. Wells that only intersect a single 

model layer were simulated using the THIEM option in the 

MNW2 Package.

Following preliminary transient simulations, it was 

observed that the desired amount of pumping was not being 

extracted from the model for current conditions (2000-2007). 

This occurred because Kskin values that resulted in reasonable 

time-varying patterns of drawdown were restricting inflow to 

wells sufficiently that wells were shutting off due to drawdown 

limits without producing enough water. This apparent 

contradiction occurred because pumping rates are controlled 

by the aquifers that are well-connected to wells (aquifers that 

can supply large amounts of water easily), but commingling 

between aquifers is limited by aquifers that are less wellconnected to wells (reduced flow occurs at restrictions). To 

split the pumping and commingling processes, pumping was 

reallocated to vertical strings of single node THIEM wells, and 

commingling was represented with multi-node wells with a 

specified Kskin and no associated pumping.

For each vertical string of single node wells, pumping 

was assigned to each node proportionally to the relative 

potential for each node to supply water. A node was defined as 

possibly contributing if a well is contained in the model cell 

and if simulated hydraulic head in the model cell was above 

the bottom of the model cell during the stress period. Relative 

potential to supply water was defined as the product of the 

hydraulic head above cell bottom (representing flow potential), 

the cell thickness (representing area available to flow), and an 

estimate of relative hydraulic conductivity (representing the 

relative ease of flow). Based on preliminary steady-state model 

calibration, hydraulic conductivity of Saddle Mountains basalt 

and Overburden was assumed to be twice the conductivity of 

the Wanapum basalt, which was twice the conductivity of the 

Grande Ronde basalt. The resulting distribution of pumping 

allows water to be removed from upper layers more easily 

than from lower layers, consistent with the typical history of 

well construction. Estimation of single node pumping wells 

required the estimation of hydraulic head; therefore, the 

transient model was first run allowing commingling, but not 

allowing pumping. Iterative application of the method allowed 

reallocation of simulated pumping based on successive 

improved estimates of head distribution over time.

The methodology to assign pumping and commingling 

wells resulted in a maximum of 20,752 commingling wells 

spanning from a minimum of 2 layers to a maximum of 

38 layers for a total of 90,938 nodes. Kskin was set to 0.1 ft/d 

for all commingling well cells. Pumping wells were assigned 

to 62,446 nodes. Final estimates and the spatial distribution of 

groundwater pumping for current conditions (2000–2007) are 

shown in figure 13.

Groundwater Recharge

Over most of the study area, natural recharge occurs 

only during the wet winters. During the dry summer months, 

recharge is dominated by irrigation inefficiency and associated 

irrigation infrastructure (for example, irrigation canals). 

Irrigated areas have higher antecedent soil moisture at the 

start of the wet season, increasing the fraction of precipitation 

that recharges groundwater. Groundwater recharge was 

estimated using three methods to accommodate three different 

periods in the simulation: (1) correlation of precipitation 

to recharge, (2) SOWAT estimates, and (3) hind casting of 

SOWAT estimates.
Numerical Simulation of Groundwater Flow 31

tac13-0856_fig 13

Mean annual groundwater pumpage,

in acre feet

 0 No

symbol

 0.1 – 25.0

 25.1 – 50.0

 50.1 – 100.0

 100.1 – 250.0

 250.1 – 500.0

 500.1 – 1,000.0

1,000.1 – 2,236.3

EXPLANATION

Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

Active model domain

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 13. Mean annual groundwater pumpage for current conditions (2000–2007), Columbia Plateau Regional Aquifer System, 

Idaho, Oregon, and Washington.
32 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Predevelopment groundwater recharge was estimated 

using the annual 1895–2007 gridded historical estimates of 

annual precipitation provided by the Parameter-elevation 

Regressions on Independent Slopes Model (PRISM) (PRISM 

Climate Group, 2004) and the Bauer-Vaccaro regression 

equation (Bauer and Vaccaro, 1990) relating annual 

precipitation to groundwater recharge. Recharge estimates for 

1920–84 were derived from the SOWAT recharge estimates 

combined with the groundwater pumping estimates of Cline 

and Collins (1992) and a timeline of construction of largescale surface-water irrigation projects. Recharge estimates 

for 1920 were used for the period 1900–19. Annual recharge 

estimates from the SOWAT modeling results of Kahle and 

others (2011) were limited to the period 1985–2007 because 

of the availability of evapotranspiration estimates from 

satellite imagery.

To account for the effect of decadal variations in natural 

recharge, the SOWAT-estimated recharge was divided into 

natural recharge and the additional recharge from irrigation. 

Additional recharge accounts for irrigation inefficiency and 

high antecedent soil-moisture content, and is computed as the 

difference between the total SOWAT-estimated recharge and 

the recharge estimated using the Bauer and Vaccaro regression 

equation. The additional recharge was divided into surfacewater and groundwater sources using the estimates of fraction 

of irrigated land in each 1-km grid cell supplied by each 

source used in the SOWAT model.

For the period 1985–2007, recharge estimates were 

extracted from the 1-km grid SOWAT modeling results of 

Kahle and others (2011). For grid cells dominated by irrigated 

agriculture, the SOWAT model used a simplified energy 

balance and soil reservoir model to estimate net recharge into 

the groundwater system from irrigation and precipitation. For 

all other grid cells, net recharge was estimated from annual 

precipitation using the Bauer and Vaccaro regression equation 

(Bauer and Vaccaro, 1990). Spatially-distributed annual 

estimates of precipitation were extracted from historical 

PRISM simulations (PRISM Climate Group, 2004).

Recharge for 1920–84 was reconstructed by hind casting 

the 1985–2007 recharge estimates. This method assumed that 

current irrigation is the most extensive and that irrigation was 

less extensive in the past. Irrigation from groundwater was 

related to gridded estimates of pumping (Cline and Collins, 

1992), and irrigation from surface water was tied to the history 

of large-scale irrigation projects.

Predevelopment recharge was smaller than recharge 

under existing conditions due to the infiltration of irrigation 

water during delivery and application to croplands. Irrigation 

projects exist throughout much of the Columbia Plateau, but 

the largest areas of intense irrigation are in the Yakima and 

Walla Walla River Basins, Umatilla subarea, and Columbia 

Basin Irrigation Project (fig. 2). Detailed annual estimates 

of historical surface-water irrigated acreage do not exist, 

so irrigation, and therefore recharge from irrigation, was 

increased linearly from the first year of irrigation water 

delivery until full irrigation deliveries were in effect.

To replicate the general timing and amount of recharge 

from large surface-water irrigation projects, recharge was 

systematically increased to 1985 SOWAT-estimated values 

using some simple rules: (1) CBIP irrigation was assumed to 

start in 1950 with irrigation increasing linearly to 100 percent 

in 1953; (2) Yakima River Basin surface-water irrigation was 

assumed to start in 1921 with irrigation increasing linearly to 

100 percent in 1930; and (3) Umatilla subarea surface-water 

irrigation was assumed to be 24 percent at the start of the 

transient simulation (1900), increasing linearly to 34 percent 

in 1961, to 75 percent in 1975, and to 100 percent by 1984 

(M. Ladd, Oregon Water Resources Department, written 

commun., 2012). Surface-water irrigation outside of the large 

projects in the Yakima River Basin, Umatilla subarea, and 

the CBIP are a relatively small fraction of total surface-water 

irrigation; these irrigated lands were assumed to exist for the 

entire period of transient simulation. Implementation of the 

previous assumptions resulted in recharge stresses that are 

considerably simpler than reality, but provided a reasonable 

timing of stresses to ensure the timing of changes in recharge 

correspond to the approximate history of water use on the 

Columbia Plateau.

The amount of estimated additional recharge from 

irrigation was examined and two anomalies were noted. First, 

in areas of sparse irrigation, where 1-km model cells were 

not dominated by irrigated agriculture, extra recharge was 

proportionally too high. This was attributed to the fact that 

the soil-water balance model assumed the entire cell was 

irrigated, leaving antecedent soil moisture too high for parts 

of the cell that are not irrigated, overestimating recharge 

from winter precipitation. Because SOWAT estimates actual 

ET, the rate of irrigation water applied and recharged from 

irrigation inefficiency should be correct. The total area 

occupied by these cells was relatively small compared to the 

entire model domain, so the error in total estimated recharge 

was negligible given the model scale. Second, recharge 

from irrigation varied in an oscillatory manner from year to 

year by about 20 percent, possibly caused by the method of 

representing when irrigation occurred. The soil-water balance 

model assumes a deep soil zone, resulting in infrequent 

addition of irrigation water (typically 4–5 times a year). 

Because SOWAT assumes that water is added to the deep soil 

zone when “maximum allowable soil depletion” is reached, 

it is possible that irrigation would be simulated on the last 

day of the simulated irrigation season or that irrigation 

would have been simulated on the next day, if the end of 

the season was not simulated, a difference of about 20 to 25 

percent (based on events occurring 4–5 times a year) would 

result. Because the purpose of the groundwater-flow model 

was to examine long-term signals, the natural component 

of recharge (Bauer and Vaccaro [1990] regression equation) 

was subtracted from total SOWAT-estimated recharge, and 

the irrigation component was smoothed and assumed to 

increase monotonically due to the increasing area of irrigation 

until 2007.
Numerical Simulation of Groundwater Flow 33

The mean annual natural groundwater recharge for the 

33,000 mi2 model domain for 1895–2007 was 3.2 in/yr 

(5.7 million acre-ft/yr) which was used as the estimate of 

recharge for the predevelopment steady-state model (fig. 14A). 

During model calibration, large interannual differences in 

natural recharge estimates caused model instabilities. This 

was attributed to the fact that the groundwater-flow model 

does not simulate a vadose zone, which would dampen the 

highly variable, precipitation-dependent recharge estimated 

using the Bauer and Vaccaro (1990) method. To add stability 

while retaining decadal trends in precipitation, annual recharge 

from precipitation was averaged with the preceding 2 years, 

resulting in a 3-year rolling average that was applied to the 

recharge estimates at each 1-km cell. Averaging with preceding 

years approximated vadose zone lag because high recharge 

years were damped and spread across several years. By 

2007, recharge associated with surface-water irrigation was 

estimated as 1.0 million acre-ft/yr and recharge associated 

with groundwater irrigation was estimated as 0.61 million 

acre-ft/yr. During model calibration, surface-water irrigation 

efficiency was increased from 50 to 75 percent, similar to 

assumed groundwater irrigation efficiency. This adjustment 

reduced recharge due to surface-water irrigation over time. 

Approximately 21 percent more recharge was entering the 

aquifer system during 2000–2007 (fig. 14B), than during the 

predevelopment period.

Streams and Surface-Water Features

The exchange of groundwater and surface water is an 

important hydrologic process in the groundwater-flow system 

and, to the extent possible for a regional model, the CPRAS 

model was constructed to simulate this process. Rivers and 

streams throughout the model domain were simulated with 

the MODFLOW River (RIV) and Drain (DRN) Packages 

(Harbaugh, 2005). Surface water was allowed to exchange with 

groundwater at all mapped surface-water features and where 

groundwater could seep out from exposed geologic units.

The RIV Package was used to simulate natural stream 

and river reaches with calculated mean annual flows from the 

National Hydrography Dataset (NHD; Simley and Carswell, 

2009) equal to or greater than 25 ft3/s throughout the model 

domain (fig. 10). The DRN Package was used to simulate 

all other reaches and areas where geologic units that are 

exposed at land surface can drain into nearby surface waters 

(fig. 10). To calibrate the model, river and drain fluxes were 

summed and compared to measured and estimated base flows. 

River parameters required by the RIV Package include river 

stage (STAGE), river bottom elevation (RBOT), and river 

conductance (computed as an adjustable multiplier times the 

parameter CONDFACT). Drain parameters required by the 

DRN package include the drain elevation (ELEVATION) 

and the drain conductance (an adjustable multiplier times the 

parameter CONDFACT). Initial values of riverbed conductance 

were based on stream length (determined using a geographic 

information system), and average depth and width for the 

river reaches were based on mean annual streamflow from the 

USGS National Hydrography Dataset (U.S. Geological Survey, 

2014) and regression equations determined by Magirl and 

Olsen (2009). The river and drain conductance parameters that 

control the effective rate of exchange of surface water with the 

aquifer system as a function of the head in the aquifer system, 

were assumed to be proportional to the area through which 

water can flow, and were adjusted during calibration using a 

multiplier. In model cells that contain multiple river and (or) 

drain reaches, those reaches were aggregated and formulated 

as a single river reach, drain reach, or both. The river and drain 

conductance parameters were assumed to be proportional to the 

cumulative area of all features within the cell that act as rivers 

or drains, respectively.

To estimate the cumulative area associated with all rivers 

and drains within each model cell, the NHD dataset was 

intersected with the model grid in plan view, and stream reach 

lengths and NHD properties were associated with the model 

grid. For each model row and column, if a reach exceeded the 

annual flow threshold of 25 ft3/s, then a RIV cell was simulated 

at the model layer where the lowest downstream reach elevation 

existed (plus a 0.1 ft upward offset to ensure no computational 

problems when rivers and drains have head controlling 

elevations too close to the NWT dry cell condition). The lowest 

elevation downstream reach was specified as the value for the 

RBOT, and STAGE was set to this elevation plus the associated 

depth of water estimated by Magirl and Olsen (2009) for the 

NHD. For each stream reach simulated using the RIV Package, 

area was estimated as the length of the wetted perimeter 

(computed from NHD estimates of top and bottom width 

and depth of water) times the length of the reach intersecting 

the model grid. All of these areas were summed to estimate 

CONDFACT for the RIV cell. If the maximum mean annual 

flow from the NHD at any MODFLOW row and column was 

less than 25 ft3/s, then the surface-water feature was simulated 

using the DRN package with CONDFACT computed in the 

same way and ELEVATION computed as the elevation of the 

lowest downstream reach plus the associated estimate of depth 

of water. For all model cells overlying these lowest cells, DRN 

cells also were simulated with ELEVATION set to the elevation 

of the top of each model cell. This allowed groundwater to 

leave the system by draining from the layers that are throughcut by stream incision. Setting the elevation at the cell top 

improved model stability by limiting dewatering and was 

consistent with the notion that each layer has an interflow zone. 

The area estimate for computation of CONDFACT assumed 

that the drains would be exposed in the canyon walls, and was 

computed for each reach as the distance from the top of the cell 

(plus 0.1-ft offset) minus the downstream reach elevation (or 

the cell bottom, whichever is greater) times the reach length 

times two (representing the two exposed faces in the canyon 

cut) times 10 percent (the estimated total thickness occupied by 

the interflow zone versus the flow interior). This formulation 

for drain area allowed that drainage was more efficient for 

stratigraphically high units intersected by more stream reaches, 

than with low units intersected less frequently.
34 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

tac13-0856_fig 14a

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

EXPLANATION Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

EXPLANATION

Active model domain

Predevelopment recharge, 

in inches per year

0

0.1 – 1.0

1.1 – 2.0

2.1 – 3.0

3.1 – 4.0

4.1 – 5.0

5.1 – 10.0

10.1 – 25.0

25.1 – 50.0

50.1 – 61.1

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

A

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 14. Mean annual groundwater recharge for (A) predevelopment and (B) current conditions (2000–2007), Columbia Plateau 

Regional Aquifer System, Idaho, Oregon, and Washington.
Numerical Simulation of Groundwater Flow 35

tac13-0856_fig 14b

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

EXPLANATION Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

EXPLANATION

Active model domain

Current condition (2000–2007) recharge,

in inches per year

0.1

0.2 to 1.0

1.1 to 2.0

2.1 to 3.0

3.1 to 4.0

4.1 to 5.0

5.1 to 10.0

10.1 to 20.0

20.1 to 40.0

40.1 to 46.5

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

B

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 14.—Continued
36 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

The final locations requiring drains were at geologic 

outcrops controlled by faulting or strong folds combined with 

erosion and where simulation of horizontal flow barriers along 

stream valleys might prevent drainage from each side of the 

valley into the streams. All of these features were represented 

by the HFB Package (see section, “Geologic Structures”), 

so during model construction, both sides of each HFB were 

checked to determine if a stratigraphic layer was exposed in 

outcrop on one side of the HFB. If this condition occurred, a 

DRN cell was simulated with ELEVATION set to the elevation 

of the top of the model cell and CONDFACT equal to the 

thickness of the model cell times the effective length of the 

outcrop times the fraction of the thickness contributing to the 

flow (for erosional cuts, this was assumed to be 10 percent). 

The product of the effective length times the thickness fraction 

was selected by iteration with an early version of the model, 

and a value of 5,000 ft provided robust and stable solutions 

that prevented anomalous hydraulic heads near fault scarps. 

Uncertainty was not evaluated, but the probability that 

drainage occurs from the broken uplifted side of a fault is 

consistent with the conceptual model of flow through CRBG 

units.

This approach added numerical stability to the model 

in areas far from regions of greatest interest and accounted 

for the generally gaining stream reaches. The drain hydraulic 

conductance is a function of the surrounding hydrogeologic 

material and the drain geometry. Information necessary 

to calculate a drain conductance, such as the distribution 

and hydraulic conductivity of material near the drain, were 

not available. Commonly, drain conductance is a lumped 

parameter that is adjusted during calibration to match 

measured flows. The CPRAS model included 43,311 drain 

cells assigned from model layers 1–64 and 1,727 river cells 

assigned to model layers 2–60.

Model Calibration and Sensitivity

Model calibration is the process in which model 

parameters, model structure, and boundary conditions are 

adjusted to obtain a reasonable fit between simulated heads 

and fluxes and measured data. The final CPRAS model used 

the principle of parsimony, where the simplest model that 

provides a good fit between observed data and simulated 

equivalents is retained. If the model fit is determined to 

be “good enough” by the modeler, taking into account the 

purpose of the model, limitations imposed by the modeling 

assumptions, and quality and distribution of the data 

supporting estimation of model parameters then the model 

is considered to be adequately calibrated. Initially, a simple 

distribution of parameters was used, followed by addition of a 

large number of parameters that allowed testing of the role of 

geologic structure on groundwater flow, followed by a process 

of simplifying the model by grouping parameters. Highly 

parameterized models require a different approach to model 

calibration and uncertainty, such as regularized inversion 

(Doherty and Hunt, 2010).

The model was calibrated using the iterative parameter 

estimation (PEST) software package (Doherty, 2010). 

PEST uses a nonlinear least-squares regression to find the 

set of parameter values that minimizes a weighted sum-ofsquared-errors objective function. Throughout the calibration 

process, no adjustments were made that conflicted with the 

general understanding of the aquifer system and previously 

documented information.

Although the results of the predevelopment steadystate and long-term transient models are presented here in 

logical order, in reality, the models were calibrated iteratively. 

Model run times associated with the transient (1900–2007) 

simulation period and the complexities added by the MNW2 

Package, increased recharge from irrigation, and addition of 

storage properties, made automated calibration of the transient 

model impractical. Parameter-estimated values and hydraulic 

heads from the predevelopment model were used as starting 

parameter values and initial conditions in the transient model 

runs. After the transient model was evaluated, changes were 

made to the steady-state model (for example, locations of 

RIV, DRN, and HFB cells) and the automated calibration of 

the steady-state model was repeated. In this way, the utility 

of automated parameter estimation and the computational 

efficiency of trial-and-error calibration were used to create 

well-constrained models.

Observations Used in Model Calibration 

Water-Level Altitudes, Water-Level Altitude Changes, and 

Associated Errors

The winter median groundwater levels described by 

Burns and others (2012b) were used for model calibration. 

For purposes of this study, the winter median level is defined 

as the median value of all water-level measurements collected 

during the first 3 months of each calendar year. These 

measurements were selected because they are not strongly 

affected by seasonal pumping and provide a record of longerscale hydraulic response to system stresses.

Because most of the water-level measurements were 

taken in or near commingling wells (Burns and others, 2012b), 

the standard approach of defining error-based observation 

weights (Hill, 1998) was not appropriate because the larger 

part of the error in much of the model domain is due to bias 

in the measurement resulting from commingling. Burns and 

others (2012b) identified measurement noise (attributed to 

commingling and geologic complexity) of about ±200 ft over 

much of the Palouse Slope. The large range in water levels 

from relatively similar locations and stratigraphic positions 

are shown in figure 15 (well groups are explained in section, 

“Transient Model Fit and Model Error—Comparison of 

Simulated and Measured Hydraulic Heads”). Discontinuous 
Numerical Simulation of Groundwater Flow 37

interflow zones, differing times of well construction, and 

thus commingling, and different screen intervals in the 

completed wells further complicated information that could 

be obtained from the water level. Most correctly, for the 

CPRAS each water-level measurement should be considered 

a complex flow-weighted average of hydraulic heads in all 

aquifers intersected by the well or neighboring wells. This 

complex interpretation of each hydraulic-head measurement 

explains why water-level measurements in nearby, similarly 

constructed wells can be hundreds of feet different (fig. 15; 

Burns and others, 2012b). Given this complexity in 

interpreting individual hydrographs, it was assumed that 

each observation had a similar amount of error, and the error 

increased over time as commingling became more prevalent.

For calibration of the steady-state model, well screen 

information was not available, so the earliest water-level 

measurement at each well was used with the value of the 

observation being assigned to the model layer with the top 

closest to the well bottom altitude. This assumes that most 

aquifers are associated with interflow zones, and that the 

driller stopped at an altitude where a productive aquifer was 

detected. In areas of significant upward hydraulic gradient, 

the simulated head likely would be higher than the measured 

commingled head, and conversely in areas of downward 

gradient. The total number of head observations used for the 

steady-state and transient model calibrations were 10,525 and 

46,460, respectively.

Examination of well hydrographs (Burns and others, 

2012b) indicated that the first widespread water-level changes 

was in response to large-scale surface-water irrigation 

projects about 1950. Next, about 1970, many areas across 

the CPRAS started to exhibit strong declines associated with 

rapid development of the groundwater resource. Because the 

earliest observation at a well was possibly during the period 

of transient hydraulic response, the uncertainty associated 

with using late-time observations to calibrate the steady-state 

groundwater model was represented by selecting lumped 

weights based on response periods identified in the data. 

Observation weights reduce the influence of observations that 

are less accurate and increase the influence of observations 

that are more accurate (Hill, 1998). A weight of 1.0 was 

assigned for pre-1950 measurements, 0.01 for pre-1970 

measurements, and 0.0001 for post-1970 measurements. This 

weighting strategy ensured that the earliest-time data in each 

area of the model domain would have the most influence on 

parameter estimates.

Winter median water-level elevation, in feet (NAVD88)

Slope eastern well groups

Stratigraphic position relative to top of the Grande Ronde basalt, in feet

Well bottom in

the Grande Ronde Basalt 

Well bottom above

top of Grande

Estimated top of

Ronde Basalt

Grande Ronde Basalt

700

800

900

1,000

1,100

1,200

1,300

1,400

1,500

1,600

1,700

1,800

1,900

-2,200 -2,000 -1,800 -1,600 -1,400 -1,200 -1,000 -800 -600 -400 -200 0 200 400 600

EXPLANATION

Well

Group 9

Group 10

Group 8

Group 11 - Deep

Group 11 - Shallow

Group 12

Figure 15. Distance of the well bottom from the estimated top of Grande Ronde Basalt 

compared with water levels measured in wells for the Palouse Slope eastern well groups, 

Washington. From Burns and others (2012b, fig. 20E).
38 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Following steady-state model calibration, a transient 

model was developed for the period 1900–2007. An initial 

1900–1919 transient stress period was simulated using pre1920 commingling wells with pumping rates increased from 

zero starting in 1920. Predevelopment steady-state model 

heads were used as initial conditions for the transient model. 

Using the transient 1900–1919 stress period provided a better 

match of simulated to measured water levels than using a 

steady-state 1920s condition model. Use of a steady-state 

stress period for 1920 resulted in significant head declines due 

to pumping in structurally complex basins, especially those 

with strong upward hydraulic gradients such as the Yakima 

and Walla Walla River Basins. Because pre-1920 pumps 

had poor lift characteristics (Cline and Collins, 1992), the 

historical use of artesian aquifers in these basins is consistent 

with the conceptual model and the drilling records.

After reasonably good steady-state model fit was 

achieved using lumped weights, preliminary transient 

simulations were completed, and computed post-2000 

drawdowns were used to adjust steady-state model calibration 

weights. If large transient simulated drawdowns were 

detected (usually in the most heavily anthropogenically 

influenced areas), then the observation target was determined 

to be less useful for steady-state calibration; therefore, 

the steady-state observation weight was reduced. An 

exponential decay function was used to decrease observation 

weight such that for every 50 ft of computed drawdown 

(magnitude), observation weight was reduced by an order 

of magnitude (fig. 16). The iterative process for the steadystate and transient model calibration improved confidence 

that higher weight measurements likely were closer to 

predevelopment conditions.

Several limitations of the steady-state model calibration 

are implied: 

• Few, if any of the water-level measurements were taken 

during true predevelopment conditions, but they might 

still represent predevelopment, or near predevelopment 

conditions.

• Assignment of the measured hydraulic head to the 

deepest model layer is not strictly correct for a 

commingled measurement.

• Predevelopment vertical gradients are poorly reflected 

in the measured data because most observations 

represent some degree of commingling, commonly 

with an apparently random distribution of high and 

low measured values in many areas (Burns and others, 

2012b). Because few if any data constrain high vertical 

gradients, the calibrated predevelopment model will 

have the tendency to overestimate Kv, matching 

measured heads on average in many regions.

Subject to these limitations, the steady-state model 

calibration is best at providing good general estimates of 

model parameters for use with the transient model, and 

adjustments in model parameters that improve transient model 

fit should take precedence over steady-state model-calibration 

parameter estimates. The transient model incorporated the 

addition of stresses that provided information about aquifer 

properties that could not be explored in a non-stressed 

system and used more recent, and presumably more accurate, 

records of stresses (pumping and recharge) and responses 

(water levels).

Use of the water-level measurements to calibrate the 

transient model is subject to the same limitations regarding 

measuring water levels in commingling wells, but the 

limitations associated with the fact that few observations 

truly represent predevelopment conditions is removed. To 

address the commingling problem, Burns and others (2012b) 

showed that locally, hydraulic gradients were frequently much 

more variable vertically than horizontally. This indicates that 

simulated heads spanning many layers of a row and column 

should produce a range of heads that, when averaged by layer, 

can reasonably produce the range of heads observed. Because 

this method of using the observation data to evaluate model 

performance is complex, because transient model run times 

were long, and because further tightening of Kv of CRBG 

units resulted in model instability, manual calibration and 

automated parameter estimation, was accomplished for the 

transient model. However, the spatial patterns of differences 

between simulated and measured water levels for selected 

periods were examined (under the assumption that measured 

head represented head at the altitude of the well bottom). 

Ideally, multi-layer observations would have been used to 

constrain the transient model calibration. However, attempting 

to simulate open screen intervals in the appropriate aquifer 

units would have introduced an unknown and unquantified 

level of uncertainty.

Streamflow Observations and Associated Errors

The emphasis of the CPRAS model is the groundwater 

system, and streamflow observations are important to augment 

the water-level information and further constrain the parameter 

estimation about groundwater-surface-water interaction. 

Correctly simulating stream base flow (the groundwater 

discharge component of streamflow) is an important way of 

ensuring the correct amount of water is moving through the 

system. The modeling effort did not attempt to simulate the 

heavily regulated and managed surface-water features of the 

study area, such as reservoir releases or diversions and returns. 

Estimated or measured base flows were compared to simulated 

groundwater discharge to selected stream reaches.
Numerical Simulation of Groundwater Flow 39

tac13-0856_fig16

Steady-state hydraulic-head observation weight

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1890 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000 2010

Year

Figure 16. Hydraulic-head observation weights used in the steady-state model, Columbia Plateau Regional Aquifer System, 

Idaho, Oregon, and Washington.
40 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Two methods were used to estimate base flows at 

50 locations (fig. 17; table 3). The most reliable estimates 

of base flow were computed from an automated hydrograph 

separation technique (HYSEP; Sloto and Crouse, 1996) used 

to evaluate the groundwater contribution to total streamflow 

(base flow) at active and inactive USGS streamflow-gaging 

stations in Washington State (applied by Sinclair and 

Pitz, 1999). The estimates by Sinclair and Pitz (1999) are 

basin-scale averages for the period of record; therefore, 

the estimates were considered appropriate observations for 

model calibration.

For streamflow-gaging stations not included in Sinclair 

and Pitz (1999), base flow was estimated as 62 percent of 

annual median streamflow. This approach was based on 

the findings by Sinclair and Pitz (1999) that groundwater 

discharge in eastern Washington averaged about 62 percent of 

total annual median streamflow. All estimates by this method 

were checked to ensure the estimate was greater than the 

minimum monthly streamflow.

Base-flow estimates by the two methods (62 percent of 

annual median streamflow and methods of Sinclair and Pitz 

[1999]) provided a general idea of groundwater contributions 

to streamflow, but were not appropriate to use as a specific 

model observation because of the variability of recharge and 

effects of human activities on streamflow. Therefore, an upper 

and lower bound around the base-flow estimate were used to 

calibrate the model. For estimates based on the hydrograph 

method of Sinclair and Pitz (1999), the bounds were an 

order of magnitude greater than and less than the base-flow 

estimate. For estimates based on 62 percent of annual median 

streamflow, the upper bound was the median annual flow 

and the lower bound was the minimum monthly flow. If the 

simulated streamflow gain was within the upper and lower 

bounds, the observation did not contribute to the objective 

function and no penalty was assessed during automated 

calibration of the steady-state model.

Model Evaluation

Parameterization and Regularization

To test the effect of geologic variability on groundwater 

flow, the model domain was divided into zones of similar 

geologic character. The division of units created the 

possibility of 24 plan-view model parameter zones (fig. 18) 

and 8 vertical zones (5 geologic units plus 3 confining units), 

but not all of the vertical zones existed in each of the 24 

plan-view zones. Within each parameter zone, Kh and Kv, 

and river and drain conductance multipliers were selected as 

independent parameters.

The model then was parameterized to test the importance 

of geologic features that had been simulated in the CP-RASA 

model (Hansen and others, 1994). The following parameter 

zones were created for testing:

• Plan-view zones: Burns and others (2012b) observed 

that in areas where data were too sparse to resolve 

individual horizontal flow barriers, the net effect of 

increasing geologic structure was a corresponding 

decrease in Kh. Zones were created by dividing the 

model domain into model cell zones with similar 

observable patterns of structure (fig. 18). To prevent 

the possibility of an ill-posed parameter estimation 

problem, Tikhonov Regularization (Doherty, 2010) 

was used to relate zones of similar structural patterns, 

indicating that properties in these zones are likely to 

have similar parameter values.

• Vertical (geologic unit) zones: Within each plan-view 

zone, it was postulated that each of the geologic units 

of Burns and others (2011) could have markedly 

different hydraulic characteristics because of changes 

in deposition. Tikhonov Regularization was used to 

inform the parameter estimation process that, unless 

observations support otherwise, each of the three basalt 

units should have similar hydraulic properties and 

each of the sedimentary units (Overburden and two 

interbeds) should have similar hydraulic properties. 

This parameterization allowed testing of the postulated 

reduction in hydraulic conductivity with increasing 

depth and testing of the role of depositional style on 

hydraulic conductivity.

• Zones with larger vertical conductance (basalt 

units only): Basalt flow interiors were postulated 

to be confining units even though jointing patterns 

frequently crosscut most of the dense flow interior. 

Although tight sealing of joints is expected where 

overlying rock compresses flow interiors, it was 

postulated that in the absence of large amounts of 

overlying rock, joints could provide enhanced vertical 

conductivity, enhancing areal recharge and discharge in 

river valleys. This mechanism could explain the pattern 

of enhanced vertical permeability estimated near major 

rivers in the work of Davies-Smith and others (1988) 

and Hansen and others (1994). To test this mechanism, 

the uppermost two basalt model layers (about 200 ft in 

thickness) were allowed to have a different Kv than the 

underlying units. Tikhonov Regularization was used to 

inform the parameter estimation process that unless the 

observation data indicated otherwise, the Kv of each 

cell should be similar to the deeper parts of the same 

geologic unit.

In addition to the spatial parameter zones, 25 discrete 

HFB categories were parameterized, generally based on 

mapped geologic structure. The adjustable parameter, 

barrier conductance, was assumed to be constant along 

the fault and with depth. Tikhonov Regularization was 

used to inform the parameter estimation process that fault 

conductance in a region should have similar values unless data 

support differences.
Numerical Simulation of Groundwater Flow 41

tac13-0856_fig 17

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

Paradise

Creek

9

7 8

6

5

4

1 3 2

46

44

41

34

30

28 27

24

22

19

12

10

33

32

36

35

39

42

40

38

43

31

26

29

18

20

37

25

23

48

21

47

49

45

15

50

11

16

13

17

14

Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

EXPLANATION

Active model domain

Base-flow observation site and No. 2

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Snake River

COLUMBIA RIVER

Umatilla River

Walla Walla River

Crab Creek

Study

area

WASHINGTON

OREGON

IDAHO

Figure 17. Locations of base-flow observation sites used to calibrate the steady-state and transient models of the Columbia 

Plateau Regional Aquifer System, Idaho, Oregon, and Washington.
42 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington Table 3. Mean annual base-flow estimates, and upper and lower bounds used to calibrate the steady-state and transient models of the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington. [Site locations are shown in figure 17. Station name: OR, Oregon; WA, Washington; ID, Idaho; NR, near; R, river, W, west; abv, above; BNDY, boundary; CR or C, creek; SO, south, FK, fork. ft3/s, cubic foot per second] Base-flow observation site Station name Base-flow estimate (ft3/s) Lower bound (ft3/s) Upper bound (ft3/s) Weight Source

1 WHITE RIVER BELOW TYGH VALLEY, OR 234.70 118.9 378.6 10 62 percent of annual median streamflow2 RHEA CREEK NEAR HEPPNER, OR 16.09 4.0 26.0 10 62 percent of annual median streamflow3 ROCK CREEK AB WHYTE PARK NR CONDON, OR 42.28 2.0 68.2 10 62 percent of annual median streamflow4 WILLOW CREEK ABV WILLOW CR LAKE, NR HEPPNER, OR 11.66 1.2 18.1 10 62 percent of annual median streamflow5 BUTTER CREEK NEAR PINE CITY, OR 22.89 5.2 36.9 10 62 percent of annual median streamflow6 MCKAY CREEK NEAR PILOT ROCK, OR 71.21 1.4 114.9 10 62 percent of annual median streamflow7 UMATILLA R AT W RESERVATION BNDY NR PENDLETON, OR 322.51 41.5 520.2 10 62 percent of annual median streamflow8 MOONSHINE CREEK NEAR MISSION, OR 1.77 0.1 2.9 10 62 percent of annual median streamflow9 UMATILLA RIVER AT PENDLETON, OR 299.91 45.6 483.7 10 62 percent of annual median streamflow

10 UMATILLA RIVER ABOVE MEACHAM CREEK, NR GIBBON, OR 138.40 44.1 223.2 10 62 percent of annual median streamflow

11 ROCK CREEK NEAR ROOSEVELT, WA 17.00 1.7 170.0 1 Sinclair and Pitz, 1999

12 SOUTH FORK WALLA WALLA RIVER NEAR MILTON, OR 104.38 102.2 168.4 10 62 percent of annual median streamflow13 SPRING CREEK NEAR BLOCKHOUSE, WA 14.00 1.4 140.0 1 Sinclair and Pitz, 199914 LITTLE KLICKITAT R NR WAHKIACUS, WA 103.00 10.3 1,030.0 1 Sinclair and Pitz, 199915 ALDER CR AT ALDERDALE, WA 3.00 0.3 30.0 1 Sinclair and Pitz, 199916 LITTLE KLICKITAT R NR GOLDENDALE, WA 35.00 3.5 350.0 1 Sinclair and Pitz, 199917 MILL CREEK NEAR BLOCKHOUSE, WA 14.00 1.4 140.0 1 Sinclair and Pitz, 199918 MILL CREEK NEAR WALLA WALLA, WA 69.00 6.9 690.0 1 Sinclair and Pitz, 199919 WALLA WALLA RIVER NEAR TOUCHET, WA 346.03 18.4 558.1 10 62 percent of annual median streamflow20 BLUE CREEK NEAR WALLA WALLA, WA 10.00 1.0 100.0 1 Sinclair and Pitz, 199921 TOUCHET R NR TOUCHET, WA 148.00 14.8 1,480.0 1 62 percent of annual median streamflow22 MILL CREEK NEAR WALLA WALLA, WA 47.99 1.1 77.4 10 62 percent of annual median streamflow23 DRY CREEK NEAR WALLA WALLA, WA 13.00 1.3 130.0 1 Sinclair and Pitz, 199924 YAKIMA RIVER AT KIONA, WA 1,977.48 1,489.4 3,189.5 10 62 percent of annual median streamflow25 TOUCHET RIVER AT BOLLES, WA 145.00 14.5 1450.0 1 Sinclair and Pitz, 199926 EAST FK TOUCHET R NR DAYTON, WA 95.00 9.5 950.0 1 Sinclair and Pitz, 199927 ASOTIN CREEK BELOW CONFLUENCE NEAR ASOTIN, WA 30.20 23.9 48.7 10 62 percent of annual median streamflow28 GRANGER DRAIN AT GRANGER, WA 21.61 19.5 34.9 10 62 percent of annual median streamflow29 TUCANNON RIVER NEAR STARBUCK, WA 127.00 12.7 1,270.0 1 Sinclair and Pitz, 199930 AHTANUM CREEK AT UNION GAP, WA 42.41 15.4 68.4 10 62 percent of annual median streamflow31 MEADOW CREEK NR CENTRAL FERRY, WA 1.90 0.2 19.0 1 Sinclair and Pitz, 199932 SO FK PALOUSE R ABV PARADISE C NR PULLMAN, WA 7.20 0.7 72.0 1 Sinclair and Pitz, 199933 PARADISE CR NR PULLMAN, WA 5.10 0.5 51.0 1 Sinclair and Pitz, 199934 PARADISE CR AT UNIVERSITY OF IDAHO AT MOSCOW, ID 4.33 0.8 7.0 10 62 percent of annual median streamflow35 SOUTH FORK PALOUSE RIVER AT PULLMAN, WA 18.13 3.8 29.2 10 62 percent of annual median streamflow
Numerical Simulation of Groundwater Flow 43

Table 3. Mean annual base-flow estimates, and upper and lower bounds used to calibrate the steady-state and transient models of the Columbia Plateau Regional Aquifer 

System, Idaho, Oregon, and Washington.—Continued

[Site locations are shown in figure 17. Station name: OR, Oregon; WA, Washington; ID, Idaho; NR, near; R, river, W, west; abv, above; BNDY, boundary; CR or C, creek; SO, south, FK, fork. ft3/s, cubic foot 

per second]

Base-flow 

observation 

site

Station name

Base-flow

estimate

(ft3/s)

Lower

bound

(ft3/s)

Upper

bound

(ft3/s)

Weight Source

36 MISSOURI FLAT CREEK AT PULLMAN, WA 3.80 0.4 38.0 1 Sinclair and Pitz, 1999

37 PALOUSE RIVER AT HOOPER, WA 333.53 37.9 538.0 10 62 percent of annual median streamflow

38 UNION FLAT CREEK NEAR COLFAX, WA 18.00 1.8 180.0 1 Sinclair and Pitz, 1999

39 FOURMILE CR AT SHAWNEE, WA 4.80 0.5 48.0 1 Sinclair and Pitz, 1999

40 PALOUSE RIVER BELOW SOUTH FORK AT COLFAX, WA 211.63 21.9 341.3 10 62 percent of annual median streamflow

41 PALOUSE RIVER NR POTLATCH, ID 152.88 11.2 246.6 10 62 percent of annual median streamflow

42 PALOUSE RIVER NEAR COLFAX, WA 175.00 17.5 1,750.0 1 Sinclair and Pitz, 1999

43 PINE CREEK AT PINE CITY, WA 20.00 2.0 200.0 1 Sinclair and Pitz, 1999

44 CRAB CREEK AT ROCKY FORD ROAD NEAR RITZVILLE, WA 22.09 11.9 35.6 10 62 percent of annual median streamflow

45 ROCKY FORD CREEK NEAR EPHRATA, WA 69.00 6.9 690.0 1 Sinclair and Pitz, 1999

46 CRAB CREEK AT IRBY, WA 20.82 4.9 33.6 10 62 percent of annual median streamflow

47 CRAB CREEK AT IRBY, WA 29.00 2.9 290.0 1 Sinclair and Pitz, 1999

48 COAL CREEK AT MOHLER, WA 1.20 0.1 12.0 1 Sinclair and Pitz, 1999

49 PARK CREEK BLW PARK LAKE NR COULEE CITY, WA 9.50 1.0 95.0 1 Sinclair and Pitz, 1999

50 DOUGLAS CREEK NEAR ALSTOWN, WA 0.92 0.1 9.2 1 Sinclair and Pitz, 1999
44 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

tac13-0856_fig 18

100

210

180

410

160

130

150

110

140

120

220

200

310

300

400

420

430

170

320

470

460

440

450 440

Pasco

Yakima

Moscow Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

Active model domain

Model parameter zone

Fault

Fold

110

EXPLANATION Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

EXPLANATION

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 18. Initial model parameter zones and geologic structure for the Columbia Plateau Regional Aquifer System, Idaho, 

Oregon, and Washington.
Numerical Simulation of Groundwater Flow 45

Optimal Parameter Estimates

Parameter estimation was iterative with 684 parameters 

being the most parameters estimated at one time for the 

predevelopment steady-state model. Results of these 

highly-parameterized PEST runs were used to evaluate 

the mechanisms identified in section, “Parameterization 

and Regularization.”

The 24 plan-view parameter zones (fig. 18) were 

incrementally combined manually into larger zones with 

constant horizontal conductivity until 4 distinct Overburden 

zones (fig. 19A) and 3 CRBG zones (fig. 19B) remained. Zones 

were combined when calibrated Kh was sufficiently close to 

the same value such that use of a single value caused only a 

minor increase in the weighted least-squares objective function. 

Final calibrated parameter values are presented in table 4.

Basalt Kh was highest near the center of the CPRAS 

where flows are the thickest and are less folded and faulted. At 

the margins of the lava flows, Kh of the relatively undeformed 

basalt units decreased slightly. As lavas continue to thin and 

are increasingly cross-cut by geologic structure on the flanks of 

the Cascades Range and the Blue Mountains, effective Kh of 

the basalts decreased greatly. From a regional perspective, the 

Overburden unit could be simulated even more simply than the 

basalt units. All Overburden in the Yakima River Basin could 

have a single homogeneous value (Kh = 99.81 ft/d) and all 

other Overburden values could range from 0.35 ft/d (zone 1) to 

100 ft/d (zone 3). The Mabton and Vantage Interbeds proved 

to be negligibly thin so they were lumped with the overlying 

basalt unit in each case. The distribution of Kv was even 

simpler, with a uniform value of 5.24 × 10-3 ft/d used for the 

upper 200 ft of the CRBG units and 1.00 × 10-4 ft/d for the 

rest of the CRBG (table 4). Overburden Kv was tied to the 

Kh at a fixed ratio of 100 times less than Kh. This relatively 

simple model of the Overburden unit conductivity was used to 

simulate the net effect of the relatively sparse overburden on 

the CRBG aquifers, and little effort was placed on refining and 

improving model fit through modifications of the overburden 

aquifer parameters.

Parameter estimation was performed for storage 

properties, but the model proved to be relatively insensitive 

to changes in those parameters. Trial and error manual 

calibration and values from previous studies were evaluated to 

provide initial parameter values for the transient calibration. 

Throughout the calibration process, the storage parameters 

occasionally were reassessed and changed to improve model 

results. Final values for specific yield for the Overburden and 

CRBG units were 1.0 × 10-1 and 2.5 × 10-2, respectively. Final 

values for specific storage for the Overburden and CRBG units 

were 2.50 × 10-5 ft-1 and 2.50 × 10-6 ft-1, respectively (table 4).

When all 684 parameters were estimated independently, 

there was no clear pattern of Kh with depth between CRBG 

units. Even though prior information was used to inform the 

automated parameter estimation process that CRBG units 

have similar hydraulic properties to other nearby CRBG units, 

some plan-view zones showed decreasing Kh with depth, 

whereas others displayed the opposite behavior. This was 

attributed to tradeoffs between river and drain parameters 

and Kh resulting in small improvements to model fit at the 

expense of highly irregular estimates of river and drain 

conductance. Based on preliminary calibration, river and drain 

conductance multipliers were combined into three major zones 

for much of the calibration: all Overburden, Yakima River 

Basin CRBG units, and all other CRBG units. Generally, the 

resulting multipliers were of similar magnitude, especially 

when comparing similar rock types, indicating that scaling 

the conductance to stream geometry factors worked well. 

Combining river and drain parameters into these three groups 

allowed the CRBG units to be combined into four groups with 

similar Kh, after which drain multipliers were allowed to vary 

by zone to fine tune the model fit. Lastly, the Overburden river 

conductance was adjusted to fine-tune river exchange in the 

deep sedimentary valleys.

When evaluating the initial estimates of Kv for the 

parameter zones “zones with larger vertical conductance 

(basalt units only)” with 684 independent parameters, 

automated parameter estimation frequently tended to make 

these units less permeable (lower Kv) than corresponding 

underlying layers, which contradicted the expected behavior. 

Following the grouping of similar zones, however, the 

calibrated Kv value for the upper 200 ft of basalt was 

approximately 1.5 orders of magnitude greater than the 

underlying basalts. Other than the upper 200 ft and unlike the 

previous CP-RASA model (Hansen and others, 1994), the 

final model parameter distribution had no vertical trend in 

hydraulic conductivity in the underlying CRBG units. This 

lack of vertical trend is likely because the current model was 

connected to surface-water features only where thin model 

layers intersected the land surface. The CP-RASA model used 

five thick layers, requiring that Kh partially account for the 

lack of connection between deep parts of CPRAS and surfacewater features. In the new model, the Kv controlled this lack 

of connectivity.

During parameter estimation, the river and drain 

conductances were varied as a function of geology, similar to 

Burns and others (2012a). For each model cell, the length of 

a stream depends on the path across the cell, and conductance 

is linearly dependent on the path length. The adjustable 

multipliers to determine final river and drain conductances are 

presented in table 4.
46 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

tac13-0856_fig 19A

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

EXPLANATION Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

Model parameter zone

1

2

3

4

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

A

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 19. Simplified final model parameter zones for the (A) Overburden unit and (B) Columbia River Basalt Group, Columbia 

Plateau Regional Aquifer System, Idaho, Oregon, and Washington.
Numerical Simulation of Groundwater Flow 47

tac13-0856_fig 19b

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

EXPLANATION Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

Model parameter zone

1

2

3

Horizontal-flow barrier

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

B

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 19.—Continued
48 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Table 4. Final calibrated model parameters for the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington.

[Unit: OB, Overburden unit; CRBG, Columbia River Basalt Group. Complex parameter zones are shown in figure 18. Simplified parameter zones are shown 

in figure 19A (OB) and figure 19B (CRBG). HFB groups are shown in figure 9. Abbreviations: HFB, horizontal-flow barrier; ft, foot; ft-1, per foot; ft/d, foot 

per day; ft2/d, foot squared per day; (ft/d)/ft, foot per day per foot]

Unit

Horizontal hydraulic conductivity (ft/d)

Complex parameter zone

Simplified 

parameter 

zone

Minimum 

allowable 

value

Parameter

value

Maximum 

allowable 

value

OB 1.00E-04 3.05 1.00E+03 100, 110, 120, 130, 140, 150, 170, 180 OB 1

OB 1.00E-04 30.15 1.00E+03 160 OB 2

OB 1.00E-04 200 1.00E+02 200, 210, 220, 300, 310, 320 OB 3

OB 1.00E-04 150.81 1.00E+03 400, 410, 420, 430, 440, 450, 460, 470 OB 4

CRBG 1.00E-04 1.07 1.00E+06 110, 140, 150, 210, 220, 310, 410, 420, 430, 440, 470 CRBG 1

CRBG 1.00E-04 8.16 1.00E+06 120, 130, 170, 180, 200, 300, 320, 400, 450, 460 CRBG 2

CRBG 1.00E-04 27.90 1.00E+06 100, 160 CRBG 3

Unit

Vertical hydraulic conductivity (ft/d)

Complex parameter zone

Simplified 

parameter 

zone

Minimum 

allowable 

value

Parameter

value

Maximum 

allowable 

value

OB 1.00E-07 3.05E-02 1.00E+02 100, 110, 120, 130, 140, 150, 170, 180 OB 1

OB 1.00E-07 3.01E-01 1.00E+02 160 OB 2

OB 1.00E-07 2.00 1.00E+02 200, 210, 220, 300, 310, 320 OB 3

OB 1.00E-07 1.51 1.00E+02 400, 410, 420, 430, 440, 450, 460, 470 OB 4

CRBG–above 200 ft 5.00E-05 5.24E-03 1.00E+02 All zones All

CRBG–below 200 ft 1.00E-10 1.00E-04 1.00E-03 All zones All

Unit

Specific yield 

(dimensionless)

Parameter

value

OB 0.10

CRBG 0.025

Unit

Specific storage 

(ft-1)

Parameter

value

OB 2.50E-05

CRBG 2.50E-06

Unit

Drain conductance (ft2/d)

Minimum

value

Mean

value

Maximum 

value

OB 1.00 84.20 125.00

CRBG 0.20 92.50 200.00

Unit

River conductivity (ft/d)

Minimum

value

Maximum 

value

OB 0.01 0.57

CRBG 0.01 0.01
Numerical Simulation of Groundwater Flow 49

Table 4. Final calibrated model parameters for the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington.—Cont.

[Unit: OB, Overburden unit; CRBG, Columbia River Basalt Group. Complex parameter zones are shown in figure 18. Simplified parameter zones are shown 

in figure 19A (OB) and figure 19B (CRBG). HFB groups are shown in figure 9. Abbreviations: HFB, horizontal-flow barrier; ft-1, per foot; ft/d, foot per day; 

ft2/d, foot squared per day; (ft/d)/ft, foot per day per foot]

Unit

HFB hydraulic characteristic (1/d)

HFB

group

Minimum 

allowable 

value

Parameter

value

Maximum 

allowable 

value

CRBG 1.00E-16 1.00E-15 1.00E+03 1

CRBG 1.00E-16 2.11E-09 1.00E+03 10

CRBG 1.00E-16 1.20E-14 1.00E+03 21

CRBG 1.00E-16 8.67E-08 1.00E+03 22

CRBG 1.00E-16 1.16E-05 1.00E+03 23

CRBG 1.00E-16 4.53E-06 1.00E+03 32

CRBG 1.00E-16 1.94E-05 1.00E+03 33

CRBG 1.00E-16 1.20E-10 1.00E+03 34

CRBG 1.00E-16 2.30E-11 1.00E+03 35

CRBG 1.00E-16 5.96E-11 1.00E+03 41

CRBG 1.00E-16 1.93E-14 1.00E+03 51

CRBG 1.00E-16 3.88E-05 1.00E+03 52

CRBG 1.00E-16 4.26E-15 1.00E+03 53

CRBG 1.00E-16 8.82E-11 1.00E+03 54

CRBG 1.00E-16 6.53E-14 1.00E+03 61

CRBG 1.00E-16 4.94E-11 1.00E+03 62

CRBG 1.00E-16 6.12E-11 1.00E+03 63

CRBG 1.00E-16 1.17E-11 1.00E+03 64

CRBG 1.00E-16 5.19E-15 1.00E+03 65

CRBG 1.00E-16 1.44E-13 1.00E+03 66

CRBG 1.00E-16 1.15E-13 1.00E+03 67

CRBG 1.00E-16 7.31E-11 1.00E+03 68

CRBG 1.00E-16 1.05E-08 1.00E+03 69

CRBG 1.00E-16 8.94E-13 1.00E+03 70

CRBG 1.00E-16 2.80E-12 1.00E+03 71
50 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Examination of groundwater-level residuals revealed a 

few water levels measured in river valleys that would have 

been more closely matched by simulated water levels if Kv 

was significantly higher. However, upon closer examination 

of the measured hydraulic head in the valleys, the head was 

determined to be typically at the same approximate altitude 

as the trace of the axis of an anticline that crosses the river 

at a downgradient location. Because the geologic model of 

Burns and others (2011) did not capture all mapped folds (a 

result of low data density in some areas), the geologic model 

(and therefore the layers of the flow model) was frequently 

smoother than the mapped geology, indicating that deeper 

geologic strata may be exposed along the river than was being 

simulated in the groundwater-flow model. This indicated that 

the flow model included layers (and associated Kv) that did 

not exist at certain locations, and explained why higher values 

of Kv would improve the simulated fit to observed values. 

To improve simulation of hydraulic head in areas where an 

anticline is eroded, possibly exposing deeper interflow zones, 

DRN cells were added to the deeper model layers. Because the 

new drains were associated with river reaches, ELEVATION 

and CONDFACT for these new DRN cells were computed in 

the same way as described in section, “Streams and SurfaceWater Features.” 

Horizontal-flow barrier hydraulic characteristic was 

highly variable with a range of 1.0 × 10-15 to 3.88 × 10-5

(fig. 9; table 4). Early in the calibration process, results were 

not sensitive to the hydraulic characteristic assigned to the 

large offset faults at flow divides, so these values were fixed at 

a low value of 1.0 × 10-15.

Statistical Measures of Model Fit

The measure of model fit can be represented with many 

statistical and graphical methods. One measure of model fit 

is based on the difference between simulated and measured 

heads and flows, or residuals. The overall magnitude of 

the residuals was considered, but the distribution of those 

residuals, both statistically and spatially, could be equally 

important. The magnitude of residuals could initially point to 

gross errors in the model, the data (measured quantity), or how 

the measured quantity is simulated (Hill, 1998). A complete 

discussion of the statistical measures discussed in this section 

is available in Hill (1998).

Steady-State Model Fit and Model Error

The calibrated predevelopment steady-state model was 

a parsimonious model in which model simplification was 

made with relatively minor degradation to model fit. The 

distribution of Kh and Kv was considerably simpler than 

the distribution used in previous models (Davies-Smith and 

others, 1988; Hansen and others, 1994; Ely and others, 2011). 

This is attributed to two features of the new model: (1) the 

use of the HFB package precluded the need to simulate very 

narrow regions of low Kh (for example, see the permeability 

distributions of Hansen and others [1994]), and (2) the use of 

many thin layers to represent the CRBG units ensured proper 

local connectivity with hydraulic features, excluding the need 

to create complex Kh and Kv zones to explain local head 

anomalies that were the result of attempting to simulate fine 

scale hydraulic features with a coarse model grid.

A relatively good fit of the steady-state predevelopment 

model was achieved with a wide range of parameters. Many 

thin layers were used in model construction and a detailed 

geologic model ensured that the head-controlling hydraulic 

features (for example, surface-water features) were well 

represented. Continued automated calibration was determined 

to further minimize the calibration objective function, but 

deviate from the conceptual understanding of the flow system. 

This occurred because of calibration data bias, in particular:

• Most water-level measurements were made in 

commingling or commingling affected-wells, which 

resulted in highly variable head observations assigned 

to most of the constructed depth interval of wells. This 

indicated that the average vertical gradient is small, 

even though the high variability of measurements 

likely is caused by commingling of wells in 

heterogeneous geology (Burns and others, 2012b). 

Predevelopment vertical gradients were seldom, if 

ever, represented by the data. Because automated 

calibration will seek to match most of the data on 

average, Kv will be erroneously increased.

• In areas with simulated strong upward gradients, 

measured commingled heads would be persistently less 

than simulated head at the well bottom (where the head 

observation is assigned). Conversely, in areas of strong 

downward gradients, measured commingled heads 

would be persistently higher than simulated head at 

the well bottom. Use of the simulated value of head at 

the well bottom would cause regional bias in residuals 

that automated parameter estimation will seek to 

remove. Because commingled heads depend on aquifer 

heterogeneity and the combination of wells constructed 

in an area, there is no a-priori way to select the model 

layer to assign reliably to the measurement.

• Little of the data represented true predevelopment 

conditions. Weights were developed to minimize the 

effect of observations that likely are compromised by 

anthropogenic activities, but the large number of latetime data and the uneven coverage of data representing 

true predevelopment conditions will result in calibrated 

parameters that fit the biased measurements.
Numerical Simulation of Groundwater Flow 51

A modified automated calibration was used to ensure that 

the final parameters resulting from the predevelopment steadystate calibration were appropriate for use with the transient 

simulation. In particular, Kv of the CRBG units was decreased 

to approximately fit the few locations where vertical gradients 

were supported by the data. Predevelopment and transient 

simulations were run using progressively lower values of Kv, 

and the final selected value (1.00 × 10-4 ft/d) was the lowest 

value above which model convergence issues began to occur 

for the transient simulation; actual Kv may be somewhat 

lower. Nevertheless, this value of Kv provided reasonable 

model fit to the few data documenting high vertical gradients, 

and allowed testing of the transient model. Final steady-state 

model calibration was achieved by fixing CRBG Kv at this 

value, finding reasonable values of CRBG Kh, and allowing 

PEST to adjust HFB conductances to obtain estimates of 

model parameters for transient model testing.

Comparison of Simulated and Measured Hydraulic Heads

Unweighted hydraulic-head residuals (fig. 20) and a 

comparison between unweighted and weighted residuals 

(fig. 21) for the steady-state, predevelopment model 

indicates that most large-error simulated values occur at 

locations with lower observation weights. Low weights 

are used for observations that likely do not represent true 

predevelopment conditions (see figure 16 and the summary 

of the weighting strategy in section, “Water-Level Altitudes, 

Water-Level Altitude Changes, and Associated Error”). 

Observations with low weights are still valuable in the model 

calibration, and therefore are not removed. In areas with data 

ranging from high to low quality (for example, water levels 

likely affected by groundwater development), high-quality 

data will dominate the calibration. However, some areas of the 

model contain only low-quality data, and their removal from 

the calibration would leave the parameter estimation process 

completely unconstrained.

Following transient model calibration, the steady-state 

calibration was rechecked by comparing the magnitude and 

direction of the unweighted residuals to simulated drawdowns. 

Because the anticipated high-error low-weight residuals 

potentially represent late-time conditions, the unweighted 

steady-state residuals should be similar in magnitude and 

direction to the late-time drawdowns (2000–2007). For 

example, the Saddle Mountains Basalt unweighted residuals 

were similar in magnitude and direction when compared 

with the late-time simulated drawdown maps (fig. 22), 

indicating that the weighting scheme was appropriate and 

that the magnitude of declines was approximately correct. 

This comparison is not always straightforward (especially for 

thick sequences of hydrogeologic unit) because the computed 

drawdown maps show average values across many cells 

vertically, and because the measurements are affected by 

commingling. However, this method of comparison shows 

that unweighted residuals generally compare favorably with 

drawdown maps, indicating that the steady-state calibration 

method provided robust results. The generally good agreement 

between the location and magnitude of the residuals compared 

to simulated drawdown shows that steady-state calibration 

error can be because many measurements do not represent 

predevelopment conditions.

Comparison of Simulated and Measured Streamflow

The CPRAS model simulated streamflow for 1,727 river 

and 43,311 drain cells. The simulation of streamflow included 

only streamflow gains (river and drain cells) and losses (river 

cells only) and did not account for diversions and returns or 

any other aspects of the regulated system. This approach is 

valid for the predevelopment steady-state model presented 

here, but would not have accurately reproduced streamflow 

in the transient simulation. A comparison of simulated and 

estimated base flow at selected sites in the study area provided 

additional information on the reliability of the CPRAS model. 

Streamflow gains and losses also are important components of 

the simulated flow-system water budget, especially related to 

the total water mass balance.

Simulated annual base flows for predevelopment 

conditions generally were within the upper and lower bounds 

of estimated observations (table 5; fig. 17). Of the 50 base flow 

sites, 14 (28 percent) were not within the observations bounds, 

but the errors (simulated base flow minus lower/upper bound) 

generally were small. The largest error, Yakima River at 

Kiona, Washington (site 24), has a large amount of uncertainty 

associated with the base-flow estimate. The Yakima River and 

its major tributaries originate outside the model domain and 

some attempt was made to subtract that amount of inflow from 

the base-flow estimate. However, deriving annual groundwater 

contributions to streamflows downstream of the headwaters 

for predevelopment conditions introduced errors.

The method of constraining the parameter estimation 

with upper and lower bounds allowed for a wide margin of 

error in some cases. Therefore, this method was considered 

appropriate for the regional nature of the model, the lack 

of predevelopment and annual base-flow estimates, and the 

simplistic approach toward simulating a heavily regulated 

surface-water system.
52 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

tac13-0856_fig 20

EXPLANATION

Active model domain

Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Unweighted hydraulic-head residual, in feet

Less than -500

-499.99 to -400.00

-399.99 to -300.00

-299.99 to -200.00

-199.99 to -100.00

-99.99 to -10.00

-9.99 to 10.00

10.01 to 100.00

100.01 to 200.00

200.01 to 300.00

300.01 to 400.00

400.01 to 500.00

Greater than 500

Figure 20. Unweighted hydraulic-head residuals (simulated minus measured) from the steady-state model, Columbia Plateau 

Regional Aquifer System, Idaho, Oregon, and Washington.
Numerical Simulation of Groundwater Flow 53

tac13-0856_fig21

EXPLANATION

Unweighted residual

Weighted residual

Hydraulic-head residual, in feet -1,500

-1,000

-500

0

500

1,000

1,500

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Observation weight

Figure 21. Comparison of unweighted and weighted hydraulic-head residuals with weight, Columbia Plateau Regional 

Aquifer System, Idaho, Oregon, and Washington.
54 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

tac13-0856_fig 22

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

EXPLANATION

-9.9 to 10

Saddle Mountains unit simulated drawdown, in feet

Active model domain

-528.6 to -350

-349.9 to -250

-249.9 to -200

-199.9 to -150

-149.9 to -100

-99.9 to -50

-49.9 to -10

10.1 to 50

50.1 to 91.4

Saddle Mountains unit 

unweighted residuals, in feet

-398.4 to -250.0

-249.9 to -100.0

-99.9 to -10.0

-9.9 to 10.0

10.1 to 100.0

100.1 to 250.0

250.1 to 417.3

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

John Day

River Deschutes River

Palouse River

Yakima Rive

K

r

lickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 22. Comparison of the unweighted residuals with simulated current conditions (2000–2007) average drawdown for the 

Saddle Mountains unit, Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington.
Numerical Simulation of Groundwater Flow 55

Table 5. Simulated steady-state base flow and errors, Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington.

[Site locations are shown in figure 17. All values are in cubic feet per second. Station name: OR, Oregon; WA, Washington; ID, Idaho; NR, near; R, river, W, 

west; abv, above; BNDY, boundary; CR or C, creek; SO, south, FK, fork. Error: Calculated as simulated base flow minus lower/upper bound. If simulated base 

flow is within the lower and upper bounds, error is 0]

Base flow 

observation 

site

Station name Simulated

base flow

Lower

bound

Upper

bound Error

1 WHITE RIVER BELOW TYGH VALLEY, OR 26.0 118.9 378.6 -92.9

2 RHEA CREEK NEAR HEPPNER, OR 11.5 4.0 26.0 0.0

3 ROCK CREEK AB WHYTE PARK NR CONDON, OR 28.0 2.0 68.2 0.0

4 WILLOW CREEK ABV WILLOW CR LAKE, NR HEPPNER, OR 11.2 1.2 18.1 0.0

5 BUTTER CREEK NEAR PINE CITY, OR 51.9 5.2 36.9 15.0

6 MCKAY CREEK NEAR PILOT ROCK, OR 118.0 1.4 114.9 3.1

7 UMATILLA R AT W RESERVATION BNDY NR PENDLETON, OR 391.5 41.5 520.2 0.0

8 MOONSHINE CREEK NEAR MISSION, OR 4.8 0.1 2.9 1.9

9 UMATILLA RIVER AT PENDLETON, OR 438.9 45.6 483.7 0.0

10 UMATILLA RIVER ABOVE MEACHAM CREEK, NR GIBBON, OR 163.2 44.1 223.2 0.0

11 ROCK CREEK NEAR ROOSEVELT, WA 58.7 1.7 170.0 0.0

12 SOUTH FORK WALLA WALLA RIVER NEAR MILTON, OR 137.1 102.2 168.4 0.0

13 SPRING CREEK NEAR BLOCKHOUSE, WA 3.4 1.4 140.0 0.0

14 LITTLE KLICKITAT R NR WAHKIACUS, WA 117.5 10.3 1,030.0 0.0

15 ALDER CR AT ALDERDALE, WA 2.5 0.3 30.0 0.0

16 LITTLE KLICKITAT R NR GOLDENDALE, WA 29.3 3.5 350.0 0.0

17 MILL CREEK NEAR BLOCKHOUSE, WA 12.3 1.4 140.0 0.0

18 MILL CREEK NEAR WALLA WALLA, WA 74.0 6.9 690.0 0.0

19 WALLA WALLA RIVER NEAR TOUCHET, WA 784.7 18.4 558.1 226.6

20 BLUE CREEK NEAR WALLA WALLA, WA 6.9 1.0 100.0 0.0

21 TOUCHET R NR TOUCHET, WA 281.0 14.8 1,480.0 0.0

22 MILL CREEK NEAR WALLA WALLA, WA 87.1 1.1 77.4 9.7

23 DRY CREEK NEAR WALLA WALLA, WA 17.6 1.3 130.0 0.0

24 YAKIMA RIVER AT KIONA, WA 1,003.8 1,489.4 3,189.5 -485.6

25 TOUCHET RIVER AT BOLLES, WA 250.7 14.5 1,450.0 0.0

26 EAST FK TOUCHET R NR DAYTON, WA 119.9 9.5 950.0 0.0

27 ASOTIN CREEK BELOW CONFLUENCE NEAR ASOTIN, WA 13.2 23.9 48.7 -10.7

28 GRANGER DRAIN AT GRANGER, WA 0.0 19.5 34.9 -19.5

29 TUCANNON RIVER NEAR STARBUCK, WA 222.2 12.7 1,270.0 0.0

30 AHTANUM CREEK AT UNION GAP, WA 78.4 15.4 68.4 10.0

31 MEADOW CREEK NR CENTRAL FERRY, WA 22.1 0.2 19.0 3.1

32 SO FK PALOUSE R ABV PARADISE C NR PULLMAN, WA 11.0 0.7 72.0 0.0

33 PARADISE CR NR PULLMAN, WA 7.2 0.5 51.0 0.0

34 PARADISE CR AT UNIVERSITY OF IDAHO AT MOSCOW, ID 1.8 0.8 7.0 0.0

35 SOUTH FORK PALOUSE RIVER AT PULLMAN, WA 20.1 3.8 29.2 0.0

36 MISSOURI FLAT CREEK AT PULLMAN, WA 1.6 0.4 38.0 0.0

37 PALOUSE RIVER AT HOOPER, WA 471.0 37.9 538.0 0.0

38 UNION FLAT CREEK NEAR COLFAX, WA 36.4 1.8 180.0 0.0

39 FOURMILE CR AT SHAWNEE, WA 0.8 0.5 48.0 0.0

40 PALOUSE RIVER BELOW SOUTH FORK AT COLFAX, WA 93.6 21.9 341.3 0.0

41 PALOUSE RIVER NR POTLATCH, ID 25.6 11.2 246.6 0.0

42 PALOUSE RIVER NEAR COLFAX, WA 63.3 17.5 1,750.0 0.0

43 PINE CREEK AT PINE CITY, WA 53.4 2.0 200.0 0.0

44 CRAB CREEK AT ROCKY FORD ROAD NEAR RITZVILLE, WA 16.9 11.9 35.6 0.0

45 ROCKY FORD CREEK NEAR EPHRATA, WA 2.4 6.9 690.0 -4.5

46 CRAB CREEK AT IRBY, WA 53.7 4.9 33.6 20.1

47 CRAB CREEK AT IRBY, WA 57.9 2.9 290.0 0.0

48 COAL CREEK AT MOHLER, WA 1.5 0.1 12.0 0.0

49 PARK CREEK BLW PARK LAKE NR COULEE CITY, WA 11.8 1.0 95.0 0.0

50 DOUGLAS CREEK NEAR ALSTOWN, WA 13.0 0.1 9.2 3.8
56 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Transient Model Fit and Model Error

Comparison of Simulated and Measured Hydraulic Heads

A traditional and intuitive assessment of model 

calibration is a simple plot of measured hydraulic heads as a 

function of simulated hydraulic heads (fig. 23). For 46,460 

water-level measurement points, the mean and median 

difference between simulated minus measured hydraulic heads 

is -10 and 4 ft, respectively. The residuals for the transient 

simulation period show that 52 percent of the simulated heads 

exceeded measured heads with a median residual value of 

43 ft, and 48 percent were less than measured heads with 

a median residual value of -76 ft. The residuals should be 

normally distributed along a line with a 1:1 line (slope equal 

to 1.0 and a y-intercept of zero), if no model bias exists. The 

weighted measurements compared with weighted simulated 

values generally are along a straight line with a slope of 1.03 

and a y-intercept of -18.0 ft.

The root-mean-square (RMS) error of the difference 

between simulated and measured hydraulic heads in the 

observation wells, divided by the total difference in measured 

hydraulic heads in the groundwater system (Anderson and 

Woessner, 1992, p. 241), should be less than 10 percent to be 

acceptable (Drost and others, 1999). The calibrated transient 

model produced an RMS error of 167, which divided by the 

total difference in water levels (5,648 ft), was 2.9 percent.

The residuals for overburden and basalt units indicated 

a reasonable fit (66 percent of the residuals were within 

±100 ft); however, the spatial distribution of the residuals less 

than -100 ft or greater than 100 ft showed definite patterns 

of bias (fig. 24). For example, 67 percent of simulated 

heads in the GWMA were less than measured heads, in 187 

instances at 24 locations, by greater than 500 ft. Most of the 

largest GWMA residuals were on the margins of the GWMA 

boundary or in the northern part of the region.

Simulated heads generally were larger than measured 

heads in the Yakima River Basin and Umatilla subarea, 63 

and 67 percent of the time, respectively. There also was a bias 

toward residuals of greater than 500 ft in areas of the Yakima 

River Basin and Umatilla subarea.

The pattern of underprediction in the GWMA and 

overprediction in the Yakima River Basin and Umatilla 

subarea was somewhat expected. In areas of significant 

upward hydraulic gradient, the simulated head likely would be 

higher than the measured commingled head, and conversely, 

in areas of downward gradient. The Yakima River Basin and 

Umatilla subarea have areas of upward gradients, and the 

GWMA has areas of downward gradients for approximately 

the upper 1,000 ft, although this pattern varies by position.

Model observations were assigned to one layer 

represented by the well depth reported on drillers’ logs. Many, 

if not most, of the wells used for this analysis are open to 

multiple transmissive interflow zones; therefore, the measured 

water level is a composite hydraulic head from across the 

open well interval, so water-level residuals (simulated minus 

measured) can be misleading. Burns and others (2012b) 

determined that in some areas of the Palouse Slope, hydraulic 

heads measured in wells open to the same formations at 

similar altitudes were highly variable (differing by hundreds of 

feet). The variability was attributed to strong vertical hydraulic 

gradients and geologic or well-construction variability 

causing different amounts of commingling, with heads in 

well-connected and more transmissive aquifers dominating 

individual boreholes. Evaluation of residuals shows bias 

associated with the use of commingled wells as observations. 

A more detailed analysis of residuals indicates that model fit 

likely would be better if all measurements were from wells 

constructed only in single aquifers.

Another way to assess the ability of the CPRAS model 

to represent groundwater flow and trends is to examine the 

vertical distribution of hydraulic heads in the Umatilla subarea 

and Palouse Slope well groups identified in Burns and others 

(2012b). The model layers containing simulated heads selected 

for plotting were identified as potentially commingled by 

using well construction records.

A total of 4,235 measured water levels from 286 wells in 

the Umatilla River Basin in the CRBG aquifers were divided 

into clustered groups of wells with similar water levels and 

trends (fig. 25). The method and justification for well zonation 

is described in Burns and others (2012b). Zones of low 

permeability may separate the groups of wells with similar 

water levels and trends. These zones represent leaky barriers 

to groundwater flow and compartmentalize the CRBG aquifer 

system. The degree of compartmentalization is variable, but 

it occurs in both the vertical and horizontal directions in the 

Umatilla area. Within each well group in the Umatilla River 

Basin, hydraulic heads can be correlated over tens of miles in 

the horizontal direction.

The transient simulation results showed a reasonably 

good match to measured water levels in the Umatilla area 

(fig. 26). The presence of large vertical gradients within 

groups can explain the large range of data values within 

many of the groups. For example, for Umatilla well group 2, 

simulated heads in layers 32–41 span the range of measured 

water levels. The timing of simulated declines agreed with 

measured water-level trends, with steeper declines occurring 

during the mid-1970s. Head differences between groups 

were preserved, ensuring that regional flow patterns were 

adequately simulated.
Numerical Simulation of Groundwater Flow 57

tac13-0856_fig23

0

1,000

2,000

3,000

4,000

5,000

6,000

7,000

0 1,000 2,000 3,000 4,000 5,000 6,000 7,000

Measured hydraulic head, in feet

Simulated hydraulic head, in feet

1:1 line

EXPLANATION

Grande Ronde

Wanapum

Saddle Mountains

Overburden

Figure 23. Measured hydraulic heads as a function of simulated heads, Columbia Plateau Regional Aquifer System, Idaho, 

Oregon, and Washington.
58 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

tac13-0856_fig 24

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

EXPLANATION

Greater than 500 feet

Residual—Difference between simulated 

 and measured water levels

Greater than 100 feet

Less than -100 feet

Less than -500 feet

Umatilla Subarea

Yakima River Basin

Ground Water Management Area

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 24. Differences between simulated and measured water levels (residuals), Columbia Plateau Regional Aquifer System, 

Idaho, Oregon, and Washington. The spatial distribution of residuals less than -100 feet or greater than 100 feet showed definite 

patterns of bias.
Numerical Simulation of Groundwater Flow 59

tac13-0856_fig 25

COLUMBIA RIVER

Umatilla River

Ione

Echo

Heppner

Umatilla

Boardman

Stanfield

Pendleton

Lexington

Hermiston

Pilot

Rock

Irrigon

UMATILLA

MORROW

YAKIMA

BENTON

KLICKITAT

WALLA

WALLA

WILLOW CREEK MONOCLINE

REITH

ANTICLINE

SERVICE ANTICLINE

AGENCY SYNCLINE

DALLES - UMATILLA SYNCLINE

Study

area

WASHINGTON

OREGON

IDAHO

EXPLANATION

Wells by group

Group 2

Group 6

Group 9

Group 1

Group 3

Group 4

Group 5

Group 7

Group 8

Group 12

Group 13

Group 11

Group 10

Group 14

Group 15

Group 16

Group 17

Group 18

Shallow-deep pairs

Fault

Fold

Yellow pair

Blue pair

Purple pair

Orange pair

Umatilla Subarea

Oregon Water Resources Department

Administrative Areas

Butter Creek

Ella Butte

Ordnance Basalt

Stage Gulch

0 10 20 MILES

0 10 20 KILOMETERS

Base map modified from USGS and other digital data, various 

scales. Coordinate system: State Plane-Washington South FIPS 

4602, Projection: Lambert Conformal Conic; North American 

Datum of 1983.

119°

46°

Figure 25. Well groups with similar hydraulic response within the Umatilla subarea, Oregon. Wells are grouped to show a general 

North–South transect (circles and squares) and a general East–West transect (triangles and squares). From Burns and others 

(2012b).
60 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Umatilla 7-17

Umatilla 7-15

Umatilla 14-31

Umatilla 14-36

Umatilla 3-18 Umatilla 3-17

Umatilla 10-33 (East)

Umatilla 10-35 (West)

Umatilla 10-39 (East)

Umatilla 10-39 (West)

Umatilla 4-40 (East)

Umatilla 4-41 (West)

Umatilla 4-42(West)

Umatilla 4-38(East)

Umatilla 5-36

Umatilla 4-37(East)

Umatilla 15-42

Umatilla 15-47

Umatilla 16-39

Umatilla 17-50

Umatilla 18-41

Umatilla 8-38

Umatilla 8-40

Umatilla 5

EXPLANATION

Umatilla 10

Umatilla 3

Umatilla 4

Umatilla 7

Umatilla 8

Umatilla 14

Umatilla 15

Umatilla 16

Umatilla 17

Umatilla 18 Umatilla 8-40

Umatilla group Model layer

Measured water level

Simulated hydrograph 200

400

600

800

1,000

1,200

1,400

1,600

1,800

2,000

1935 1940 1945 1950 1955 1960 1965 1970 1975

Year

1980 1985 1990 1995 2000 2005 2010

Groundwater level, in feet above North American Vertical Datum of 1988

0

A. North-South group

Figure 26. Simulated water levels for the (A) North-South; and (B) East-West well groups near the Oregon Water Resources 

Department administrative areas in the Umatilla subarea, Oregon. Measured water levels depict the winter median water 

level for individual wells within each well group. Locations of wells are shown in figure 25. The numbered lines corresponding 

to each data series show transient simulated water levels from layers intersected by nearby potentially commingling wells.
Numerical Simulation of Groundwater Flow 61

tac13-0856_fig26b

Umatilla 13

Umatilla 2

EXPLANATION

Umatilla 3

Umatilla 6

Umatilla 7

Umatilla 8

Umatilla 9

Umatilla 11

Umatilla 12

Umatilla 14

Umatilla 9-38

Umatilla group Model layer

Simulated hydrograph

Groundwater level, in feet above North American Vertical Datum of 1988

1935 1940 1945 1950 1955 1960 1965 1970 1975

Year

1980 1985 1990 1995 2000 2005 2010

100

200

300

400

500

600

700

800

900

1,000

B. East-West group

Umatilla 11-17

Umatilla 2-41

Umatilla 2-39

Umatilla 2-32

Umatilla 3-17

Umatilla 3-18

Umatilla 6-16

Umatilla 7-15

Umatilla 7-16

Umatilla 8-40

Umatilla 8-38

Umatilla 9-38

Umatilla 9-40

Umatilla 9-42

Umatilla 12-34

Umatilla 12-32

LAY 34

Umatilla 13-36

Umatilla 13-35

Umatilla 13-32

Umatilla 14-36

Umatilla 14-29

Umatilla 14-31

Figure 26.—Continued
62 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

During model calibration, it was determined that wells 

near groups 4, 8, and 10 had steeper water-level declines 

than simulated declines, so modifications to the model were 

made to better match the observations. Two mechanisms 

were identified as possible causes: (1) flow barriers that 

exist in shallow CRBG units were not present in deep units, 

allowing commingling to drain water away through the deeper 

aquifer system, or (2) aquifer storage is much lower in some 

mid-slope areas, causing much steeper water-level declines 

in response to the relatively modest amount of simulated 

pumping. To evaluate the efficacy of the first mechanism, the 

final model was calibrated by removing horizontal-flow barrier 

segments 61, 63, and 65 below model layer 40 and horizontalflow barrier segments 70 and 71 from all CRBG units except 

Saddle Mountains (fig. 9). This mechanism was deemed the 

more likely cause of the observations, because (1) the final 

heads in several different groups apparently were trending 

toward the same final head value and (2) storage parameter 

adjustments had a negligible effect. Mechanistically, the flow 

barriers in the shallower units are attributed to depositional 

variability that occurred as the younger lava flows onlapped 

the anticline that was contemporaneously forming the 

Blue Mountains.

The GWMA encompasses much of the Palouse Slope 

and the eastern part of the Yakima Fold Belt, which forms 

a transition area between the two physiographic provinces 

(fig. 1). Burns and others (2012b) examined 8,622 measured 

water levels from 1,202 wells within the Palouse Slope and 

eastern Yakima Fold Belt in the CRBG aquifers and divided 

them into groups of wells exhibiting similar changes in 

hydraulic head over time (fig. 27).

Similar to the analysis of the Umatilla subarea well 

groups, a comparison of measured and simulated heads 

by group was performed for the Palouse Slope. Unlike the 

Umatilla subarea, discrete horizontal-flow barriers were 

much less influential in the Palouse Slope area, consistent 

with the lower degree of structural deformation in this gently 

folded area. The lack of horizontal-flow barriers explains 

the significant overlap of measured heads between groups, 

with the large range in measured heads attributed to using 

commingling water levels in an area with large vertical 

gradients (fig. 28).

The magnitude, timing, and range of simulated heads 

were in good general agreement with measured water levels. 

Comparison between the groups of wells representing different 

flow paths (figs. 28A–D), indicated that the simulated regional 

flow field matched the measured flow field, including the 

reversal of groundwater-flow direction in the shallow aquifers 

that accompanied application of large quantities of CBIP 

irrigation water (Burns and others, 2012b).

Comparison of Simulated and Measured Streamflow

Average long-term base-flow estimates of streamflow 

were helpful in calibrating the steady-state model during 

automated calibration, but uncertainties in those estimates 

made them insufficient to fully evaluate the transient 

model. The addition of large-scale agricultural irrigation 

and groundwater pumping has had competing effects on 

streamflow in the CPRAS. Surface-water irrigation, and to 

a lesser degree groundwater irrigation added large amounts 

of water to shallow groundwater storage throughout 

the period of development. This additional water has 

increased streamflow in areas of irrigation, even resulting 

in groundwater flooding and landslides. However, in areas 

away from major surface-water irrigation, pumping has 

reduced groundwater in storage, lowering water levels 

and streamflows. These patterns and trends were largely 

reproduced by the transient model, as shown by the change in 

base flows from predevelopment to current conditions (2000–

2007; fig. 29). The Yakima River Basin and the Columbia 

River along the eastern boundary of the Yakima River Basin 

consistently showed an increase in base flow compared 

to predevelopment conditions. Groundwater pumping has 

captured and thereby reduced annual streamflow in the 

Yakima River system by about 200 ft3/s (Ely and others, 

2011), but that reduction is small compared to increases 

from irrigation, so there is a net gain in flow. Surface-water 

irrigation in the Yakima River Basin has been extensive, but 

simulating the Yakima River headwaters and the basin-wide 

regulation (reservoirs, diversions, and returns) was beyond 

the scope of this study.

The vast majority of stream reaches in the CPRAS 

model area have decreased simulated base flows compared 

to predevelopment conditions (fig. 29). This decrease is 

caused by two factors. First, predevelopment recharge from 

precipitation is greater than current condition recharge from 

precipitation. Second, groundwater pumping has reduced 

water levels and captured streamflow. One area of interest is 

Crab Creek (sites 44, 46, and 47, fig. 17). Base flow along 

the upper reaches of Crab Creek has decreased over recent 

time, whereas the model simulated modest increases. The 

model did simulate the decreased flow in the lower reaches 

of Crab Creek. This inconsistency shows the difficulty of 

simulating site-specific stresses in a regional model with 

coarse horizontal discretization.
Numerical Simulation of Groundwater Flow 63

tac13-0856_fig 27

Pasco

FERRY STEVENS

OKANOGAN

CHELAN

DOUGLAS

SPOKANE GARFIELD

LINCOLN

GRANT

KITTITAS

WHITMAN

ADAMS

YAKIMA

FRANKLIN

BENTON

COLUMBIA

FRENCHMAN HILLS

SADDLE MOUNTAINS

Moses

Lake

WALLA

WALLA

COLUMBIA

Banks Lake

RIVER

East Low Canal

McCarteney

Creek

West Canal

Lower Crab Creek

COLUMBIA

Palouse River

Snake River

Yakima River

Potholes

Reservoir

Crab Creek

EASTERN

SOUTHERN

WESTERN

120° 119°30' 119° 118°30' 118°

48°

47°

30'

47°

46°

30'

0 5 10 15 20 25 MILES

0 5 10 15 20 25 KILOMETERS

Base map modified from USGS and other digital data, various 

scales. Coordinate system: State Plane-Washington South FIPS 

4602, Projection: Lambert Conformal Conic; North American 

Datum of 1983.

EXPLANATION EXPLANATION

508-14

Odessa

Quincy

Geologic model extent of Columbia River

Basalt Group

Columbia Basin Ground Water Management Area

Washington Department of Ecology Ground 

Water Management Subareas 

Fault

Fold

Generalized groundwater 

flow path under current 

conditions

Well

Group 6

Group 7

Group 5

Middle groups

Group 13

Group 14 - 1

Group 14 - 2

Group 15

Southern groups

Group 11 - Deep

Group 11 - Shallow

Group 8

Group 10

Group 9

Group 12

Eastern groups

Group 2

Group 1

Group 3

Group 4 - Deep

Group 4 - Shallow

Western groups

Palouse Slope

Clearwater Embayment

Blue Mountains

Yakima Fold Belt

B

C

P

Y

Structural regions

OR ID

WA

B

C

P

Y

Figure 27. Well groups with similar hydrologic responses and generalized groundwater-flow paths under 2000–10 conditions near 

the Washington State Department of Ecology administrative areas in parts of the Palouse Slope/eastern Yakima Fold Belt and the 

Columbia Basin Ground Water Management Area, Washington. From Burns and others (2012b).
64 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

tac13-0856_fig28a

Palouse 1

Layer 34

Layer 36

Palouse 2

Layer 31

Layer 33

Layer 34

Palouse 3

Layer 18

Layer 34

Palouse 4 Shallow

Palouse 4 Deep

Layer 29

Layer 31

Layer 36

Well group and layer No.

EXPLANATION

Year

0

200

400

600

800

1,000

1,200

1,400

1,600

1,800

2,000

1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000 2010

Water-level altitude, in feet above North American Vertical Datum of 1988

A. Western flow path

Figure 28. Simulated water levels for the (A) western flow path; (B) eastern flow path; (C) middle flow path; and 

(D) southern flow path well groups near the Washington State Department of Ecology administrative areas in parts of the 

Palouse Slope/eastern Yakima Fold Belt and the Columbia Basin Ground Water Management Area, Washington. Locations 

of wells are shown in figure 27. The lines under each well group show corresponding transient model simulation results from 

layers intersected by nearby potentially commingling wells.
Numerical Simulation of Groundwater Flow 65

tac13-0856_fig28b

Well group and layer No.

EXPLANATION

Palouse 8

Layer 32

Layer 33

Layer 34

Layer 42

Palouse 9

Layer 31

Layer 34

Layer 40

Palouse 10

Layer 31

Layer 34

Layer 40

Palouse 11 Shallow

Layer 31

Layer 33

Layer 36

Palouse 11 Deep

Layer 41

Layer 42

Palouse 12

Layer 31

Layer 44

Layer 55

Year

0

200

400

600

800

1,000

1,200

1,400

1,600

1,800

2,000

1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000 2010

Water-level altitude, in feet above North American Vertical Datum of 1988

B. Eastern flow path

Figure 28.—Continued
66 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

tac13-0856_fig28c

Well group and layer No.

EXPLANATION

Palouse 5

Layer 30

Layer 32

Layer 34

Layer 39

Palouse 6

Layer 29

Layer 31

Layer 34

Palouse 7

Layer 31

Layer 33

Layer 36

Year

0

200

400

600

800

1,000

1,200

1,400

1,600

1,800

2,000

1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000 2010

Groundwater altitude, in feet above North American Vertical Datum of 1988

C. Middle flow path

Figure 28.—Continued
Numerical Simulation of Groundwater Flow 67

tac13-0856_fig28d

Well group and layer No.

EXPLANATION

Layer 31

Layer 55

Layer 44

Layer 39

Layer 45

Layer 36

Layer 31

Layer 30

Palouse 12

Palouse 13

Palouse 14.1

Palouse 14.2

Palouse 15

Layer 24

Layer 17

Year

D. Southern flow path

0

200

400

600

800

1,000

1,200

1,400

1,600

1,800

2,000

1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000 2010

Groundwater altitude, in feet above North A

merican Vertical Datum of 1988

Figure 28.—Continued
68 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

tac13-0856_fig 29

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

EXPLANATION Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

Change in simulated base flow, in cubic feet per second

-144.95 to -100.00

-99.99 to -10.00

-9.99 to -1.00

-0.99 to -0.10

-0.09 to 0.10

0.10 to 1.00

1.01 to 10.00

10.01 to 100.00

100.01 to 615.04 No

symbol

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane R

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 29. Change in simulated base flow from predevelopment (pre-1900) to current conditions (2000–2007) for the Columbia 

Plateau Regional Aquifer System, Idaho, Oregon, and Washington. Red indicates current condition base flow is less than 

predevelopment base flow. Blue indicates current condition base flow is greater than predevelopment base flow.
Numerical Simulation of Groundwater Flow 69

Simulated Potentiometric Surfaces

Each layer of the model can be conceptually viewed 

as a (potential) aquifer. The presence of an actual aquifer is 

controlled by the presence of sufficient local porosity and 

permeability, which varies as a function of geologic deposition. 

Previous modeling studies simulated significantly fewer 

model layers. To allow comparison with previous studies 

and create summary maps of conditions, composite head and 

drawdowns were computed (figs. 30 and 31, respectively) for 

the Overburden, Saddle Mountains, and Wanapum units, and 

the upper 10 model layers of Grande Ronde Basalt (as much 

as 2,000 ft near the center of the study area). Composite heads 

were computed as the thickness-weighted average of heads in 

the layers being averaged. Dry cells (simulated head below 

the cell bottom) were not used in the computation. Composite 

drawdown was computed as predevelopment composite head 

minus the late-time (2000–2007) average of the composite 

heads for each year. In all cases where cells were dry for only a 

few of the years, the bottom elevation of the lowest model layer 

used in the average was used for the dry years. If all cells were 

dry for all years, the final drawdown value was 0.

Model Uncertainty and Limitations

The CPRAS model is a set of mathematical equations 

designed to represent an extremely complex natural system. 

Furthermore, this natural system has been perturbed by 

human activities in ways only generally understood and 

not fully quantified. Intrinsic to the model is the error and 

uncertainty associated with the approximations, assumptions, 

and simplifications that must be made. In addition to those 

intrinsic errors, hydrologic modeling errors typically are the 

consequence of a combination of errors in the (1) input data, (2) 

representation of the physical processes by the algorithms of 

the model, and (3) parameter estimation during the calibration 

procedure (Troutman, 1985). These three types of model errors 

limit application of the CPRAS model as follows:

• Data on extents and thicknesses of mapped 

hydrogeologic units, location and nature of structural 

features, water levels, pumpage, recharge, and hydraulic 

properties were taken from Kahle and others (2009), 

Burns and others (2011, 2012b), and other previous 

investigations. Most of the data were concentrated along 

major river valleys and populated areas. This means 

that for some of the study area, information was not 

available to calibrate the model, especially for the areas 

lacking water-level data. Additionally, the methods used 

to estimate important stresses of the flow system, such 

as groundwater pumping and recharge, were limited to 

recent years. Reconstructing model stresses for the first 

half of the model simulation period likely led to larger 

error and uncertainty.

• A numerical model cannot completely represent all 

physical processes within a flow system. Determining 

if a weakness in a simulation is attributable simply 

to input data error or to shortcomings in how the 

model represents the governing physical processes 

is often not possible. The model inevitably relies on 

simplifying assumptions and generalizations that affect 

the results of the simulation in complex ways. The 

CPRAS model was not designed to represent some 

details of the hydrologic processes. For example, the 

extensive network of surface-water diversions, canals, 

wasteways, and drains were not directly simulated. 

Model drain cells were placed in any model cell that 

might have streamflow based on the NHD (Simley 

and Carswell, 2009). Reservoirs and reservoir releases 

on the Columbia River and its tributaries were not 

included in the CPRAS model. Furthermore, modeldiscretization errors resulted from (1) the effects of 

averaging altitude information over the model cell size, 

(2) the time averaging of observed values inherent 

in an annual simulation stress period, and (3) the 

inaccuracies in translating mapped hydrogeologic units 

into orthogonal model cells.

• Errors in parameter estimates occur when improper 

values are selected during the calibration process. 

Various combinations of parameter values can result in 

low residual error, yet improperly represent the natural 

system. An acceptable degree of agreement between 

simulated and measured values does not guarantee 

that the estimated model parameter values uniquely 

and reasonably represent the actual parameter values. 

The use of automated parameter estimation techniques 

and associated statistics, such as composite scaled 

sensitivities and correlation coefficients, removes 

(or allows consideration of) some of the effects of 

non-uniqueness, but certainly does not eliminate the 

problem entirely. Ensuring that calibrated values are 

comparable to a reasonable range of independentlyderived or literature values also can reduce error 

caused by parameter estimation. Limitations of the 

model observations used in this study include errors in 

simulated heads for multilayer wells, uncertainty in the 

model hydrogeologic unit designations and depths, and 

streamflow observations based on base-flow estimates 

from different sources and methods.

If the coarse horizontal discretization, annual stress 

periods, and estimated stresses are considered, the effects of 

the simplifications and other potential errors can be limited. 

If the model is used for simulations beyond which it was 

designed, however, the generalizations and assumptions used 

could significantly affect the results.
70 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

tac13-0856_fig 30a

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

EXPLANATION Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

Active model domain

Overburden simulated hydraulic head, in feet

Dry

501 to 1,000

1,001 to 1,500

1,501 to 2,000

2,001 to 2,500

2,501 to 3,000

3,001 to 3,500

0 to 500 3,501 to 4,000

4,001 to 4,500

4,501 to 5,000

5,001 to 5,500

5,501 to 5,952

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

A

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 30. Simulated composite hydraulic head for (A) Overburden unit, (B) Saddle Mountains unit, (C) Wanapum unit, and 

(D) upper Grande Ronde unit for the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington.
Numerical Simulation of Groundwater Flow 71

tac13-0856_fig 30b

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

EXPLANATION Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

Active model domain

Saddle Mountains simulated hydraulic head, in feet

Dry

501 to 1,000

1,001 to 1,500

1,501 to 2,000

2,001 to 2,500

2,501 to 3,000

3,001 to 3,500

0 to 500 3,501 to 3,786

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

B

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 30.—Continued
72 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

tac13-0856_fig 30c

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

EXPLANATION Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

Wanapum simulated hydraulic head, in feet

Active model domain

Dry

501 to 1,000

1,001 to 1,500

1,501 to 2,000

2,001 to 2,500

2,501 to 3,000

3,001 to 3,500

0 to 500 3,501 to 4,000

4,001 to 4,500

4,501 to 5,000

5,001 to 5,500

5,501 to 5,734

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

C

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 30.—Continued
Numerical Simulation of Groundwater Flow 73

tac13-0856_fig 30d

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

EXPLANATION Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

Active model domain

Upper Grande Ronde simulated hydraulic head, in feet

0 to 500

501 to 1,000

1,001 to 1,500

1,501 to 2,000

2,001 to 2,500

2,501 to 3,000

Dry 3,001 to 3,500

3,501 to 4,000

4,001 to 4,500

4,501 to 5,000

5,001 to 5,500

5,501 to 6,000

6,001 to 6,535

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

D

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 30.—Continued
74 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

tac13-0856_fig 31c tac13-0856_fig 31a

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

EXPLANATION Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 10.1 to 50

50.1 to 113.5

-504 to -300

-299.9 to -250

-249.9 to -200

-199.9 to -150

-149.9 to -100

-99.9 to -50

-49.9 to -10

No

symbol -9.9 to 10

Overburden simulated drawdown, in feet

Active model domain

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

A

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 31. Simulated composite drawdowns from predevelopment to current conditions (2000–2007) for (A) Overburden unit, 

(B) Saddle Mountains unit, (C) Wanapum unit, and (D) upper Grande Ronde unit for the Columbia Plateau Regional Aquifer 

System, Idaho, Oregon, and Washington.
Numerical Simulation of Groundwater Flow 75

tac13-0856_fig 31c tac13-0856_fig 31b

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

EXPLANATION

Saddle Mountains simulated drawdown, in feet

Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. -528.6 to -350

-349.9 to -250

-249.9 to -200

-199.9 to -150

-149.9 to -100

-99.9 to -50

-49.9 to -10

10.1 to 50

50.1 to 91.4

No

symbol -9.9 to 10

Active model domain

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

B

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 31.—Continued
76 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

EXPLANATION

10.1 to 50

50.1 to 100

100.1 to 150

150.1 to 200

200.1 to 250

250.1 to 300

300.1 to 350

350.1 to 400

400.1 to 450

-204.3 to -150

-149.9 to -100

-99.9 to -50

-49.9 to -10

No

symbol -9.9 to 10

Wanapum simulated drawdown, in feet

Active model domain

Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

C

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 31.—Continued
Numerical Simulation of Groundwater Flow 77

tac13-0856_fig 31d

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Moses

Lake

Walla

Walla

EXPLANATION

Upper Grande Ronde simulated drawdown, in feet

Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. -361 to -150

-149.9 to -100

-99.9 to -50

-49.9 to -10

No

symbol -9.9 to 10

10.1 to 50

50.1 to 100

100.1 to 150

150.1 to 200

200.1 to 250

250.1 to 300

300.1 to 532.3

Active model domain

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

D

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 31.—Continued
78 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Model Application

The CPRAS model was constructed to simulate regionalscale groundwater flow and it can be used to help answer 

questions regarding groundwater-flow issues at that scale. For 

example, the interactions of stresses, such as recharge and 

groundwater pumping, on the groundwater or surface-water 

system can be evaluated for major basins or subbasins, such as 

the Yakima and Umatilla River Basins or the Palouse Slope. 

The CPRAS model includes some localized features and sitespecific data, where available, but it is not meant to evaluate 

issues or hydrologic processes at local scales.

The CPRAS model also can be used to evaluate 

alternative conceptualizations of the flow system that are 

likely to have a regional effect. These might include the effects 

of climate, different interpretations of geologic structures, or 

the role of commingling wells on hydraulic heads.

Many aspects of the CPRAS affect the groundwaterflow system, including groundwater pumping, irrigation 

application, geologic structure (faults and folds), and the 

complex, heterogeneous hydraulic properties of the CRBG and 

Overburden unit. Throughout some of the study area, data are 

limited and the CPRAS model has been the appropriate tool to 

understand the role of these aspects of the system over time.

The annual stress period used in the model was selected 

to provide adequate temporal resolution for analyses of 

variations in pumping and recharge rates corresponding to 

annual changes in precipitation, pumping, and irrigation. 

Effects of management decisions on time scales less than 

multi-year are unlikely to be adequately simulated with the 

CPRAS model; however, the model can be used appropriately 

to analyze long-term changes in water-use practices or 

potential future climate.

The basic structure of the calibrated model allows for 

alternative uses of the model. Sensitivities of stresses on the 

flow system can provide information for directing additional 

data analyses and (or) data collection. Cause and effect also 

can be assessed. For example, assessing the dual effects of 

increased irrigation recharge and groundwater development on 

the water levels would be an appropriate use of the model.

The model can be used to examine the effects of 

continued or increased pumping on the regional groundwaterflow system to evaluate the efficacy of various groundwater 

resources management alternatives. For increased pumping, 

the model should not be used for assessing the effects of 

one or two new wells on the flow system, but could be used 

for such applications as estimating the quantity of pumpage 

in an area that leads to specified water-level declines and 

(or) streamflow capture. With increasing demand for water 

throughout the Columbia Plateau, the CPRAS model could 

be further developed to test optimization strategies for the 

conjunctive use of surface water and groundwater given 

specified constraints.

Groundwater Budget

The CPRAS model was used to derive components 

of the groundwater budget for the simulation period of 

predevelopment to 2007. During this period, the distribution 

and amount of pumpage changed and recharge varied, and a 

cumulative or mean annual water budget would not highlight 

these variations. Thus, simulated annual water budgets are 

presented for predevelopment to 2007 (fig. 32).

Water-budget components are presented as the difference 

between inflows and outflows to the groundwater system. 

Inflows are fluxes into the aquifer system. Outflows are fluxes 

out of the aquifer system. For example, recharge is considered 

an inflow and net well withdrawals are considered an outflow. 

Discharge from the groundwater system to streamflow (rivers 

and drains) is presented as an outflow. Storage is shown so 

negative numbers represent outflows from the groundwater 

system into storage, so a net gain in storage is represented by 

negative numbers.

The annual groundwater budgets show several trends 

over the simulation period (fig. 32). Groundwater pumping 

was negligible until the 1950s and began to increase 

significantly during the 1970s and 1980s. Recharge was highly 

variable due to the interannual variability of precipitation, but 

began to increase in the late 1940s because of the increase 

in irrigation projects. Streamflow gains and net storage 

followed recharge closely, but a loss of groundwater in 

storage increased and streamflow decreased as groundwater 

pumping increased.

Included in the simulation period were examples of a 

wet (1997), average (2000), and dry (2001) year (fig. 33). 
Numerical Simulation of Groundwater Flow 79

tac13-0856_fig32

Storage

Streams

Recharge

Wells

Year

1920 1930 1940 1950 1960 1970 1980 1990 2000 2010

Outflows Inflows -10

Predevelopment

12

-8

-6

-4

-2

0

2

4

6

8

10

Water budget flux, in millions of acre-feet per year

Figure 32. Simulated annual water-budget flux for predevelopment to 2007 for the Columbia Plateau Regional Aquifer 

System, Idaho, Oregon, and Washington.

Storage

Streams

Recharge

Wells

EXPLANATION

-10

-8

-6

-4

-2

0

2

4

6

8

10

12

Predevelopment 1997 (Wet) 2000 (Average) 2001 (Dry) 2000–2007

Outflows Inflows

Water budget flux, in millions of acre-feet per year

Figure 33. Simulated annual water-budget flux for predevelopment, wet (1997), average (2000), and dry (2001) years, and 

current conditions (2000–2007), Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington.
80 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

These 3 years captured the hydrologic variability present 

in the model domain and were representative of the current 

groundwater pumping and recharge from irrigation. The 

net annual budgets for each of the three representative 

years showed that recharge was the largest water-budget 

component and dominated the inflows to the system (fig. 33). 

Groundwater recharge ranged from about 6.2 million acre-ft 

in 2001 to 10.9 million acre-ft in 1997, an increase of about 

76 percent. For the average year (2000), recharge was within 

2 percent of the mean annual recharge for 1950–2007. For 

the wet year (1997), outflows to storage greatly exceeded 

inflows from storage, resulting in an increase of groundwater 

in storage of 1.5 million acre-ft. For the dry year (2001), with 

less recharge and additional effects of pumping (wells), the 

difference is an outflow of groundwater in storage of about 

1.2 million acre-ft.

The simulated water budgets for predevelopment, 

select wet, average, and dry years, and current conditions 

(2000–2007) demonstrate the role of human development 

of the water resources on the groundwater system (fig. 33). 

The predevelopment stress period was steady-state and 

therefore, no change in storage was simulated. Recharge 

in predevelopment conditions was derived solely from 

precipitation and was balanced by outflows from the 

groundwater system to streams. Although recharge increased 

from predevelopment to current conditions due to irrigation, 

streamflow gradually decreased as groundwater extraction 

from wells captured more streamflow and induced inflows 

from storage.

Simulated Effects of Pumping, Commingling 

Wells, and Irrigation Recharge

In order to evaluate the influence of different drawdown 

mechanisms, simulated changes in water levels were divided 

into three components: (1) commingling, (2) pumping, and 

(3) irrigation recharge. The effects of commingling wells, 

pumping wells, and irrigation for the Wanapum unit are shown 

in figures 34A–C. The changes due to commingling (fig. 34A) 

were simulated using commingling wells with zero pumping 

and recharge computed using the Bauer-Vaccaro regression 

equation (no extra recharge from irrigation sources). The 

changes due to pumping (fig. 34B) were estimated as the total 

drawdown (fig. 31) minus the drawdown simulated assuming 

zero pumping (recharge and commingling wells were the 

same for the two runs). The changes due to irrigation recharge 

(fig. 34C) were computed as the drawdown simulated assuming 

zero pumping minus the drawdown-from-commingling wells. 

This method of computation ensures that the sum of the 

component drawdowns equals the total drawdown (fig. 31).

In figures 34A–C, negative values indicate areas of 

buildup or rising water levels and positive numbers indicate 

areas of declining water levels. Assessment of the three 

scenarios shows a few general findings:

(1) The pumping component was usually positive 

(declining water levels) and recharge from the irrigation 

component was usually negative (raising water levels). Minor 

departures from this pattern are associated with the method of 

computing drawdown and water-level rise (thickness weighted 

averaging of cells that are not dry).

(2) The effects of commingling wells show water-level 

declines in some parts of the aquifers and rises in others. 

Commingling wells would mostly cause water-level rises in the 

Wanapum unit, if there was no pumping or change in recharge. 

This pattern is consistent with the downward gradient in the 

upper CRBG on the Palouse Slope.

(3) Generally, pumping dominated water-level changes, 

followed by irrigation enhanced recharge as facilitated by 

commingling pathways. Commingling was a larger factor 

in upland, structurally complex areas where hydraulic-head 

gradients are naturally high. This pattern is similar to the 

pattern determined by Burns and others (2012a) in Mosier, 

Oregon, and could be an important factor on the Palouse Slope.
Numerical Simulation of Groundwater Flow 81

tac13-0856_fig 34a

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

No

symbol -9.9 to 10

-129 to -100

-99.9 to -10

Active model domain

EXPLANATION

10.1 to 100

100.1 to 200

200.1 to 300

300.1 to 400

400.1 to 409.6

Wanapum simulated drawdown, in feet, 

 due to commingling

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

A

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 34. Simulated effects of (A) commingling wells only (no pumping); (B) pumping wells only (no commingling); and 

(C) groundwater recharge from irrigation in the Wanapum basalt unit, Columbia Plateau Regional Aquifer System, Idaho, Oregon, 

and Washington.
82 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

tac13-0856_fig 34b

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

No

symbol -9.9 to 10

Active model domain

EXPLANATION

Wanapum simulated drawdown, in feet, 

 due to pumping

-51.0 to -10

10.1 to 100

100.1 to 200

200.1 to 244.1

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

B

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 34.—Continued
Numerical Simulation of Groundwater Flow 83

tac13-0856_fig 34c

Pasco

Yakima

Moscow

Pullman

Spokane

Lewiston

Wenatchee

The

Dalles

Ellensburg

Walla

Walla

Moses

Lake

Base map modified from U.S. Geological Survey 

and other digital data, various scales. Coordinate 

Reference System: State Plane, Washington 

State, South. Horizontal datum is North 

American Datum of 1983. 

No

symbol -9.9 to 10.0

Active model domain

EXPLANATION

Wanapum simulated drawdown, in feet, 

 due to irrigation

-297.5 to -200

-199.9 to -100

-99.9 to -10

121° 120° 119° 118° 117°

48°

47°

46°

45°

0 20 40 60 MILES

0 20 40 60 KILOMETERS

C

John Day

River Deschutes River

Palouse River

Yakima River

Klickitat River

COLUMBIA

RIVER

Spokane River

Snake River

Lake

Coeur

d'Alene

COLUMBIA RIVER

Umatilla River

Walla Walla River

Study

area

WASHINGTON

OREGON

IDAHO

Figure 34.—Continued
84 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Major Findings From Numerical 

Simulation of Groundwater Flow

The goals of the study of the CPRAS were to 

(1) characterize the hydrologic status of the system, 

(2) identify trends in groundwater storage and use, and 

(3) quantify groundwater availability. The simulation model 

was designed for evaluating and testing the conceptual model 

of the system defined in previous phases of the study and 

for evaluation of groundwater availability. In the process of 

accomplishing these goals, a new or improved understanding 

of certain aspects of the flow system was reached. Major 

findings from the groundwater-flow model include:

• The use of many thin layers improve representation 

of vertical gradients and discrete connectivity 

with surface features. This allowed for a simple 

representation of the permeability distribution and 

fewer barriers to horizontal flow. For example, the 

major hydraulic barrier on the Palouse Slope possibly 

is a linear drainage feature (erosional intersection of 

deeper CRBG lavas at a geologic fold) causing local 

changes in potentiometric surface.

• Hydraulic heads in individual aquifers are 

controlled by drainage features. Frequently, CRBG 

units dip at a steeper angle than the land surface. As 

a result, CRBG aquifers that are not exposed at land 

surface near the center of the basin can be exposed 

in upland river and stream cuts. When the rate of 

groundwater leakage from parts of the aquifer below 

the stream cut contact is lower than the recharge rate, 

then the aquifer fills to the altitude of the contact. In 

this way, a single aquifer may have old, slow moving 

groundwater below the discharge altitude and much 

younger faster-moving water filling and draining the 

altitudes above the spill points. The altitude controlling 

the head in an aquifer commonly is associated with an 

anticlinal fold or other geologic structure that allows 

connection of the aquifer with an erosional cut.

• Horizontal-flow barriers are important for 

representing some hydraulic features. These 

barriers can result from folds, faults, and depositional 

variability. Horizontal-flow barriers can vary with 

depth, being present in either shallow or deep units 

only. Deep barriers are formed by geologic structures 

that predate deposition of younger units. Shallow 

barriers may result from thrust faulting, depositional 

variability, and complex geometrical deformation 

occurring during folding and faulting. In order to 

understand flow barriers, more spatially and temporally 

dense head data may be required.

• Water-level changes in CRBG aquifers are the 

result of pumping, commingling, and changes in 

recharge. In areas with large-scale surface-water 

irrigation projects and commingling wells, model 

simulation results demonstrate that water-level declines 

are partially mitigated as recharge flows into deeper 

aquifers through commingling wells. In areas with 

little surface-water irrigation recharge, commingling 

is expected to result in water-level declines in higher 

head units as water flows from these units into lower 

head units.

• The simulation results identify the effects of the 

simulated processes, but because this analysis was 

not accompanied by an uncertainty analysis, the 

magnitude of the effects is not well constrained.

• Preliminary analyses indicate that commingling may 

vary by region.

• Although not investigated as part of this 

groundwater availability study, there may be 

water-quality implications related to mitigation of 

water-level declines by surface-water recharge (for 

example, transport of nutrients and pesticides into 

deep aquifers).

• Sufficient data do not exist to determine the 

possibility that hydraulic conductivity decreases 

with depth. Hansen and others (1994) suggested 

that hydraulic conductivity decreases at depth due to 

overburden pressure and secondary mineralization. 

Because most of the data are collected above 2,000 ft 

below land surface, and because most of the active 

flow system is above this depth, the simulation 

results were similar, indicating that current data 

are insufficient to distinguish between the two 

competing models. Future research could be directed 

toward distinguishing between the two models, if 

desired predictions are determined to be sensitive to 

the difference.

• Automated parameter estimation tends to 

overestimate vertical hydraulic conductivity of 

the CRBG units. In areas where commingling has 

substantially equilibrated heads between aquifers, few 

water-level measurements represent predevelopment 

conditions, so hydraulic conductivity is overestimated 

(also identified by Burns and others, 2012a). In areas 

where vertical gradients exist, but commingling results 

in highly-variable water-level measurements , the 

use of a least-squares objective function resulted in 

overestimation of vertical hydraulic conductivity as the 

parameter estimation process seeks to fit the centerof-mass of the data, rather than the outliers that likely 
Summary 85

represent the true magnitude of the vertical gradients 

largely unaffected by commingling. Future research 

should include novel approaches to constrain estimates 

of the vertical hydraulic conductivity.

• Groundwater pumping increased substantially since 

the 1970s–1980s; this increase resulted in declining 

water levels at depth and decreased base flows 

over much of the study area. The effects of pumping 

are mitigated somewhat by the increase of surfacewater irrigation, especially in the shallow Overburden 

unit (Konikow, 2013), and commingling wells in 

some areas.

• During dry to average years, groundwater pumping 

causes a net loss of groundwater in storage. Current 

levels of groundwater pumping exceed recharge during 

all but the wettest of years.

Summary

A three-dimensional numerical model of groundwater 

flow was constructed for the Columbia Plateau Regional 

Aquifer System (CPRAS), Idaho, Oregon, and Washington, 

to evaluate and test the conceptual model of the system and 

to evaluate groundwater availability. The model described in 

this report can be used as a tool by water-resource managers 

and other stakeholders to quantitatively evaluate proposed 

alternative management strategies and assess the longterm availability of groundwater. The numerical simulation 

of groundwater flow in the CPRAS was conducted with 

support from the Groundwater Resources Program of the 

U.S. Geological Survey Office of Groundwater.

The model was constructed using the U.S. Geological 

Survey modular three-dimensional finite-difference 

groundwater-flow model, MODFLOW-NWT. The model 

used 3 kilometer (9,842.5 foot) grid cells that subdivided the 

model domain by 126 rows and 131 columns. Vertically, the 

model domain was subdivided into six geologic model units: 

the Overburden, Saddle Mountains Basalt, Mabton Interbed, 

Wanapum Basalt, Vantage Interbed, and Grande Ronde 

Basalt. The geologic units in the model were represented with 

100 model layers. Soil-water balance modeling was used to 

estimate irrigation pumping during 1985–2007 as another 

component of the CPRAS availability study. These estimates 

were supplemented with pumping estimates for municipal, 

industrial, residential, and all other uses. The pumping 

estimates from previous studies were used to estimate 

pumping during 1920–84. Use of the multi-node well package 

for MODFLOW allowed representation of both groundwater 

pumping and cross-connection of aquifers through open well 

boreholes (commingling). Predevelopment groundwater 

recharge was estimated using gridded estimates of annual 

precipitation during 1895–2007. Recharge estimates from 

surface-water and groundwater irrigation application were 

extracted from soil-water balance modeling. The Columbia 

River, major tributaries, and all other surface-water features 

are included in the model as either river cells or drain cells.

Two separate models were constructed to simulate 

groundwater flow in the CPRAS: a steady-state 

predevelopment model representing conditions before 

large-scale pumping and irrigation altered the system, 

and a transient model representing the period 1900–2007. 

Calibration used an iterative approach of automated 

parameter-estimation techniques (steady-state predevelopment 

model) and traditional trial-and-error (transient model) 

methods. The total number of observations used in the 

steady-state and transient model calibrations was 10,525 and 

46,460 water levels, respectively, and 50 base-flow estimates 

for both models.

The model simulates the shape, slope, and trends of 

the potentiometric surface that generally is consistent with 

mapped water levels. For the transient model, the mean 

and median difference between simulated minus measured 

hydraulic heads is -10 and 4 ft, respectively, with a standard 

deviation of 164 ft over a 5,648 ft range of measured heads. 

The residuals for the simulation period show that 52 percent of 

the simulated heads exceeded measured heads with a median 

residual value of 43 ft, and 48 percent were less than measured 

heads with a median residual value of -76 ft.

The CPRAS model was constructed to derive 

components of the groundwater budget and help understand 

the interactions of stresses, such as recharge, groundwater 

pumping, and commingling wells on the groundwater and 

surface-water system. Through these processes, the model can 

be used to identify trends in groundwater storage and use, and 

quantify groundwater availability. The annual groundwater 

budgets show several trends over the simulation period. 

Groundwater pumping was negligible until the 1950s and 

began to increase significantly during the 1970s and 1980s. 

Recharge was highly variable due to the interannual variability 

of precipitation, but began to increase in the late 1940s due to 

the increase in irrigation projects. Streamflow gains and net 

storage followed recharge closely, but the loss of groundwater 

in storage increased and stream base flow decreased as 

groundwater pumping increased.

The simulation model was used to evaluate and test 

the conceptual understanding of the system defined in 

previous phases of the study and to evaluate groundwater 

availability. In the process of accomplishing these goals, a 

new or improved understanding of certain aspects of the flow 

system was reached. Some of the major findings from the 

groundwater-flow model include: 

1. The use of many thin layers improved representation 

of vertical gradients and discrete connectivity with 

surface features; 
86 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

2. Hydraulic heads in individual aquifers are controlled by 

drainage features; 

3. Horizontal flow barriers are important for representing 

some hydraulic features;

4. Water-level changes in Columbia River Basalt aquifers 

are the result of commingling, pumping, and irrigation 

recharge, with groundwater pumping having the 

greatest effect on water levels, followed by irrigation 

enhanced recharge;

5. Groundwater pumping has increased substantially 

since the 1970s–1980s and this increase has resulted in 

declining water levels at depth and decreasing base flows 

over much of the study area; and 

6. During dry to average precipitation years, groundwater 

pumping causes a net loss of groundwater in storage. 

Groundwater pumping presently exceeds recharge during 

all but the wettest of years.

References Cited

Anderman, E.R., and Hill, M.C., 2000, MODFLOW-2000—

The U.S. Geological Survey Modular Ground-Water 

Model-Documentation of the Hydrogeologic-Unit Flow 

(HUF) Package: U.S. Geological Survey Open-File Report 

2000‐342, 89 p., http://pubs.er.usgs.gov/usgspubs/ofr

ofr00342.

Anderson, M.R., and Woessner, W.W., 1992, Applied 

groundwater modeling simulation of flow and advective 

transport: San Diego/New York/Boston/London/Sydney/

Tokyo/Toronto, Academic Press, Inc., 381 p.

Barker, R.A., 1979, Computer simulation and geohydrology 

of a basalt aquifer system in the Pullman-Moscow Basin, 

Washington and Idaho: Washington Department of Ecology 

Water-Supply Bulletin 48, 119 p.

Bauer, H.H., and Hansen, A.J., 2000, Hydrology of the 

Columbia Plateau Regional Aquifer System, Washington, 

Oregon, and Idaho: U.S. Geological Survey WaterResources Investigations Report 96-4106, 61 p., http://pubs.

er.usgs.gov/usgspubs/wri/wri964106.

Bauer, H.H., and Vaccaro, J.J., 1990, Estimates of groundwater recharge to the Columbia Plateau regional 

aquifer system, Washington, Oregon, and Idaho, for 

pre development and current land-use conditions: 

U.S. Geological Survey Water-Resources Investigations 

Report 88-4108, 37 p., http://pubs.er.usgs.gov/usgspubs/wri/

wri884108.

Bergeron, M.P., Freshly, M.D., Last, G.V., and Mitchell, P.J., 

1986, Development of the Hanford three-dimensional 

ground-water flow model: Richland, Wash., Battelle Pacific 

Northwest Laboratory, draft report, 14 p.

Bureau of Reclamation, 1999, Yakima River Basin Water 

Enhancement Project, Washington, Final Programmatic 

Environmental Impact Statement: Yakima, Wash., U.S. 

Department of Interior, Bureau of Reclamation, Pacific 

Northwest Region, Upper Columbia Area Office, 197 p.

Burns, E.R., Morgan, D.S., Lee, K.K., Haynes, J.V., and 

Conlon, T.D., 2012a, Evaluation of long-term waterlevel declines in basalt aquifers near Mosier, Oregon: 

U.S. Geological Survey Scientific Investigations 

Report 2012–5002, 134 p.

Burns, E.R., Morgan, D.S., Peavler, R.S., and Kahle, S.C., 

2011, Three-dimensional model of the geologic framework 

for the Columbia Plateau Regional Aquifer System, Idaho, 

Oregon, and Washington: U.S. Geological Survey Scientific 

Investigations Report 2010-5246, 44 p., http://pubs.usgs.

gov/sir/2010/5246/.

Burns, E.R., Snyder, D.T., Haynes, J.V., and Waibel, M.S., 

2012b, Groundwater status and trends for the Columbia 

Plateau Regional Aquifer System, Washington, Oregon, 

and Idaho: U.S. Geological Survey Scientific Investigations 

Report 2012–5261, 52 p.

Calkins, F.V., 1905, Geology and water resources of a portion 

of east-central Washington: U.S. Geological Survey 

Water‐Supply Paper 118, 96 p

Cline, D.R., 1984, Ground-water levels and pumpage in 

east-central Washington, including the Odessa - Lind area, 

1967 to 1981: Washington State Department of Ecology, 

Water‐Supply Bulletin No. 55, 34 p

Cline, D.R., and Collins, C.A., 1992, Ground-water pumpage 

from the Columbia Plateau, Washington and Oregon, 

1945 to 1984: U.S. Geological Survey Water-Resources 

Investigations Report 90-4085, 31 p., 5 pls., http://pubs.

er.usgs.gov/publication/wri904085.

Davies-Smith, A., Bolke, E.L., and Collins, C.A., 1988, 

Geohydrology and digital simulation of the ground-water 

flow system in the Umatilla Plateau and Horse Heaven 

Hills area, Oregon and Washington: U.S. Geological Survey 

Water-Resources Investigations Report 87-4268, 72 p., 

http://pubs.er.usgs.gov/usgspubs/wri/wri874268.

Doherty, J.E., 2010, PEST, Model-independent parameter 

estimation—User manual (5th ed., with additions): 

Brisbane, Australia, Watermark Numerical Computing, 

variously paged.
References Cited 87

Doherty, J.E., and Hunt, R.J., 2010, Approaches to highly 

parameterized inversion—A guide to using PEST for 

groundwater-model calibration: U.S. Geological Survey 

Scientific Investigations Report 2010–5169, 59 p., http://

pubs.usgs.gov/sir/2010/5169/.

Drost, B.W., Ely, D.M., and Lum, W.E., 1999, Conceptual 

model and numerical simulation of the ground-water-flow 

system in the unconsolidated sediments of Thurston County, 

Washington: U.S. Geological Survey Water-Resources 

Investigations Report 99-4165, 106 p., http://pubs.usgs.gov/

wri/wri994165/.

Drost, B.W., Whiteman, K.J., and Gonthier, J.B., 1990, 

Geologic framework of the Columbia Plateau aquifer 

system, Washington, Oregon, and Idaho: U.S. Geological 

Survey Water-Resources Investigations Report 87-4238, 

10 p., http://pubs.er.usgs.gov/publication/wri874238.

Ely, D.M., Bachmann, M.P., and Vaccaro, J.J., 2011, 

Numerical simulation of groundwater flow for the Yakima 

River Basin aquifer system, Washington: U.S. Geological 

Survey Scientific Investigations Report 2011-5155, 90 p., 

http://pubs.usgs.gov/sir/2011/5155/.

Faunt, C.C., Blainey, J.B., Hill, M.C., D’Agnese, F.A., and 

O’Brien, G.M., 2010, Transient numerical model, chapter F 

of Belcher, W.R. ed., Death Valley regional groundwater 

flow system, Nevada and California—Hydrogeologic 

framework and transient groundwater flow model: 

U.S. Geological Professional Paper 1711, p. 257–344. 

(Also available at http://pubs.usgs.gov/pp/1711/.)

Faunt, C.C., Hanson, R.T., Belitz, Kenneth, Schmid, 

Wolfgang, Predmore, S.P., Rewis, D.L., and McPherson, 

Kelly, 2009, Numerical model of the hydrologic landscape 

and groundwater flow in California’s Central Valley, 

chapter C of Faunt C.C., ed., Groundwater availability of 

the Central Valley aquifer, California: U.S. Geological 

Survey Professional Paper 1766, p. 121–212, http://pubs.

usgs.gov/pp/1766/.

Gannett, M.W., and Lite, K.E., Jr., 2004, Simulation of 

regional ground-water flow in the upper Deschutes 

Basin, Oregon: U.S. Geological Survey Water-Resources 

Investigations Report 03–4195, 84 p.

Hansen, A.J., Vaccaro, J.J., and Bauer, H.H., 1994, Groundwater flow simulation of the Columbia Plateau regional 

aquifer system, Washington, Oregon, and Idaho: 

U.S. Geological Survey Water-Resources Investigations 

Report 91-4187, 101 p., 15 pls., http://pubs.er.usgs.gov/

usgspubs/wri/wri914187.

Harbaugh, A.W., 2005, MODFLOW-2005—The U.S. 

Geological Survey modular groundwater model—The 

ground-water flow process: U.S. Geological Survey 

Techniques and Methods 6-A16, variously paged, http://

pubs.usgs.gov/tm/2005/tm6A16/. 

Hill, M.C., 1998, Methods and guidelines for effective model 

calibration: U.S. Geological Survey Water-Resources 

Investigations Report 98-4005, 90 p.

Hsieh, P.A., and Freckleton, J.R., 1993, Documentation of a 

computer program to simulate horizontal-flow barriers using 

the U.S. Geological Survey’s modular three-dimensional 

finite-difference ground-water flow model: U.S. Geological 

Survey Open-File Report 92-477, 32 p., http://pubs.er.usgs.

gov/publication/ofr92477.

Jones, M.A., and Vaccaro, J.J., 2008, Extent and depth to top 

of basalt and interbed hydrogeologic units, Yakima River 

basin aquifer system, Washington: U.S. Geological Survey 

Scientific Investigations Report 2008-5045, 22 p., http://

pubs.usgs.gov/sir/2008/5045/.

Jones, M.A., Vaccaro, J.J., and Watkins, A.M., 2006, 

Hydrogeologic framework of sedimentary deposits in 

six structural basins, Yakima River basin, Washington: 

U.S. Geological Survey Scientific Investigations Report 

2006-5116, 24 p., http://pubs.usgs.gov/sir/2006/5116/.

Kahle, S.C., Morgan, D.S., Welch, W.B., Ely, D.M., Hinkle, 

S.R., Vaccaro, J.J., and Orzol, L.L., 2011, Hydrogeologic 

framework and hydrologic budget components of the 

Columbia Plateau Regional Aquifer System, Washington, 

Oregon, and Idaho: U.S. Geological Survey Scientific 

Investigations Report 2011–5124, 66 p., http://pubs.usgs.

gov/sir/2011/5124/.

Kahle, S.C., Olsen, T.D., and Morgan, D.S., 2009, Geologic 

setting and hydrogeologic units of the Columbia Plateau 

Regional Aquifer System, Washington, Oregon, and Idaho: 

U.S. Geological Survey Scientific Investigations Map 3088, 

1 sheet, http://pubs.usgs.gov/sim/3088/.

Kinnison, H.B., and Sceva, J.E., 1963, Effects of hydraulic 

and geologic factors on streamflow of the Yakima River 

Basin, Washington: U.S. Geological Survey Water-Supply 

Paper 1595, 134 p., http://pubs.er.usgs.gov/usgspubs/wsp/

wsp1595.

Konikow, L.F., 2013, Groundwater depletion in the United 

States (1900–2008): U.S. Geological Survey Scientific 

Investigations Report 2013−5079, 63 p., http://pubs.usgs.

gov/sir/2013/5079.

Konikow, L.F., Hornberger, G.Z., Halford, K.J., and Hanson, 

R.T., 2009, Revised multi-node well (MNW2) package for 

MODFLOW ground-water flow model: U.S. Geological 

Survey Techniques and Methods 6–A30, 67 p.

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

Lava Plateau, in Back, William, Rosenshein, J.S., and 

Seaber, P.R., eds., Hydrogeology, v. O-2—The geology 

of North America: Boulder, Colo., Geological Society of 

America, p. 37–50.
88 Numerical Simulation of Groundwater Flow in the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington

Lum, W.E., II, Smoot, J.L., and Ralston, D.R., 1990, 

Geohydrology and numerical model analysis of groundwater flow in the Pullman-Moscow area, Washington 

and Idaho: U.S. Geological Survey Water-Resources 

Investigations Report 89-4103, 73 p., http://pubs.er.usgs.

gov/usgspubs/wri/wri894103.

Magirl, C.S., and Olsen, T.D., 2009, Navigability potential 

of Washington rivers and streams determined with 

hydraulic geometry and a geographic information system: 

U.S. Geological Survey Scientific Investigations Report 

2009–5122, 22 p., http://pubs.usgs.gov/sir/2009/5122/.

Newcomb, R.C., 1965, Geology and ground-water resources 

of the Walla Walla River Basin, Washington-Oregon: 

Washington Division of Water Resources, Water-Supply 

Bulletin 21, 151 p., accessed January 31, 2011, at http://

pubs.er.usgs.gov/publication/70047774.

Newcomb, R.C., 1969, Effect of tectonic structure on the 

occurrence of ground water in the basalt of the Columbia 

River Group of The Dalles area, Oregon and Washington: 

U.S. Geological Survey Professional Paper 383-C, 33 p., 

1 pl., http://pubs.er.usgs.gov/publication/pp383C.

Niswonger, R.G., Panday, Sorab, and Ibaraki, Motomu, 

2011, MODFLOW-NWT, A Newton formulation for 

MODFLOW-2005: U.S. Geological Survey Techniques and 

Methods 6–A37, 44 p.

Packard, F.A., Hansen, A.J., Jr., and Bauer, H.H., 1996, 

Hydrogeology and simulation of flow and the effects of 

development alternatives on the basalt aquifers of the Horse 

Heaven Hills, south-central Washington: U.S. Geological 

Survey Water-Resources Investigations Report 94-4068, 

92 p., 2 pls., http://pubs.er.usgs.gov/publication/wri944068.

Parker, G.L., and Storey, F.B., 1916, Water powers of the 

Cascade Range, Part III—Yakima River Basin: U.S. 

Geological Survey Water Supply Paper 369, 169 p., 16 pls., 

http://pubs.er.usgs.gov/publication/wsp369.

Piper, A.M., 1932, Geology and ground-water resources 

of the Dalles region, Oregon: U.S. Geological Survey 

Water‐Supply Paper 659-B, p. 107–189

Porcello, J., Kindsey, K., Tonkin, M., Muffels, C., and 

Karanovic, M., 2010, Groundwater model development 

process for the Columbia Basin Ground Water Management 

Area of Adams, Franklin, Grant, and Lincoln Counties, 

Washington: Prepared by the Columbia Basin Ground Water 

Management Area of Adams, Franklin, Grant and Lincoln 

Counties, variously paged.

PRISM Climate Group, 2004, PRISM Spatial Climate 

Datasets: Corvallis, Oregon State University, accessed 

October 1, 2009, at http://www.prismclimate.org.

Reidel, S.P., Johnson, V.G., and Spane, F.A., 2002, Natural 

gas storage in basalt aquifers of the Columbia Basin, 

Pacific Northwest USA—A guide to site characterization: 

Richland, Wash., Pacific Northwest National Laboratory, 

variously paged, accessed November 3, 2010, at http://

www.pnl.gov/main/publications/external/technical_reports/

PNNL-13962.pdf.

Reilly, T.E., Dennehy, K.F., Alley, W.M., and Cunningham, 

W.L., 2008, Ground-water availability in the United States: 

U.S. Geological Survey Circular 1323, 70 p., http://pubs.

usgs.gov/circ/1323/.

Schwennsen, A.T., and Meinzer, O.E., 1918, Ground water 

in Quincy Valley, Washington: U.S. Geological Survey 

Water‐Supply Paper 425-E, p. 131–161

Senay, G.B., Budde, M., Verdin, J.P., and Melesse, 

A.M., 2007, A coupled remote sensing and simplified 

surface energy balance approach to estimate actual 

evapotranspiration from irrigated fields: Sensors, v. 7, 

p. 979–1000.

Simley, J.D., and Carswell, Jr., W.J., 2009, The National 

Map—Hydrography: U.S. Geological Survey Fact 

Sheet 2009-3054, 4 p.

Simons, W.D., 1953, Irrigation and streamflow depletion in 

Columbia River basin above The Dalles, Oregon: U.S. 

Geological Survey Water Supply Paper 1220, 126 p., 1 pl., 

http://pubs.er.usgs.gov/usgspubs/wsp/wsp1220.

Sinclair, K.A., and Pitz, C.F., 1999, Estimated baseflow 

characteristics of selected Washington rivers and streams: 

Washington State Department of Ecology, Water-Supply 

Bulletin No. 60, 181 p.

Sloto, R.A., and Crouse, M.Y., 1996, HYSEP—A computer 

program for streamflow hydrograph separation and analysis: 

U.S. Geological Survey Water Resources Investigation 

Report 96-4040, 46 p.

Smith, G.O., 1901, Geology and water resources of a portion 

of Yakima County, Washington: U.S. Geological Survey 

Water-Supply Paper 55, 68 p.

Snyder, D.T., and Haynes, J.V., 2010, Groundwater conditions 

during 2009 and changes in groundwater levels from 1984 

to 2009, Columbia Plateau Regional Aquifer System, 

Washington, Oregon, and Idaho: U.S. Geological Survey 

Scientific Investigations Report 2010–5040, 12 p., http://

pubs.usgs.gov/sir/2010/5040/.

State of Washington Office of Financial Management, 2014, 

Washington’s Rank in the Nations Agriculture—2013 

Data Book: State of Washington Office of Financial 

Management, accessed March 21, 2014, at http://www.ofm.

wa.gov/databook/pdf/databook.pdf.
References Cited 89

Stearns, R.T., 1942, Hydrology of lava-rock terranes, 

chap. XV, in Meinzer, O.E., ed., Hydrology: New York, 

McGraw‐Hill, 712 p

Tanaka, H.H., Hansen, A.J., Jr., and Skrivan, J.A., 1974, 

Digital-model study of ground-water hydrology, Columbia 

Basin Irrigation Project area, Washington: State of 

Washington Department of Ecology Water-Supply 

Bulletin 40, 60 p.

Taylor, G.C., Jr., 1948, Ground water in the Quincy Basin, 

Wahluke Slope, and Pasco Slope subareas of the Columbia 

Basin Project, Washington: U.S. Geological Survey 

Open‐File Report, 182 p

Tolan, T.L., Reidel, S.P., Beeson, M.H., Anderson, J.L., Fecht, 

K.R., and Swanson, D.A., 1989, Revisions to the estimates 

of the areal extent and volume of the Columbia River Basalt 

Group, in Reidel, S.P., and Hooper, P.R., eds., Volcanism 

and tectonism in the Columbia River flood-basalt province: 

Boulder, Colo., Geological Society of America Special 

Paper 239, p. 1-20, doi: 10.1130/SPE239-p1.

Troutman, B.M., 1985, Errors and parameter estimation 

in precipitation-runoff modeling 2—Case Study: Water 

Resources Research, v. 21, no. 8, p. 1,214–1,222.

U.S. Department of Agriculture, 2007, Census of 

agriculture—2007 census publications: U.S. Department 

of Agriculture, accessed June 16, 2009, at http://www.

agcensus.usda.gov/Publications/2007/Online_Highlights/

Rankings_of_Market_Value/Washington/index.asp.

U.S. Geological Survey, 2004, National hydrography 

dataset—Watershed boundary dataset: U.S. Geological 

Survey database, accessed March 13, 2014, at http://nhd.

usgs.gov/index.html.

Vaccaro, J.J., 1986, Plan of study for the regional aquifersystem analysis, Columbia Plateau, Washington, northern 

Oregon, and northwestern Idaho: U.S. Geological Survey 

Water-Resources Investigations Report 85-4151, 25 p.

Vaccaro, J.J., 1999, Summary of the Columbia Plateau 

Regional Aquifer-System Analysis, Washington, 

Oregon, and Idaho: U.S. Geological Survey Professional 

Paper 1413-A, 51 p.

Vaccaro, J.J., Jones, M.A., Ely, D.M., Keys, M.E., Olsen, 

T.D., Welch, W.B., and Cox, S.E., 2009, Hydrogeologic 

framework of the Yakima River Basin aquifer system, 

Washington: U.S. Geological Survey Scientific 

Investigations Report 2009–5152, 106 p., http://pubs.usgs.

gov/sir/2009/5152/.

Vermeul, V.R., Bergeron, M.P., Cole, C.R., Murray, C.J., 

Nichols, W.E., Scheibe, T.D., Thorne, P.D., Waichler, S.R., 

and Xie, Y., 2003, Transient inverse calibration of the sitewide groundwater-flow model (ACM-2)—FY03 progress 

report: Richland, Wash., Pacific Northwest National 

Laboratory, PNNL-14398, variously paged.

Vermeul, V.R., Cole, C.R., Bergeron, M.P., Thorne, P.D., 

and Wurstner, S.K., 2001, Transient inverse calibration 

of site-wide groundwater model to Hanford Operational 

Impacts from 1943 to 1996—Alternative conceptual model 

considering interaction with uppermost basalt confined 

aquifer: Richland, Wash., Pacific Northwest National 

Laboratory, PNNL-13623, variously paged.

Waring, G.A., 1913, Geology and water resources of a portion 

of south-central Washington: U.S. Geological Survey WaterSupply Paper 316, 46 p.

Weiss, Emanuel, 1982, A computer program for calculating 

relative-transmissivity input arrays to aid model calibration: 

U.S. Geological Survey Open-File Report 82-447, 18 p.

Whiteman, K.J., Vaccaro, J.J., Gonthier, J.B., and Bauer, H.H., 

1994, The hydrogeologic framework and geochemistry 

of the Columbia Plateau Regional Aquifer System in 

Washington, Oregon, and Idaho: U.S. Geological Survey 

Professional Paper 1413-B, 73 p.

Wurstner, S.K., Thorne, P.D., Chamness, M.A., Freshley, 

M.D., and Williams, M.D., 1995, Development of a 

three-dimensional groundwater model of the Hanford 

Site unconfined aquifer system—FY 1995 status 

report: Richland, Wash., Pacific Northwest Laboratory, 

PNL‐10886, variously paged

Publishing support provided by the U.S. Geological Survey

Publishing Network, Tacoma Publishing Service Center

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

 Director, Washington Water Science Center

U.S. Geological Survey

934 Broadway, Suite 300

Tacoma, Washington 98402 

http://wa.water.usgs.gov
— Ely and others Numerical Simulation of Groundwater Flow, Columbia Plateau Regional Aquifer System, Idaho, Oregon, and Washington—Scientific Investigations Report 2014–5127, ver. 1.1

ISSN 2328-0328 (online)

http://dx.doi.org/10.3133/sir20145127