U.S. Department of the Interior

U.S. Geological Survey

Scientific Investigations Report 2015–5149

Prepared in cooperation with the U.S. Environmental Protection Agency

Particle Tracking for Selected Groundwater Wells in the 

Lower Yakima River Basin, Washington
Cover: Photograph of Sunnyside Canal looking north, northwest approximately 5 miles southeast of the Sunnyside Diversion Dam on the 

Yakima River, Yakima, Washington, 2002. Photograph by David Walsh, Bureau of Reclamation.
Particle Tracking for Selected Groundwater 

Wells in the Lower Yakima River Basin, 

Washington

By Matthew P. Bachmann

Prepared in cooperation with the U.S. Environmental Protection Agency

Scientific Investigations Report 2015–5149

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: 2015

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

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:

Bachmann, M.P., 2015, Particle tracking for selected groundwater wells in the lower Yakima River Basin, Washington: 

U.S. Geological Survey Scientific Investigations Report 2015-5149, 33 p., http://dx.doi.org/10.3133/sir20155149.

ISSN 2328-0328 (online)
iii

Contents

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

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

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

Description of Study Area ...................................................................................................................3

Simulation of Particle Tracking....................................................................................................................6

Stress Periods .......................................................................................................................................6

Discretization.........................................................................................................................................6

Model Inputs........................................................................................................................................10

Calibration............................................................................................................................................14

Particle Tracking .................................................................................................................................14

Delineation of Contributing Areas....................................................................................................21

Model Uncertainties and Limitations ..............................................................................................26

Summary........................................................................................................................................................28

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

Appendix A. Summarized Particle-Tracking Information.....................................................................31
iv

Figures

1. Map showing the Yakima River Basin, Washington ...............................................................2

2. Map showing location of structural basins and geologic structures, Yakima River

Basin, Washington........................................................................................................................4

3. Map showing simplified surficial geology in the Yakima River Basin, Washington..........5

4. Map showing location and extent of the groundwater-flow model grid used to

simulate the Yakima River Basin aquifer system, Washington.............................................7

5. Map showing location of model cross section through the Toppenish/Benton

subbasin and hydrogeologic section showing relation of model hydrogeologic units

to model computational layers, Yakima River Basin aquifer system, Washington............8

6. Map showing location of the groundwater-flow model boundary conditions for

drains, general and specified heads, and streamflow routing cells, Yakima River

Basin aquifer system, Washington..........................................................................................11

7. Map showing spatial distribution of mean annual recharge simulated for

1960–2001, Yakima River Basin aquifer system, Washington..............................................13

8. Map showing location of model structures and flow barrier conductivity, Yakima

River Basin aquifer system, Washington................................................................................15

9. Graph showing weighted measured values as a function of weighted simulated

values, Yakima River Basin aquifer system, Washington ....................................................17

10. Map showing locations of 121 groundwater samples with nitrate 

concentrations exceeding 10 milligrams per liter, lower Yakima River Valley, 

Washington, 2008–2011..............................................................................................................18

11. Map showing simulated particle histories for particles that ended the 42-year

simulation period at the location of model cells containing selected groundwater

wells in the lower Yakima River Basin, Washington ............................................................19

12. Maps showing modeled structural barriers to groundwater flow in some locations

cause particle flow paths to change direction sharplyand seasonally varying

recharge patterns and associated changes in hydraulic head gradients cause

some particle paths to temporarily reverse direction, Yakima River Basin,

Washington..................................................................................................................................20

13. Maps showing recharge capture areas defined by particle ending locations

outlining the area contributing recharge, and all area circumscribed by particle

pathways outlining the areal extent of the zone of contribution for a group of six

groundwater sampling locations near Zillah, Washington..................................................22

14. Maps showing zones of contribution for a dense area of wells near Granger,

Washington, divided into two groups to make individual particle pathlines more

visible in a location where areas contributing recharge to specific wells overlap

heavily...........................................................................................................................................24

15. Map showing simulated zones of contribution for the model cells containing the

121 contaminated wells in the lower Yakima River Basin, Washington............................27

Tables

1. Descriptive statistics for the thickness of groundwater model layers used to

simulate the Yakima River Basin aquifer system, Washington...........................................10

2. Calibrated model parameter values for the Yakima River Basin aquifer system,

Washington..................................................................................................................................16
v

Conversion Factors

Inch/Pound to International System of Units

Multiply By To obtain

Length

foot (ft) 0.3048 meter (m)

mile (mi) 1.609 kilometer (km)

Area

acre 4,047 square meter (m2)

acre 0.004047 square kilometer (km2)

square foot (ft2) 0.09290 square meter (m2)

square mile (mi2) 2.590 square kilometer (km2) 

Volume

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

acre-foot (acre-ft) 0.001233 cubic hectometer (hm3) 

Flow rate

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

acre-foot per year (acre-ft/yr) 0.001233 cubic hectometer per year (hm3/yr)

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

square foot per day (ft2/d) 0.09290 square meter per day (m2d)

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

International System of Units to Inch/Pound

Multiply By To obtain

Volume

liter (L) 1.057 quart (qt)

Mass

milligram (mg) 0.00003527396 ounce, avoirdupois (oz)

Concentrations of chemical constituents in water are given in milligrams per liter (mg/L).

Datums

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

(NAVD 88).

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

and the North American Datum of 1983 (NAD 83).

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

Particle Tracking for Selected Groundwater Wells in the 

Lower Yakima River Basin, Washington

By Matthew P. Bachmann

Abstract

The Yakima River Basin in south-central Washington 

has a long history of irrigated agriculture and a more recent 

history of large-scale livestock operations, both of which 

may contribute nutrients to the groundwater system. Nitrate 

concentrations in water samples from shallow groundwater 

wells in the lower Yakima River Basin exceeded the 

U.S. Environmental Protection Agency drinking-water 

standard, generating concerns that current applications of 

fertilizer and animal waste may be exceeding the rate at 

which plants can uptake nutrients, and thus contributing to 

groundwater contamination.

The U.S. Geological Survey (USGS) recently completed 

a regional scale transient three-dimensional groundwater-flow 

model of the Yakima River Basin using MODFLOW-2000. 

The model was used with the USGS particle-tracking code 

MODPATH to generate advective flowpaths and associated 

travel times. Analyses used particle backtracking in time 

from September 2001 through 504 monthly stress periods 

to October 1959 or until pathlines terminated at a model 

boundary. The particle starting locations were assigned to 

1,000 foot square computational model cells containing one 

or more of the 121 sampling locations with measured nitrate 

concentrations greater than the U.S. Environmental Protection 

Agency drinking-water standard for nitrate (10 milligrams per 

liter [mg/L]). Of the 2,403 particles, the simulated pathlines 

for 2,080 reached the water table within the 42-year simulation 

period, thus identifying the predicted recharge areas for those 

particles. The median horizontal straight-line distance was 

13,194 feet between starting and ending locations for these 

particles and the median time-of-travel for particles that 

intersected the water table was 984 days. Well to water-table 

travel times for 75.4 percent of the particles were less than 

the average travel time of 3,749 days. Predicted recharge 

locations for all particles, including those that did not reach 

the water table in 42 years, were between 50 feet and 34 miles 

horizontal distance from their starting locations, with a median 

distance of less than 3 miles away. 

Generalized groundwater-flow directions in 

unconsolidated basin-fill deposits were towards the Yakima 

River, which acts as a local sink for shallow groundwater, 

and roughly parallel to topographic gradients. Particles 

backtracked from more shallow aquifer locations traveled 

shorter distances before reaching the water table than particles 

from deeper locations. Flowpaths for particles starting at 

wells completed in the basalt units underlying the basin-fill 

deposits sometimes were different than for wells with similar 

lateral locations but more shallow depths. In cases where 

backtracking particles reached geologic structures simulated 

as flow barriers, abrupt changes in direction in some particle 

pathlines suggest significant changes in simulated hydraulic 

gradients that may not accurately reflect actual conditions. 

Most groundwater wells sampled had associated zones of 

contribution within the Toppenish/Benton subbasin between 

the well and the nearest subbasin margin, but interpretation of 

these results for any specific well is likely to be complicated 

by the assumptions and simplifications inherent in the model 

construction process. Delineated zones of contribution for 

individual wells are sensitive to the depths assigned to the 

screened interval of the well, resulting in simulated areal 

extents of the zones of contribution to a discharging well that 

are elongated in the direction of groundwater flow. 

Introduction

The Yakima River Basin in south-central Washington 

covers 6,200 mi2 and includes parts of Yakima, Kittitas, and 

Benton Counties, and part of the Yakama Indian Reservation 

(fig. 1). Several water-quality investigations in recent 

decades have measured nitrate concentrations exceeding 

the U.S. Environmental Protection Agency (EPA) drinking 

water standard of 10 mg/L in shallow drinking-water wells, 

primarily in the lower valley Toppenish/Benton subbasin 

(U.S. Environmental Protection Agency, 2012, Washington 

State Department of Ecology 2010). The contaminated 

drinking-water wells tend to be shallow, and in many cases 

are primarily used by Tribal members of the Yakama Nation 

and Spanish-speaking families participating in the local 

agricultural economy (Washington State Department of 

Ecology, 2010).
2 Particle Tracking for Selected Groundwater Wells in the Lower Yakima River Basin, Washington

tac15-1019_fig01

82

22

22

97

97

97

12

12

82

82

240

24

410

240

90

90

90

24

410

12

97

12

241

82

82 12

KITTITAS COUNTY

YAKIMA COUNTY

COUNTY

BENTON COUNTY

KLICKITAT

COUNTY

YAKIMA

EXPLANATION

Boundary of Yakima River basin

Canal

Roza

Sunnyside Canal

Main Canal

Extension

Kennewick

Main Canal

Marion Drain

North

Branch

Canal

Kittitas Main Canal

Bumping

Lake N Fk

S Fk

S Fk

N Fk

S Fk

N Fk

N Fk

Cle Elum River

Swauk Creek

Caribou Cr

Coleman Cr

Naneum Cr

Cooke Cr

Crow Creek

Bumping River

American River

Rattlesnake Cr

Little Rattlesnake Cr

Nile Creek

Manastash

Badger Creek

Lmuma Cr

Ahtanum Creek

Oak Creek

Tieton River

Wide Hollow Cr

Wenas Creek

Selah Creek

North Fork

Dry Creek

Simcoe Creek

Agency Creek

Creek

Toppenish

Corral Cr

RIVER

Spring Cr

Logy

Creek

Satus

Creek

Mule Dry

Creek

Cle

Elum

Lake

Rimrock

Lake

Kachess Lake

Keechelus Lake

Lake

Wallula

Priest

Rapids

Lake

Teanaway River Reecer Cr

Wilson Creek

Park Cr

Taneum Creek

Umtanum Creek

Little Naches River

Cr

Cowiche Creek

Cowiche Creek

River

Naches

Snipes Cr

Sulphur Cr Wasteway

Cherry

Cr

S Fk

N Fk

YAKIMA

YAKIMA

RIVER

Cle Elum

Parker

Wapato

Ellensburg

Umtanum

Kittitas

Thrall

Toppenish

Moxee City Union Gap

Yakima

Selah

Naches

Sunnyside Richland

Kiona

Kennewick Prosser

Grandview

Mabton

Satus

Zillah

Granger

Roza

Easton

Teanaway

Tieton Rimrock

Snoqualmie

Pass

Van Epps

Pass

Esmeralda

Peaks

KITTITAS

VALLEY

M A N A S T A S H

KACHESS RIDGE

R I D G E

UNION GAP

CLEMAN

MOUNTAIN

TEANAWAY

NORTH RIDGE

U M T A N U M

R I D G E

Deception

Pass

YAKIMA

 VALLEY

AHTANUM RIDGE

COWICHE MTN

AMERICAN RIDGE

CLE ELUM

RIDGE

RIDGE

TOPPENISH RIDGE

R A T T L E S N A K E H I L L S

Y A K I M A R I D G E

H O R S E

H E A V E N

H I L L S

C

A

S C A

D E R

A

N

G

E

W E N AT C H E E

M O U N TA I N S

Gilbert

Peak

Yakama Nation 

Reservation

121°00'

47°

00'

47°

30'

46°

00'

121°30' 30' 120°00' 119°30'

30'

0 5 10 20 30 40 MILES

0 5 10 20 30 40 KILOMETERS

COLUMBIA

RIVER

Base from U.S. Geological Survey and/or

Washington Division of Geology and Earth Resources

digital data, 1983, 1:100,000

Universal Transverse Mercator projection, Zone 10

Horizontal Datum: North American Datum of 1927 (NAD 27)

Yakima

River

Basin

WASHINGTON

Section of the Yakima basin 

where particle tracking was done

Figure 1. Yakima River Basin, Washington.
Introduction 3

Substantial amounts of nitrogen have been applied to 

the land in the Yakima River Basin for decades, primarily 

as fertilizer for row crops, but in recent years also as 

animal waste generated by more than 60 animal feeding 

operations (Washington State Department of Ecology, 2010). 

Although water-quality sampling for nutrients in the basin 

has a long history (Van Winkle, 1914), there has been little 

work attempting to link specific contaminated wells to any 

identified land uses or locations of contaminant loading. The 

completion of a transient three-dimensional groundwaterflow model of the Yakima River Basin (Ely and others, 2011), 

which documents groundwater-flow directions and hydraulic 

conductivities for each mapped hydrogeologic unit in the 

basin, has made possible the identification of simulated travel 

times and contributing areas for specific groundwater wells. 

For this purpose, the particle-tracking code MODPATH 5.0 

(Pollock, 1994) was used for particle backtracking from 

the estimated point of extraction for selected groundwater 

sampling locations to simulated source areas. For pathlines 

that terminated at the water table, the particle tracking analysis 

connects water from a groundwater well to an area of land 

surface that contributed recharge to that well, and is not 

specific to any groundwater contaminant or land use practice.

Purpose and Scope

The purpose of this report is to describe the results of 

using an existing transient groundwater-flow model of the 

Yakima River Basin in conjunction with a particle-tracking 

program to document localized groundwater-flow directions 

and velocities, to estimate time-of-travel for groundwater to 

flow from the water table to selected wells, and to identify 

locations where nitrates applied to land surface may contribute 

to nitrate concentrations measured in groundwater samples 

that exceed the EPA drinking-water standard of 10 mg/L 

(U.S. Environmental Protection Agency, 2012).

The MODFLOW-2000 groundwater-flow model of 

the Yakima Basin developed by Ely and others (2011) was 

used to generate flow directions for every computational 

model cell. That model was constructed for the purpose of 

evaluating regional scale water-management scenarios and 

exchanges between groundwater and surface-water features 

such as rivers, streams, and tile drains. The size of the model 

cells and duration of the model time steps were formulated to 

address regional scale system responses and features and are 

not ideally suited for the analysis of individual wells. Thus, 

the results of this particle-tracking analysis only approximates 

actual groundwater-flow paths.

Although the groundwater sampling locations (wells) 

used to determine particle-tracking starting locations in this 

analysis were selected because of their measured nitrate 

concentrations, no part of the nitrogen cycle is explicitly 

simulated by either MODFLOW-2000 or MODPATH. 

Particle-tracking starting, ending, and intermediate locations 

are provided based on groundwater flowpath travel time 

but concentrations of nitrate or any other solute are not 

simulated or reported. The particle-tracking results can be 

used to identify the direction and rate of groundwater flow at 

various depths in the aquifer system, the land-surface areas 

that contribute recharge to a specific simulated model cell 

containing one or more groundwater wells, and the time-oftravel required for groundwater to move along the simulated 

flowpaths. This analysis considers only nonreactive advective 

transport and as such macro-scale hydrodynamic dispersion, 

diffusion and other processes that may attenuate or retard 

contaminant concentrations, such as adsorption or chemical 

transformations, are not considered. Thus, the results are 

most applicable to chemically conservative groundwater 

contaminants that are fully miscible in groundwater. EPA 

measured nearly 200 different analytes in these samples, 

including major ions and trace metals, perchlorate, fecal 

coliform, pesticides, and wastewater pharmaceuticals, and 

identified elevated concentrations in multiple drinking-water 

sources (U.S. Environmental Protection Agency, 2012). 

Description of Study Area

The modeled area of the Yakima River Basin in 

Washington extends from the crest of the Cascade Mountains 

east to the Columbia River. Unregulated streamflow averages 

4.1 million acre-ft/yr, and is subject to five managed storage 

reservoirs and diversions to irrigated agriculture of about 

2.8 million acre-ft/yr (Ely and others, 2011). The basin is 

separated into several broad valleys by east-west trending 

anticlinal ridges that separate valley floors that slope gently 

toward the Yakima River. These valleys contain most of 

the population and economic activity in the basin. The 

Toppenish/Benton subbasin, the largest and most southerly of 

the five structural subbasins, encompasses almost 1,500 mi2

between Ahtanum Ridge and Rattlesnake Hills on the north 

and Toppenish Ridge and the Horse Heaven Hills on the 

south (fig. 2).

The basin-fill deposits (fig. 3) consist predominantly 

of alluvial fan, loess, dune sand, alluvial terrace, continental 

sedimentary, and Touchet Formation and Ellensburg 

Formation deposits. These sediments reach a maximum 

documented thickness of more than 1,200 ft west of the 

town of Toppenish and gradually decrease in thickness to 

0 ft at the northern end of the basin where the underlying 

basalt outcrops, but remain about 1,000 ft deep south of the 

Yakima River until the basin abruptly ends against the steeply 

emergent Toppenish Ridge (Jones and others, 2006).

About 34 percent of the basin’s population is not served 

by public supply systems and must rely on private wells 

(Washington State Department of Ecology, 2010). Residential 

drinking-water wells in the area typically are shallow, with 

some logged wells as shallow as 10 ft deep. Wells drilled for 

commercial operations, which typically require a higher flow 

rate, are more likely to be drilled through the entire thickness 

of basin sediments to penetrate the Wanapum or Saddle 

Mountain basalts that underlie the basin.
4 Particle Tracking for Selected Groundwater Wells in the Lower Yakima River Basin, Washington

tac15-1019_fig02

121°00'

47°

00'

47°

30'

46°

00’

121°30' 30' 120°00' 119°30'

30'

KITTITAS COUNTY

YAKIMA COUNTY

COUNTY

BENTON COUNTY

KLICKITAT

COUNTY

YAKIMA

Bumping

Lake

S Fk

N Fk

Cle Elum River

Bumping River

Ahtanum Creek

Tieton River Wenas Creek

North Fork

Dry Creek

To

Creek ppenish

RIVER

YAKIMA

Satus

Creek

Cle

Elum

Lake

Rimrock

Lake

Kachess Lake

Keechelus Lake Teanaway River

Umtanum Creek

Little Naches River

Cowichee Creek

Cowichee Creek

River

Naches

RIVER

YAKIMA

RIVER

YAKIMA

COLUMBIA

RIVER

COLUMBIA RIVER

So. Fork

No. Fork

Cle Elum

Parker

Wapato

Ellensburg

Umtanum

Kittitas Thrall

Toppenish

Moxee

City Union Gap

Yakima

Selah

Naches

Sunnyside

Richland

Kiona

Kennewick

Prosser

Grandview

Mabton

Satus

Zillah

Granger

Roza

Easton

Teanaway

Tieton Rimrock

0 5 10 20 30 MILES

0 5 10 20 30 KILOMETERS

Fault

Fold

EXPLANATION

Structural basins

Roslyn

Kittitas

Selah

Yakima

Toppenish/Benton

Base from U.S. Geological Survey and (or) Washington Division

of Geology and Earth Resources digital data, 1983, 1:100,000

Universal Transverse Mercator projection, Zone 10

Horizontal Datum: North American Datum of 1927 (NAD 27)

Moxee

Ahtanum

Syncline

Toppenish Ridge

Wahluke Syncline

Cold Creek Syncline

Kittitas Valley Syncline

Naneum Ridge Anticline

Manastash Anticline

Umtanum Anticline

Horse Heaven Hills Structure

Ainsley Canyon Anticline

Laurel Hill Monocline

Taneum Monocline

Cleman Mtn. Anticline

Ahtanum Ridge

Wapato Syncline

Wenas Syncline

Rattlesnake

Hills

Structure

Yakima Ridge

Hog

Ranch Structure

Medicine Valley Syncline

Wilson Creek Monocline

Reeser Canyon

Monocline

Naneum Creek

Monocline

Section of the Yakima basin 

where particle tracking was done

Figure 2. Location of structural basins and geologic structures, Yakima River Basin, Washington. (From Ely and others, 

2011.)
Introduction 5

tac15-1019_fig03

Kiona

Naches

Yakima

Wapato

Prosser

Cle Elum

Sunnyside Richland

Union Gap

Ellensburg

Toppenish

121°00'

47°

00'

47°

30'

46°

00'

121°30' 30' 120°00' 119°30'

30'

0 5 10 20 30 40 MILES

0 5 10 20 30 40 KILOMETERS

EXPLANATION

Simplified surficial geology

Quaternary flood deposits

Quaternary deposits, loess

Quaternary and Pliocene volcanic rock

Columbia River Basalt Group rock

Miocene and older volcanic rock

Tertiary granitic and intermediate intrusive rock

Non-marine sedimentary rock

Marine sedimentary rock

Pre-Tertiary metamorphic and intrusive rock

Section of the Yakima basin 

where particle tracking was done

Base from U.S. Geological Survey and/or

Washington Division of Geology and Earth Resources

digital data, 1983, 1:100,000

Universal Transverse Mercator projection, Zone 10

Horizontal Datum: North American Datum of 1927 (NAD 27)

Figure 3. Simplified surficial geology in the Yakima River Basin, Washington. (From Fuhrer and others, 1994.)
6 Particle Tracking for Selected Groundwater Wells in the Lower Yakima River Basin, Washington

Irrigated agriculture has been the primary economic 

activity in the basin since the first non-Indian settlers arrived 

in 1848. Significant development of irrigation infrastructure 

began in 1880 and, by 1902, about 120,000 acres were being 

irrigated, already significantly altering the local groundwaterflow system. As of 2010, about 80 percent of the land area in 

the basin was used for agriculture, including apples, pears, 

cherries, peaches, vegetables, hay, mint, and hops (Washington 

State Department of Ecology, 2010). Dairy operations 

were greatly expanded starting in the 1980s, and are now 

(2015) mostly concentrated around the cities of Sunnyside, 

Grandview, Mabton, and Granger. Current exemptions in State 

water permitting rules place no restrictions on the amount of 

water that can be withdrawn by a well owner for the purposes 

of raising livestock (RCW 90.44.050), so areas with high 

concentrations of Confined Animal Feeding Operations may 

be contributing returned water as groundwater recharge at 

rates greatly exceeding natural recharge for their proportional 

area representation. At present (2015), there are few 

regulations governing how dairy operators must construct or 

manage manure lagoons and no oversight of how animal waste 

is distributed to surrounding fields (173-303 WAC), leading 

to concerns that dairy operations are contributing nitrates to 

the groundwater system in areas immediately upgradient of 

shallow residential drinking-water wells. The recent growth of 

the local dairy industry adds a new potential source for nitrate 

in groundwater, in addition to the application of fertilizers 

during the 135-year local history of irrigated agriculture. Prior 

to the completion of the Yakima River Basin groundwater-flow 

model, there was no rigorous quantitative means to analyze 

hydrologic connections between specific aquifer-recharge 

areas and groundwater wells.

Simulation of Particle Tracking

MODPATH (Pollock, 1994, 2012) is a particle-tracking 

post-processing package that was developed to compute 

three-dimensional flow paths using output from steady-state 

or transient groundwater-flow simulations by MODFLOW, 

the U.S. Geological Survey numerical (finite-difference) 

groundwater-flow model. In this study, MODPATH version 

5.0 was used for the backtracking of particles in time based 

on a previously calibrated transient MODFLOW-2000 

(Harbaugh and others, 2000) groundwater-flow model 

of the Yakima River Basin (Ely and others, 2011). The 

regional Yakima River Basin groundwater-flow model was 

constructed for the purposes of evaluating regional scale 

water-management scenarios and surface-water exchanges, not 

analyses of local-scale impacts at individual wells, so careful 

and defensible interpretation of the particle-tracking results 

requires a detailed understanding of how this model was 

constructed.

Stress Periods

The transient simulation period extends from October 1, 

1959, to September 30, 2001, for a total of 42 water years 

(1960–2001) and 504 monthly stress periods, within each 

of which specified model stress data such as precipitation 

or pumping are held constant. The number of time steps in 

a stress period is equal to the number of days in the month. 

The monthly stress period and daily time step create efficient 

model run times while still capturing relevant flow-system 

characteristics such as seasonal variations in recharge that 

reflect precipitation and irrigation patterns. Recharge derived 

from precipitation primarily occurs in the winter at high 

elevations, and recharge derived from irrigation primarily 

occurs in the summer on the floors of the sedimentary basins.

The particle backtracking analysis specified particle 

“starting” locations at the end of the simulation period, and 

calculated particle movement backwards in time through 

all 504 monthly stress periods until October 1959 or until 

reaching a model boundary. Over this period of time, changes 

to land use, climatic variations, and rapid expansion of 

groundwater use all significantly influenced groundwaterflow directions. The particle pathlines accordingly are 

representative of the flow histories of water withdrawn from 

the sampled wells in September 2001. A more accurate 

representation of flowpaths to locations sampled in 2011 

would require extending the flow model simulation period 

to include an additional decade beyond 2001, but the results 

of such an extension are unlikely to qualitatively alter the 

findings because infrequent installation of new groundwater 

pumping wells during the 2001–2011 decade would limit 

additional pumpage data to use in the extended model 

simulation. The primary benefit of such an extension would 

be to determine the fate of some additional particles that did 

not reach land surface during the 42-year backward tracking 

simulation period.

Discretization

In MODFLOW, the groundwater-flow system is 

subdivided laterally and vertically into rectilinear blocks 

referred to as model cells. The hydraulic properties of the 

material in each cell are assumed to be homogeneous.

The entire Yakima River Basin was simulated using 

a model grid of 600 columns and 600 rows with 1,000-ft2

horizontal cells, but pathlines were not simulated in some of 

this active domain because it was external to the Toppenish/

Benton subbasin (fig. 4). No cells included in the previously 

calibrated and published model were deactivated for 

any simulation, to avoid any possible disruption of local 

groundwater head gradients. The part of this grid relevant to 

particle tracking was only about 300 columns and 200 rows. 

Vertically, the model domain was subdivided 

into 24 model layers, scaled vertically so model 

cells become increasingly thicker with depth (fig. 5). 
Simulation of Particle Tracking 7

tac15-1019_fig04 1020304050

60708090

100200300400500

10

600

1020304050607080901020304050607080901020304050607080901020304050607080902030405060708090

90

80

70

60

50

40

30

20

90

80

70

60

50

40

30

20

90

80

70

60

50

40

30

20

90

80

70

60

50

40

30

20

10

100

10

200 

10

300

10

400

10

500

10

600

90

80

70

60

50

40

30

20

90

80

70

60

50

40

30

20

121°00'

47°

00'

47°

30'

46°

00'

121°30' 30' 120°00' 119°30'

30'

0 5 10 20 30 40 MILES

0 5 10 20 30 40 KILOMETERS

COLUMBIA

RIVER

Base from U.S. Geological Survey and/or

Washington Division of Geology and Earth Resources

digital data, 1983, 1:100,000

Universal Transverse Mercator projection, Zone 10

Horizontal Datum: North American Datum of 1927 (NAD 27)

82

22

22

97

97

97

12

12

82

82

240

24

410

240

90

90

90

24

410

12

97

12

241

82

82 12

CLEMAN

MOUNTAIN

H O R S E

H E A V E N

H I L L S

U M T A N U M

R I D G E

Y A K I M A

 V A L L E Y

Y A K I M A R I D G E

UNION GAP

AHTANUM RIDGE

COWICHE MTN

TOPPENISH RIDGE

K I T T I T A S

V A L L E Y

M A N A S T A S H

R I D G E

NORTH RIDGE

AMERICAN RIDGE

C

A

S C A

D E R

A

N

G

E

KACHESS RIDGE

TEANAWAY

CLE ELUM

RIDGE

W E N AT C H E E

M O U N TA I N S

RIDGE

Hanford Site

Snipes

Mountain

R A T T L E S N A K E H I L L S

Cle Elum

Parker

Wapato

Ellensburg

Umtanum

Kittitas

Thrall

Toppenish

Union Moxee City Gap

Yakima

Selah

Naches

Sunnyside

Richland

Kiona

Kennewick Prosser

Grandview

Mabton

Satus

Zillah

Granger

Roza

Easton

Teanaway

Tieton

Rimrock

KITTITAS COUNTY

YAKIMA COUNTY

COUNTY

BENTON COUNTY

KLICKITAT

COUNTY

YAKIMA

Bumping

Lake

Bumping River

Tieton River

Dry Creek

RIVER

YAKIMA

Cle Elum Lake

Rimrock

Lake

Kachess

Lake Keechelus Lake

Teanaway

River

Little Naches River River

Naches

RIVER

YAKIMA

RIVER

YAKIMA

Cl e Elum

River

EXPLANATION

Active model domain

10

11

101120

1

1

20

Figure 4. Location and extent of the groundwater-flow model grid used to simulate the Yakima River Basin aquifer system, 

Washington. The inset depicts the detailed horizontal discretization for the first 20 rows and columns of the model domain. 

(Modified from Ely and others, 2011.)
8 Particle Tracking for Selected Groundwater Wells in the Lower Yakima River Basin, Washington

tac15-1019_fig05a

Parker

Wapato

Umtanum

Thrall

Toppenish

Moxee

Union Gap City

Yakima

Selah

Naches

Sunnyside

Grandview

Mabton

Satus

Zillah

Granger

Roza

Tieton

Ahtanum Creek

Tieton

River

Wenas Creek

North Fork

Dry Creek

Creek

Toppenish

Satus Creek

Umtanum Creek

Cowichee

Creek

Cowichee

Creek

River

Naches

RIVER

YAKIMA

No . Fork

COLUMBIA

RIVER

KITTITAS COUNTY

YAKIMA COUNTY

KLICKITAT COUNTY

YAKIMA COUNTY

0 5 10 20 MILES

0 5 10 20 KILOMETERS

A A'

121°00' 30' 120°00'

46°

30'

46°

00'

EXPLANATION

Basin-fill deposits

Hydrogeologic Unit

Line of hydrogeologic

section

A A'

HGU

Section

location

Base from U.S. Geological Survey and/or Washington Division of Geology and Earth Resources digital data, 1983, 1:100,000

Universal Transverse Mercator projection, Zone 10. Horizontal Datum: North American Datum of 1927 (NAD 27)

A

Figure 5. Location of model cross section through the Toppenish/Benton subbasin (A) and hydrogeologic section 

showing relation of model hydrogeologic units to model computational layers (B), Yakima River Basin aquifer 

system, Washington. (Modified from Ely and others, 2011.)
Simulation of Particle Tracking 9

tac15-1019_fig05b

-3,500

-4,000

-3,000

-2,500

-2,000

-1,500

-1,000

1000

FEET FEET

WEST EAST WEST EAST

METERS

-500

500

0 0

-4,500

-5,000

-1,300

-1,400

-1,500

-1,200

-1,100

-1,000

-900

-800

-700

-600

-500

-400

-300

-200

200

300

100

-100

-3,500

-4,000

-3,000

-2,500

-2,000

-1,500

-1,000

1000

-500

500

0

-4,500

-5,000

A A' A A'

1

2

3

4

5

6

7

10

11

12

13

14

15

16

17

18

19

20

21

22

23

8

9

24

Model

layer

HGU 10

HGU 11

HGU 12

HGU 13

HGU 14

HGU 21

HGU 22

HGU 23

HGU 24

HGU 25

HGU 26

HGU 27

HGU 28

HGU 29

HGU 30

HGU 31

HGU 32

HGU 33

HGU 34

HGU 35

HGU 36

HGU 37

HGU 38

HGU 39

HGU 40

HGU 41

HGU 42

HGU 43

HGU 44

HGU 45

HGU 46

Overburden Saddle Mountain Basalt

Maston

Interbed Wanapum Basalt

Vantage

Interbed Grande Ronde Basalt

0 5 10 20 MILES

0 5 10 20 KILOMETERS

Bedrock

B

Figure 5.—Continued
10 Particle Tracking for Selected Groundwater Wells in the Lower Yakima River Basin, Washington

Total model domain thickness ranged from 1,780 to 8,800 ft. 

The percentage of thickness per layer increased parabolically 

with depth. Each cell in the top layer was specified as 1.00 

percent of the total thickness at that model row and column. 

Each successively lower layer was specified as a slightly 

larger percentage of total thickness, with the lowest and 

thickest layer, 24, being 9.82 percent of the total thickness. 

Mean layer thickness is 43 ft (maximum 88 ft) for layer 1 and 

419 ft (maximum 866 ft) for layer 24 (table 1). The paraboliclayering approach was used to provide relatively finer vertical 

discretization for shallow layers of the model (Ely and others, 

2011). Hydraulic property values for each hydrogeologic unit 

in the model were calibrated and then translated into effective 

properties for each of the 24 computational layers using the 

Hydrogeologic-Unit Flow package in MODFLOW-2000 

(Anderman and Hill, 2000). 

Table 1. Descriptive statistics for the thickness of groundwater 

model layers used to simulate the Yakima River Basin aquifer 

system, Washington.

[From Ely and others (2011)]

Model 

layer

Percentage of 

total model 

thickness

Layer thickness, in feet

Minimum Mean Maximum

1 1.00 18 43 88

2 1.06 19 45 93

3 1.15 20 49 101

4 1.26 22 54 111

5 1.41 25 60 124

6 1.59 28 68 140

7 1.80 32 77 159

8 2.03 36 87 179

9 2.30 41 98 203

10 2.59 46 110 228

11 2.92 52 125 257

12 3.27 58 139 288

13 3.66 65 156 323

14 4.07 72 174 359

15 4.51 80 192 398

16 4.98 89 212 439

17 5.48 97 234 483

18 6.02 107 257 531

19 6.58 117 281 580

20 7.17 127 306 632

21 7.79 139 332 687

22 8.43 150 359 743

23 9.11 162 388 803

24 9.82 175 419 866

Model Inputs

Boundary conditions define the locations and manner 

in which water enters and exits the active model domain. 

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 surface-water irrigation, and (2) exits the system as 

streamflow, evapotranspiration, and groundwater pumpage. 

The specified boundaries of the model coincide as much as 

possible with natural hydrologic boundaries. Four types of 

model boundaries were used: 

1. No-flow boundaries (groundwater divides),

2. Head-dependent flux boundaries (drains, streams, 

reservoirs),

3. Specified-head boundaries, and

4. Specified-flux boundaries (recharge, pumpage, and 

streamflow inflow).

No-flow boundary conditions were specified at the major 

topographic divides that coincide with the lateral model 

boundaries other than the Columbia River, which are assumed 

to be groundwater divides. The topographic divides are either 

exposed bedrock or basalt. These divides are the crest of the 

Cascade Range to the west, the drainage-basin divide of the 

Yakima and Wenatchee Rivers to the north, and Horse Heaven 

Hills to the south (fig. 6).

Surface-water features including streams and canals were 

simulated with the MODFLOW Drain (DRN) (Harbaugh and 

others, 2000) or Streamflow Routing (SFR) packages (Prudic 

and others, 2004; Niswonger and Prudic, 2005). Drains 

were used primarily in the headwaters of the Yakima River 

or for streams that either were ephemeral or had very small 

annual discharge. This approach added numerical stability to 

the model in areas far from regions of greatest interest and 

accounted for the generally gaining upland stream reaches. 

The simulated quantity of water exiting the system at a drain 

cell (MODFLOW only allows simulated groundwater flow 

into a drain cell) is equal to the product of a user-specified 

drain conductance and the difference between the simulated 

hydraulic head in the drain cell and the specified elevation 

of the drain (stream stage). Drain elevation was set at land 

surface (sampled from the U.S. Geological Survey [USGS] 

10-meter Digital Elevation Model) minus 10 ft. 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 unavailable. Drain conductances were divided 

into groups based on drain size, with values determined by 

the iterative model calibration process to most accurately 

reflect measured groundwater heads and streamflows (Ely and 

others, 2011).
Simulation of Particle Tracking 11

tac15-1019_fig06

121°00'

47°

00'

47°

30'

46°

00'

121°30' 30' 120°00' 119°30'

30'

0 5 10 20 30 40 MILES

0 5 10 20 30 40 KILOMETERS

COLUMBIA

RIVER

Base from U.S. Geological Survey and/or

Washington Division of Geology and Earth Resources

digital data, 1983, 1:100,000

Universal Transverse Mercator projection, Zone 10

Horizontal Datum: North American Datum of 1927 (NAD 27)

CLEMAN

MOUNTAIN

H O R S E

H E A V E N

H I L L S

U M T A N U M

R I D G E

Y A K I M A

 V A L L E Y

Y A K I M A R I D G E

UNION GAP

AHTANUM RIDGE

COWICHE MTN

TOPPENISH RIDGE

K I T T I T A S

V A L L E Y

M A N A S T A S H

R I D G E

NORTH RIDGE

AMERICAN RIDGE

C

A

S C A

D E R

A

N

G

E

Gilbert

Peak

Snoqualmie Pass

Van Epps

Pass

Esmeralda

Peaks KACHESS RIDGE

TEANAWAY

Deception

Pass

CLE ELUM

RIDGE

W E N AT C H E E

M O U N TA I N S

RIDGE

Hanford Site

Bumping

Lake

Bumping River

Tieton River

Dry Creek

RIVER

YAKIMA

Cle Elum

Lake

Rimrock

Lake

Kachess

Lake Keechelus Lake

Teanaway

River

Little Naches River River

Naches

RIVER

YAKIMA

RIVER

YAKIMA

Cl e Elum

River

EXPLANATION

Groundwater-flow model

boundary conditions

USGS gaging station and No.

Streamflow-routing cell

General-head cell

Drain cell

Specified-head cell

12472800

12472800

Figure 6. Location of the groundwater-flow model boundary conditions for drains, general and specified heads, and streamflow 

routing cells, Yakima River Basin aquifer system, Washington. (Modified from Ely and others, 2011.)
12 Particle Tracking for Selected Groundwater Wells in the Lower Yakima River Basin, Washington

Excluding drain cells described above, the Yakima River 

and its tributaries were simulated using the MODFLOW 

Streamflow-Routing (SFR2) package (Niswonger and 

Prudic, 2005) to route streamflow and calculate river-aquifer 

exchanges. The model has 250 SFR2 segments and 8,533 

uniquely simulated reaches, and exchange of water between 

streams and the shallow aquifer system is controlled by the 

differences between groundwater levels and stream stage for 

each reach. Some inaccuracy was introduced in the simulation 

of groundwater flow to and from the streams by using 

average stream stages and simulating average groundwater 

elevations within model cells. This uncertainty was not 

considered a problem in areas of gentle relief, such as the 

areas containing the groundwater sampling locations used for 

particle tracking in this analysis, but may contribute to model 

errors in places where streams are deeply incised. For the 

Yakima River, stream depth was computed using Manning’s 

equation assuming an eight-point cross section. Average depth 

and width for the cross sections 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).

Subsurface exchange with the five managed reservoirs 

in the basin was simulated using the MODFLOW generalhead boundary (GHB) package, using stages derived from 

USGS 1:24,000-scale topographic maps and held constant 

throughout the simulation period. This approach neglects the 

changes in reservoir exchange due to seasonally managed 

reservoir drawdowns, but the effect is localized to the 

high-elevation sites near the reservoirs and should have no 

measureable effect on particle-tracking simulations in the 

Toppenish/Benton subbasin.

The MODFLOW time variant specified-head (CHD) 

package was used to simulate the Columbia River on the 

eastern boundary of the model (fig. 6). Monthly streamflow 

for the USGS streamgage at Columbia River below Priest 

Rapids Dam (12472800) during the simulation period ranged 

from 52,500 to 461,000 ft3/s. Assigned specified head at each 

cell was based on the stage-discharge relation for the Hanford 

Reach of the Columbia River at mean annual flow simulated 

with a two-dimensional hydraulic model (T. Waddle, U.S. 

Geological Survey, written commun., 2008). Upstream of 

the several dams along this boundary, head was specified 

as a constant value equal to the mean pool elevation. Mean 

differences between monthly pool elevations over the 

simulation period were less than 5 ft, and most differences 

were less than 3 ft.

Groundwater pumping in the Yakima River Basin was 

estimated as part of this study for eight categories of use 

during 1960–2001 (Vaccaro and Sumioka, 2006). The eight 

categories of pumping were:

1. Public water supply (including wells for Group 

B systems, defined as supplies with less than 15 

connections),

2. Self-supplied domestic (permit-exempt wells), 

3. Irrigation,

4. Frost protection,

5. Livestock and dairy operations,

6. Industrial and commercial processes, 

7. Fish and wildlife propagation, and

8. Groundwater claims.

Methods of data collection, pumpage estimates, reliability of 

the estimates, and a comparison with appropriated quantities 

are described by Vaccaro and Sumioka (2006) and Vaccaro 

and others (2009). By the end of calendar year 1960, total 

annual pumpage in the basin, excluding standby/reserve 

pumpage, was about 116,000 acre-ft (160 ft3/s). By 2000, 

total annual pumpage (excluding standby/reserve rights) was 

estimated to be about 317,000 acre-ft (438 ft3/s). Irrigation 

accounted for about 60 percent of the pumpage, followed 

by public water supply at about 12 percent. The smallest 

category of pumpage was for livestock and dairy operations, 

with pumpage in 2000 estimated to be about 6,700 acre-ft 

(9.2 ft3/s). However, total pumpage for this category is 

not well understood because under Washington law, it is 

permissible to pump any quantity necessary to meet stockwatering needs (RCW 90.44.050), and no information is 

available on the quantities withdrawn for this purpose. For 

this reason, local livestock and dairy operations groundwater 

use may be underestimated by an unknown amount. Pumping 

locations were simulated in the model using the MODFLOW 

Well (WEL) package, which assigns volumetric fluxes to a 

computational model cell and not a specific point within that 

cell. Most of the wells used to determine starting locations 

for the particle backtracking analysis were domestic wells not 

large enough to have been included in the model as unique 

pumping centers. Such wells were grouped together by quarter 

section and simulated as a single point withdrawal. As a result 

of this approach, the particle pathlines presented here are 

indicative of background groundwater-flow vectors and do not 

reflect the particle pathline divergence that would be expected 

for backtracking particles from within a pumping center cone 

of depression.

Surface recharge to the aquifer system from precipitation 

and irrigation was simulated with the MODFLOW Recharge 

(RCH) package. The values used were derived from the 

Precipitation-Runoff Modeling System (PRMS) (Leavesley 

and others, 1983) and Deep Percolation Model (Bauer and 

Vaccaro, 1987; Vaccaro, 2007). The two models were used to 

simulate two different hydrologic regimes of the Yakima River 

Basin (Vaccaro and Olsen, 2007), and to estimate groundwater 

recharge for predevelopment and 2007 land-use and landcover conditions for those regimes (fig. 7). Seasonal variations 

in recharge due to both precipitation and agricultural 

applications result in annual changes in groundwater-flow 

velocities and directions in the shallow sediments of the lower 

Yakima River Basin.
Simulation of Particle Tracking 13

tac15-1019_fig07

121°00'

47°

00'

47°

30'

46°

00'

121°30' 30' 120°00' 119°30'

30'

0 5 10 20 30 40 MILES

0 5 10 20 30 40 KILOMETERS

COLUMBIA

RIVER

Base from U.S. Geological Survey and/or

Washington Division of Geology and Earth Resources

digital data, 1983, 1:100,000

Universal Transverse Mercator projection, Zone 10

Horizontal Datum: North American Datum of 1927 (NAD 27)

CLEMAN

MOUNTAIN

H O R S E

H E A V E N

H I L L S

U M T A N U M

R I D G E

Y A K I M A

 V A L L E Y

Y A K I M A R I D G E

UNION GAP

AHTANUM RIDGE

COWICHE MTN

TOPPENISH RIDGE

K I T T I T A S

V A L L E Y

M A N A S T A S H

R I D G E

NORTH RIDGE

AMERICAN RIDGE

C

A

S C A

D E R

A

N

G

E

KACHESS RIDGE

TEANAWAY

CLE ELUM

RIDGE

W E N AT C H E E

M O U N TA I N S

RIDGE

Hanford Site

Snipes

Mountain

R A T T L E S N A K E H I L L S

KITTITAS COUNTY

YAKIMA COUNTY

COUNTY

BENTON COUNTY

KLICKITAT

COUNTY

YAKIMA

Bumping

Lake

Bumping River

Tieton River

Dry Creek

RIVER

YAKIMA

Cle Elum

Lake

Rimrock

Lake

Kachess

Lake Keechelus Lake

Teanaway

River

Little Naches River River

Naches

RIVER

YAKIMA

RIVER

YAKIMA

Cl e Elum

River

Less than 0.5

25 to 48

Greater than 48

Mean annual recharge, in inches 

13 to 24

7 to 12

0.6 to 2

3 to 6

EXPLANATION

Figure 7. Spatial distribution of mean annual recharge simulated for 1960–2001, Yakima River Basin aquifer system, Washington. 

Recharge estimates for earlier years reflected historical land-use patterns. (Modified from Ely and others, 2011.)
14 Particle Tracking for Selected Groundwater Wells in the Lower Yakima River Basin, Washington

Hydrologic properties including hydraulic conductivity 

and storage coefficients are documented in Ely and others 

(2011), but one additional model property of specific 

importance to particle-tracking simulations is the use of 

hydrologic flow barriers (HFB) implemented with the 

MODFLOW HFB package (Hsieh and Freckleton, 1993). 

Previous work indicates that the lateral hydraulic conductivity 

of members of the Columbia River Basalt Group is reduced 

by 1–2 orders of magnitude in areas of intense folding and 

faulting, especially fault-associated anticlines such as the 

Yakima Fold and Thrust Belt (Hansen and others, 1994; 

Packard and others, 1996; Reidel and others, 2002). The 

HFB package was used to simulate flow barriers for various 

geologic structures, including Ahtanum and Toppenish 

Ridges, Rattlesnake Hills, and the Horse Heaven Hills, all 

of which influence groundwater-flow directions near the 

Toppenish/Benton subbasin (fig. 8). The degree to which 

such structures were considered significant to groundwater 

flow was determined by analysis of measured water-level 

variations across each structure, in conjunction with the 

results of model calibrations that included or excluded such 

effects. Barriers to flow were not simulated as continuous in 

cases where seasonal changes to groundwater head suggested 

a hydraulic connection. Where buried beneath basin-fill 

sediments, the structures were assumed to be absent in the 

overlying sediments and no barriers were simulated, so such 

barriers only affected particle pathlines in cases where such 

pathlines intersected mapped basalt units. In several cases, the 

simulated particle pathlines make sharp changes in direction 

around the edges of mapped HFBs, despite poor constraints 

on the known locations of these features at depth. For this 

reason, particle pathlines that appear to have sharp corners 

are likely to be artifacts of a model construct that represents a 

simplified version of the complex natural system, and may not 

be indicative of true particle histories.

Calibration

The model was calibrated with the iterative parameter 

estimation software package Parallel PEST (Doherty, 2010) 

using 20,279 water-level measurements from 2,196 wells 

in the USGS National Water Information System database 

and 7,056 monthly streamflow measurements from 

14 streamgages. Calibrations using groups of hydraulic heads 

and streamflows were done with observation weights adjusted 

to ensure equal contribution by these two groups to the 

model objective function in accordance with the guidelines of 

Doherty and Hunt (2010). Calibrated model parameters and 

their values are shown in table 2, and the comparison between 

measured and simulated groundwater heads is presented in 

figure 9. In the Toppenish Basin, head residuals (simulated 

minus measured head) commonly were less than 10 ft at all 

locations except those locations near the distal edge of the 

basin, adjacent to outcropping basalts. Further information 

about the model calibration procedure, weighting scheme, and 

sensitivity analysis is included in Ely and others (2011).

Particle Tracking

A primary objective of this study was to identify the areal 

extent of the zones of contribution for groundwater wells that 

EPA identified in April 2011 as having exceeded the nitrate 

maximum contaminant level of 10 mg/L. EPA collected 

samples from various groundwater sources in the lower 

Yakima River Valley and identified excess nitrate in samples 

from 121 locations (fig. 10).

In most cases, samples were collected from surface 

access points such as outdoor or kitchen faucets providing 

groundwater from wells without an identified well log that 

specified the depth in the aquifer system from which the 

groundwater was drawn. Because particle-tracking analyses 

are sensitive to the starting depths of the particles, incorrect 

identification of the screened intervals of these wells would 

result in particle pathlines that are consistent with simulated 

flow directions but may be too long (if the simulated well 

depth was too deep) or too short (if the simulated well 

depth was too shallow). In practice, this problem is partially 

mitigated by the assignment of particle starting locations to a 

range of depths determined to be appropriate for each well.

To determine these depths, every well log in the 

Washington State Department of Ecology (2015) well database 

within a 1-mi radius of each recorded sample location was 

compiled. For each sample location, well depths were plotted 

to identify the common depth to which similar wells were 

screened, on the assumption that boreholes drilled to construct 

the wells were not drilled any deeper than was necessary to 

supply sufficient water for the specified use. In cases where 

a sample location could be identified as being derived from a 

specific well in the list of well logs, the actual screened depth 

was used. Assigned well depths for each sample location are 

provided in appendix A.

With screen depths either from a well log or estimated 

based on an average of nearest neighbor values, particle 

starting locations were specified at the edges and center of 

each 1,000-ft computational model cell that contained an 

identified sampling location and the known or estimated depth 

of the well screen. By this process, the 121 identified wells 

were assigned to 89 unique model cells; some cells contained 

more than one well. For each of these 89 rectilinear cells, 

starting particle locations were specified at the middle of each 

cell, at each of the 8 cell corners, at the midpoint of each of 

the 12 cell edges, and at the center of each of the 6 cell faces, 

for a total of 27 particles per cell (9 particles in a 3-by-3 array 

in each of 3 layers). These particle locations encompass the 

entire volume of the model cells containing the groundwater 

that was sampled by EPA. The starting particle locations were 

specified to exist at these cell boundaries at the end of 42-year 

simulation period, and were tracked backwards in time from 

these locations. Thus 27 groundwater flowpaths connecting to 

each of the 89 model cells containing the 121 sample locations 

were simulated, generating 2,403 unique particle tracks. 
Simulation of Particle Tracking 15

tac15-1019_fig08

121°00'

47°

00'

47°

30'

46°

00'

121°30' 30' 120°00' 119°30'

30'

0 5 10 20 30 40 MILES

0 5 10 20 30 40 KILOMETERS

COLUMBIA

RIVER

Base from U.S. Geological Survey and/or

Washington Division of Geology and Earth Resources

digital data, 1983, 1:100,000

Universal Transverse Mercator projection, Zone 10

Horizontal Datum: North American Datum of 1927 (NAD 27)

82

22

22

97

97

97

12

12

82

82

240

24

410

240

90

90

90

24

410

12

97

12

241

82

82 12

R A T T L E S N A K E H I L L S

CLEMAN

MOUNTAIN

H O R S E

H E A V E N

H I L L S

U M T A N U M

R I D G E

Y A K I M A

 V A L L E Y

Y A K I M A R I D G E

UNION GAP

AHTANUM RIDGE

COWICHE MTN

TOPPENISH RIDGE

K I T T I T A S

V A L L E Y

M A N A S T A S H

R I D G E

NORTH RIDGE

AMERICAN RIDGE

C

A

S C A

D E R

A

N

G

E

KACHESS RIDGE

TEANAWAY

CLE ELUM

RIDGE

W E N AT C H E E

M O U N TA I N S

RIDGE

Hanford Site

Cle Elum

Parker

Wapato

Ellensburg

Umtanum

Kittitas

Thrall

Toppenish

Moxee City Union Gap

Yakima

Selah

Naches

Sunnyside Richland

Kiona

Kennewick Prosser

Grandview

Mabton

Satus

Zillah

Granger

Roza

Easton

Teanaway

Tieton Rimrock

KITTITAS COUNTY

YAKIMA COUNTY

COUNTY

BENTON COUNTY

KLICKITAT

COUNTY

YAKIMA

Bumping

Lake

Bumping River

Tieton River

Dry Creek

RIVER

YAKIMA

Cle Elum Lake

Rimrock

Lake

Kachess

Lake Keechelus Lake

Teanaway

River

Little Naches River River

Naches

RIVER

YAKIMA

RIVER

YAKIMA

Cl e Elum

River

Horizontal-flow barrier

characteristic, in feet per day

divided by feet (1/d)

Low 2.0×10-6

2.7×10-5

4.5×10-5

Medium

High

EXPLANATION

Figure 8. Location of model structures and flow barrier conductivity, Yakima River Basin aquifer system, Washington. (Modified 

from Ely and others, 2011.)
16 Particle Tracking for Selected Groundwater Wells in the Lower Yakima River Basin, Washington

Table 2. Calibrated model parameter values for the Yakima River Basin aquifer system, Washington.

[Abbreviations: HGU, hydrogeologic unit; SM, Saddle Mountains unit; WN, Wanapum unit, GR, Grande Ronde; IFZ, interflow zone; INT, flow interior; 

ft/d, foot per day; ft2/d; square foot per day]

Hydrogeolgic unit

Horizontal hydraulic conductivity (ft/d) Vertical 

anisotropy

Specific storage (ft-1)

Minimum Mean Maximum Constant Minimum Mean Maximum Constant

HGU 1 93.37 8.66 1.39E-02

HGU 2 1.24 133.32 3.45E-04

HGU 3 35.73 115.24 2.92E-04

HGU 4 215.71 150.51 4.61E-02

HGU 5 11.13 644.11 4.40E-05

HGU 6 9.67 282.53 1.02E-04

HGU 7 86.44 96.05 8.84E-04

HGU 8 69.21 485.42 4.62E-05

HGU 9 12.89 2,302.48 1.92E-05

HGU 10 91.74 246.19 8.33E-05

HGU 11 1,672.30 959.54 1.56E-03

HGU 12 24.34 11.78 2.22E-04

HGU 13 1.46 1,449.37 3.28E-04

HGU 14 99.92 60.19 1.03E-04

HGU 15 9.25 31.58 1.52E-04

HGU 16 1,490.37 98.53 1.57E-03

HGU 17 4.19 45.83 3.02E-04

HGU 18 36.07 104.19 1.04E-04

HGU 19 2.51 76.11 8.77E-05

HGU 20 98.24 106.43 1.15E-04

SM_IFZ 35.78 119.15 261.43 4,336.83 1.51E-04 3.39E-04 3.73E-04

SM_INT 4.91E-04 1.66E-03 3.56E-03 0.19 1.68E-05 1.96E-05 2.15E-05

HGU31 0.36 0.88 1.57 157.39 7.70E-03

WN_IFZ 25.67 129.94 278.32 530.33 2.83E-04 5.99E-04 6.76E-04

WN_INT 6.92E-05 3.51E-04 7.49E-04 1.66 2.68E-05 5.60E-05 6.33E-05

HGU38 9.31 461.85 3.22E-04

GR_IFZ 4.28 22.99 90.97 2,177.27 1.11E-04 2.19E-04 2.65E-04

GR_INT 7.61E-05 3.95E-04 1.54E-03 0.32 2.43E-05 4.64E-05 5.76E-05

HGU 47 24.88 112.90 1.35E-04

HGU 48 0.01 0.03 1.09 1.21 1.03E-04

Streamflowrouting cells Conductivity (ft/d)

0.45

0.75

Yakima River

Tributaries

Drain cells Conductance (ft2

/d)

Headwater streams 323,166.00

Hydraulicflow barrier Hydraulic characteristic (1/d)

Low 5.18E-06

3.59E-05

4.20E-04

Medium

High
Simulation of Particle Tracking 17

tac15-1019_fig09

Weighted simulated values, in feet

0

-5,000 0 5,000 10,000 15,000 20,000 25,000 30,000 35,000 40,000 45,000

5,000

-5,000

10,000

15,000

20,000

25,000

30,000

35,000

40,000

45,000

Weighted measured values

y = 1.02x + 88.8

Figure 9. Weighted measured values as a function of weighted simulated values, 

Yakima River Basin aquifer system, Washington. (From Ely and others, 2011.)

In the absence of any reliable information about the 

porosity of the sedimentary deposits simulated, MODPATH’s 

default porosity value of 0.01 was used. The uncertainty 

associated with this value is large, and the time of travel 

computed from the simulated pathlines may differ from actual 

time of travel for the same particle starting locations.

MODPATH implements both strong and weak sinks as 

potential exit points for simulated particles. In the case of 

weak sinks, some particles can be allowed to pass through a 

model cell containing a pumping well that does not capture 

all water in that model cell, but in the particle backtracking 

analysis all particles were started in model cells containing 

the pumping well and then backtracked from there. The 

MODPATH input files were constructed to specify that 

particles were allowed to pass through weak sinks.

Over the duration of the simulation, most particles 

eventually intersected the water table and were thus 

considered to have entered the model as recharge at the water 

table in less than the 42-year simulation period (fig. 11). Using 

the aquifer properties described above, the average time-oftravel for all particles was about 10 years and 3 months. This 

result does not suggest that the contaminated wells of interest 

are currently pumping water with this age, because some of 

the particles had not yet reached the water table or other model 

boundary within the 42-year simulation duration. Although 

travel times for some individual particles from shallow 

starting locations were as short as 22 days, the 1,950-day 

average time-of-travel for particles that did reach the water 

table suggests that at least 5 years should be expected to elapse 

before any modifications to current land‐use practices that 

might change nitrate loading rates would manifest as changes 

in measured nitrate concentrations at the average well location. 
18 Particle Tracking for Selected Groundwater Wells in the Lower Yakima River Basin, Washington

tac15-1019_fig10

Sunnyside

Grandview

Zillah

Toppenish

Granger

Wapato

Mabton

Harrah

Yakima River

097

082

223

022

241

120°00' 10' 20' 120°30'

46°

30'

46°

20'

46°

10'

0 4 6 8 2 MILES

0 4 8 6 2 KILOMETERS

Base modified from U.S. Geological Survey data and other digital sources. Lambert Conformal Conic projection 

Horizontal Datum: North American Datum of 1983 (NAD 83)

300

400

310

390

380

370

360

350

340

330

320

290

280

270

260

250

410

420

410

420430440450460

470

480490

500

510520

530

Column

Row

EXPLANATION

Model grid (10 × 10 cells)Model cell used in particle 

tracking analysis

Figure 10. Locations of 121 groundwater samples with nitrate concentrations exceeding 10 milligrams per liter, lower Yakima River Valley, Washington, 2008–2011.
Simulation of Particle Tracking 19 tac15-1019_fig11

Sunnyside

Union Gap

Prosser

Grandview

Zillah

Toppenish

Moxee

Granger

Wapato

Mabton

Harrah

Yakima River 097

082

024

023

022

024

200

300

400

310

390

380

370

360

350

340

330

320

210

290

280

270

260

250

240

230

220

190

180

170

410

460

450

440

430

420

390

400

410420430440450460470480490

500

510520530540550560570580590

Base modified from U.S. Geological Survey data and other digital sources. Lambert Conformal Conic projection 

Horizontal Datum: North American Datum of 1983 (NAD 83)

120°00' 20' 120°40'

46°

20'

46°

00'

0 4 6 8 2 MILES

0 4 8

10

10 6 2 KILOMETERS

EXPLANATION

Model grid (10 × 10 cells)Water-particle pathline—Shows 

path of groundwater flow. Each line color represents different source domestic well

Column

Row

Figure 11. Simulated particle histories for particles that ended the 42-year simulation period at the location of model cells containing selected groundwater wells in the lower Yakima River Basin, Washington.
20 Particle Tracking for Selected Groundwater Wells in the Lower Yakima River Basin, Washington

Ahtanum Creek

Tuyuck Creek

Red Willow

Creek

120°42' 40' 38' 36' 120°34'

46°

30'

28'

46°

26'

Base modified from U.S. Geological Survey data and other digital sources. 

Lambert Conformal Conic projection 

Horizontal Datum: North American Datum of 1983 (NAD 83)

EXPLANATION

Model grid

Model cell used in particle tracking 

analyses

Water-particle pathline—Shows path 

of groundwater flow. Each line 

color represents different source 

domestic well

Hydrologic flow barrier

Pathline endpoint—Particle tracking 

ending location

0 1 2 MILES

0 1 2 KILOMETERS

210

220230240

400

410

430

420

Column

Row

A

Figure 12. Modeled structural barriers to groundwater flow in some locations cause particle flow paths to change direction 

sharply (A) and seasonally varying recharge patterns and associated changes in hydraulic head gradients cause some particle 

paths to temporarily reverse direction (B), Yakima River Basin, Washington.

By contrast, the 323 particles that did not reach a model 

boundary in the 42-year simulation period are associated with 

groundwater sampling locations for which changes to nitrate 

loading in the contributing recharge area are decades away 

from having a measurable effect at the point of withdrawal, 

and at which future nitrate concentrations may continue 

to increase as a result of the increasing trend of nutrient 

application rates over the last 42 years. 

Long delays between changes in nitrate loading rates and 

measurable changes in nitrate concentrations at distant 

groundwater sampling locations do not imply an aquifer 

has not yet been affected, only that the changes have not yet 

propagated to the sampling location. Particle tracking results 

are summarized in appendix A.

In some locations, structural features such as anticlinal 

ridges present barriers to groundwater flow, causing sharp 

changes in hydraulic gradients and groundwater flow 

directions (fig. 12A). For areas with low hydraulic gradients 

and high recharge rates derived from seasonal irrigation, 

fluctuating groundwater levels may cause flow directions and 

particle pathlines to nearly reverse direction on an annual 

timescale (fig. 12B).
Simulation of Particle Tracking 21

tac15-1019_fig12b

119°55' 54' 53' 52' 51' 50' 119°49'

21'

46°

22'

20'

46°

19'

0 1 MILE

0 1 KILOMETER

Base modified from U.S. Geological Survey data and other digital sources. 

Lambert Conformal Conic projection 

Horizontal Datum: North American Datum of 1983 (NAD 83)

EXPLANATION

Model grid

Model cell used in particle 

tracking analyses

Water-particle pathline—Shows 

path of groundwater flow. 

Each line color represents 

different source domestic well

Pathline endpoint—Particle 

tracking ending location

B

410420430

470

460

Column

Row

Figure 12.—Continued

In all cases, particle backtracking results record the 

time for groundwater to travel through the saturated zone to 

the model cell containing a groundwater well, and neglect 

additional time that may be required for water to percolate 

through the unsaturated zone between land surface and the 

water table. Within the area of interest in the Yakima Basin, 

unsaturated zone transport is reasonably assumed to be 

vertical and relatively brief in irrigated areas with shallow 

water tables.

Delineation of Contributing Areas

Particle pathlines from each of the 89 individual 

computational cells containing at least one nitrate-affected 

well were traced backward by following transient 

groundwater-flow directions until they reached the water 

table or other model boundary or the model simulation period 

ended (fig. 13). For cells with particles that terminated at the 

water table, the simulated particle ending locations outline the 

approximate area of land surface that contributed water to that 

cell at the end of the simulation period (fig. 13A). This area 

typically is small for cases where the cells contained wells 

in which simulated pumping rates were too small to cause 

significant simulated drawdowns throughout the entire model 

cell. The absence of a cone of depression around the particle 

starting locations prevents divergence of the particle pathlines. 

For this reason, the reported areas contributing recharge 

should be interpreted as the calculated centroid of the true 

capture area, which may be larger than reported here.

Sources of nitrate loading to the aquifer that are not 

applied at land surface might contribute to the contaminated 

wells if they are located anywhere along the particle pathlines 

(fig. 13B). Examples of potential non-surface nutrient loads 

include septic systems or improperly sealed or abandoned 

wellheads that allow surface contamination to rapidly migrate 

to relevant aquifer depths. Such sources, if present, would 

expand the area contributing recharge to any specified model 

cell containing a well to include the entire area between the 

well location and the most distant particle endpoints, referred 

to as the areal extent of the zone of contribution (Alley, 2003).

Some parts of the study area include a high density of 

sampled wells. In those areas, such as northeast of Granger, 

overlapping particle pathlines are difficult to interpret and 

have been separated into pathlines representing multiple 

groups of wells for the purpose of easing visual interpretation 

(figs. 14A and 14B). 

Pathlines originating from model cells near the town of 

Mabton (fig. 11) show significant divergence depending on 

depth. Groundwater flow directions in the shallow basin-fill 

sediments are from the basin margins toward the Yakima 

River, but flow directions in the deeper basalt units are 

controlled by groundwater head gradients at depth rather 

than local surface topography, resulting in flowpaths that 

appear to cross under the Yakima River as they flow north to 

contributing recharge areas east of Sunnyside.
22 Particle Tracking for Selected Groundwater Wells in the Lower Yakima River Basin, Washington

tac15-1019_fig13a

A

Zillah

Toppenish

Yakima River

082

Base modified from U.S. Geological Survey data and other digital sources. 

Lambert Conformal Conic projection 

Horizontal Datum: North American Datum of 1983 (NAD 83)

120°18' 17' 16' 15' 14' 13' 120°12'

46°

28'

27'

26'

25'

24'

46°

23'

0 1 MILE

0 1 KILOMETER

EXPLANATION

Model grid

Model cell used in particle tracking 

analyses

Model cell included in the area 

contributing recharge

Water-particle pathline—Shows path 

of groundwater flow. Each line 

color represents different source 

domestic well

320330

450

440

430

Column

Row

Figure 13. Recharge capture areas defined by (A) particle ending locations outlining the area contributing recharge, and (B) 

all area circumscribed by particle pathways outlining the areal extent of the zone of contribution for a group of six groundwater 

sampling locations near Zillah, Washington.
Simulation of Particle Tracking 23

tac15-1019_fig13b

B

Zillah

Toppenish

Yakima

River

082

Base modified from U.S. Geological Survey data and other digital sources. 

Lambert Conformal Conic projection 

Horizontal Datum: North American Datum of 1983 (NAD 83)

0 1 MILE

0 1 KILOMETER

120°18' 16' 14' 120°12'

46°

28'

26'

46°

24'

EXPLANATION

Model grid

Model cell used in particle tracking 

analyses

Model cell included in the areal 

extent of the zone of contribution 

to the well

Water-particle pathline—Shows path 

of groundwater flow. Each line 

color represents different source 

domestic well

320330

450

440

430

Column

Row

Figure 13.—Continued
24 Particle Tracking for Selected Groundwater Wells in the Lower Yakima River Basin, Washington

tac15-1019_fig14a

A

A

Zillah

Granger

Toppenish

Yakima River

082

022

023

320330340350360

470

460

450

440

430

420

Base modified from U.S. Geological Survey data and other digital sources. 

Lambert Conformal Conic projection 

Horizontal Datum: North American Datum of 1983 (NAD 83)

0 1 2 MILES

0 1 2 KILOMETERS

120°16' 14' 12' 10' 08' 120°06'

46°

28'

26'

24'

22'

46°

20'

EXPLANATION

Model grid

Model cell used in particle tracking 

analyses

Model cell included in the areal 

extent of the zone of contribution 

of the well

Water-particle pathline—Shows path 

of groundwater flow. Each line 

color represents different 

source domestic well

Pathline endpoint—Particle tracking 

ending location

Column

Row

Figure 14. Zones of contribution for a dense area of wells near Granger, Washington, divided into two groups (A and B) to make 

individual particle pathlines more visible in a location where areas contributing recharge to specific wells overlap heavily.
Simulation of Particle Tracking 25

tac15-1019_fig14b

Zillah

Granger

Toppenish

Yakima River

082

022

023

Base modified from U.S. Geological Survey data and other digital sources. 

Lambert Conformal Conic projection 

Horizontal Datum: North American Datum of 1983 (NAD 83)

0 1 2 MILES

0 1 2 KILOMETERS

120°16' 14' 12' 10' 08' 120°06'

46°

28'

26'

24'

22'

46°

20'

EXPLANATION

Model grid

Model cell used in particle tracking 

analyses

Model cell included in the areal 

extent of the zone of contribution 

of the well

Water-particle pathline—Shows path 

of groundwater flow. Each line 

color represents different 

source domestic well

Pathline endpoint—Particle tracking 

ending location

B

320330340350360

470

460

450

440

430

420

Column

Row

Figure 14.—Continued
26 Particle Tracking for Selected Groundwater Wells in the Lower Yakima River Basin, Washington

Solute transport in groundwater systems is governed 

by advection (simulated by MODPATH) and hydrodynamic 

dispersion (mechanical dispersion and diffusion; not simulated 

by MODPATH). When a contaminant reaches the water table, 

the contaminant is advected in the direction of and at the rate 

of groundwater flow, but the contaminant also mechanically 

disperses and diffuses longitudinally and transversely to the 

predominant flow direction. Longitudinal hydrodynamic 

dispersion tends to smear the contaminant distribution along 

the longitudinal axis of the flow path and tends to increase 

the apparent advective flow velocity based on an earlier first 

arrival time of the contaminant at a given distance. Transverse 

hydrodynamic dispersion spreads the contaminant plume in 

the horizontal and vertical directions, roughly perpendicular to 

the flow direction. Taken together, hydrodynamic dispersion 

tends to reduce (attenuate) the maximum concentration of 

contaminant at any one point, but increase the volume of 

contaminated aquifer material as compared to advective 

transport. In the context of the particle pathlines presented 

here, the notable implication of simulations that do not 

consider hydrodynamic dispersion is that contaminant 

sources outside of the simulated zones of contribution also 

may contribute to contaminant concentrations measured 

at groundwater sampling locations. The degree to which 

the extent of any zone of contribution is underestimated 

varies principally with the length of the particle pathlines 

and variable groundwater flow velocities imposed by 

heterogeneities of the aquifer material; shorter particle 

pathlines and relatively less heterogeneous aquifer material 

are more likely to accurately reflect true zones of contribution, 

whereas actual zones of contribution associated with longer 

pathlines are more likely to be larger than simulated. A more 

complex analysis, such as simulations using the MODFLOWbased contaminant transport code MT3D (Zheng and others, 

2001), would be needed to examine the effects of advective 

and dispersive transport and thus more accurately identify 

contaminant capture areas for specific model cells containing 

groundwater wells. 

Figure 15 shows the areal extent of the combined zone 

of contribution simulated for all 89 model cells containing 

wells with reported nitrate concentrations exceeding 10 mg/L. 

The combined zone of contribution should not be considered 

the only potential sources of nitrate contamination in the 

basin; it represents only the zone that contributed water to 

these 89 model cells. Contaminant sources outside these areas 

may contribute to measured concentrations by dispersive 

transport. The process for selecting wells for nitrate sampling 

was non-random, so additional areas of high nitrate loading 

may exist in areas not included in the zone of contribution 

determined by the model. A more accurate analysis of 

contaminant presence and transport throughout the study area 

would include specific concentration loading maps for all 

sources and would simulate nitrate concentrations for each 

model cell.

Model Uncertainties and Limitations

The groundwater model represents an extremely 

complex natural system with a set of mathematical equations 

that describe the relevant physical processes governing 

groundwater flow. Intrinsic to the model is the error and 

uncertainty associated with the approximations, assumptions, 

and simplifications inherent in this approach. In addition to 

those intrinsic errors, hydrologic modeling errors may be the 

consequence of a combination of errors in the (1) input data, 

(2) representation of the physical processes by the algorithms 

and discretization of the numerical model, and (3) parameter 

estimation during the calibration procedure (Troutman, 1985).

Input data errors include incorrect or poorly known 

mapped extents or thicknesses of geologic layers. Most 

measurements used to generate model inputs were done in 

populated areas near the Yakima River, but even in areas with 

dense coverage of well logs, variations in lithology or basic 

sediment properties are sure to exist between the locations of 

reported values. The recharge estimates used were derived 

from physical process models that preserve water balance, 

but uncertainties about the inputs to those models carry over 

into uncertainties in this model where groundwater-flow 

velocities, and thus particle time-of-travel, are directly linked 

to recharge fluxes.

Numerical algorithm errors in a simulation are in part 

the result of shortcomings in the way the model represents 

the governing physical processes, and the result of numerical 

errors in the solutions of the governing partial differential 

equations. The model inevitably relies on simplifying 

assumptions and generalizations of the physical processes 

that complexly affect the results of the simulation. Numerical 

errors are minimized in the model through use of careful 

mass-balance and numerical-solution closure criteria. 

Model discretization errors related to the averaging of 

surface elevations over a 1,000-ft2 model grid cell affect the 

distribution of aquifer properties and the effective depth of 

simulated particle starting locations, which are referenced 

to estimated land surface elevations. For this reason, 

interpretations of simulation results should be limited to scales 

several times greater than the model spatial and temporal 

resolutions of 1,000 ft and 1 month. The semi-analytic nature 

of the MODPATH particle-tracking results suggests a level 

of precision that is does not reflect the spatial and temporal 

averaging inherent in the model construction.

The model calibration process can result in parameterestimation errors of several different types. Inappropriate 

weighting of observation networks may result in calibrated 

values outside the accepted range of natural variability if 

observation types are given unequal cumulative weights. 

Close agreement between simulated and measured values does 

not guarantee that the estimated parameter values represent a 

unique parameter set, and alternate parameter combinations 

may generate model calibration results that are equally valid, 

as assessed by the composite-scaled sensitivities, but that 
Simulation of Particle Tracking 27

tac15-1019_fig15

200

300

400

310

390

380

370

360

350

340

330

320

210

290

280

270

260

250

240

230

220

190

180

170

160

410

460

450

440

430

420

400

410420430440450460470480490

500

510520530540550560570580590

Sunnyside

Prosser

Grandview

Zillah

Toppenish

Granger

Wapato

Union Gap

Mabton

Harrah

097

082

120°00' 20' 120°40'

46°

20'

46°

00'

0 4 6 8 2 MILES

0 4 8

10

10 6 2 KILOMETERS

Yakima River

024

241

023

022

Base modified from U.S. Geological Survey data and other digital sources. Lambert Conformal Conic projection 

Horizontal datum: North American Datum of 1983 (NAD 83)

EXPLANATION

Model grid (10 × 10 cells)Model cell included in the areal 

extent of the zone of contribution of the well

Column

Row

Figure 15. Simulated zones of contribution for the model cells containing the 121 contaminated wells in the lower Yakima River Basin, Washington
28 Particle Tracking for Selected Groundwater Wells in the Lower Yakima River Basin, Washington

result in locally different values for important parameters such 

as hydraulic conductivity that would alter particle velocities. 

Limitations of the observation weighting scheme used to 

calibrate this model include equal weighting of groundwater 

head measurements without accounting for relative 

measurement error, uncertainty in the model layer assignments 

of some pumping wells, errors in heads used as calibration 

targets(for example, heads associated with possible perched 

zones above the saturated, regional groundwater-flow system), 

and streamflow weights not based on streamflow variance.

These limitations suggest that using this model to 

simulate the locations of particles derived from specific 

wells is an analysis best done in aggregate. Location errors 

on the average of at least one grid cell (1,000 ft) should be 

expected, with larger errors anticipated for longer flow paths. 

Average values of particle time-of-travel and direction for 

each group of 27 particles representing one starting location 

are representative of all wells contained within that model cell, 

but distinguishing attribution of surface recharge from specific 

geographic locations within that area to individual wells 

contained within that model cell would require a more detailed 

analysis than is possible at this resolution. Domestic wells 

were not individually simulated because of the large grid size, 

and a much finer grid would be needed to more accurately 

simulate pumping effects at these wells. Including the effects 

of hydrodynamic dispersion, not simulated by MODPATH, 

would expand in all directions the effective area from which 

any given well might capture surface-applied contaminant 

loads beyond the zones of contribution reported here as 

sources of water to those wells. Simulation of nitrate transport 

through the aquifer would require a more detailed approach 

that accounts for these complexities.

Summary

Concerns about the health risks associated with high 

concentrations of nitrate measured in residential drinkingwater wells in the lower Yakima River Basin from 2008 to 

2011 have focused renewed attention on land management 

practices that may be contributing excess nitrate to the 

shallow groundwater system and drinking-water aquifers. 

Inorganic fertilizer applied to areas of irrigated agriculture 

and manure infiltration associated with livestock operations 

are likely the two largest sources of nitrogen to the system, 

but their contribution can vary depending on application rate 

and location. No previous work has attempted to identify the 

parts of the basin where nutrient application from any source 

might result in increased nitrate concentrations measured at 

specific wells because knowledge of the groundwater-flow 

system required to do this kind of analysis was previously 

unavailable. Particle tracking analyses can potentially 

identify where specific nitrate loading rates may result in 

contamination of drinking-water wells, and at what locations 

increased nitrate loading rates to groundwater might be 

possible without affecting drinking-water wells.

Groundwater flow travel-times in the lower Yakima 

River Basin are shorter in the shallow basin-fill sediments. 

These sediments consist primarily of alluvial fan, loess, dune 

sand, alluvial terrace, continental sedimentary, and Touchet 

Formation and Ellensburg Formation deposits, and typically 

are several hundred feet thick in the area of interest. The basin 

sediments are contained by anticlinal ridges of the upper 

members of the Columbia River Basalt Group, which also may 

contribute significant groundwater flow to deeper wells. In 

the absence of strong local pumping centers or high irrigation 

rates, groundwater flow within the sediments typically is away 

from the basin margins and towards the Yakima River, and 

in the direction of downstream river flow. Sudden changes in 

river stage, such as from reservoir releases in the upper basin, 

may temporarily disrupt natural groundwater head gradients 

and alter groundwater-flow directions or velocities.

Groundwater flow in the Yakima River Basin was 

simulated using the transient three-dimensional finite 

difference model MODFLOW-2000. Unique flow vectors 

for every active model cell over a 42-year simulation period 

were used with the particle tracking post-processing package 

MODPATH version 5.0 to trace the simulated path of particles 

that arrived at specified model cells containing groundwater 

wells. Using MODPATH, these imaginary particles were 

backtracked from the model cells that contain the wells of 

interest to the beginning of the model simulation period or to 

the water table. For each model cell containing one or more 

wells of interest, the distribution of particle locations at the 

water table was used to define the area contributing recharge 

to the cell. Particle starting locations were selected based 

on the locations and estimated depths of wells with nitrate 

concentrations exceeding the U.S. Environmental Protection 

Agency drinking water standard of 10 milligrams per liter.

Particle tracks originating from shallow aquifer locations 

generally had shorter flowpaths than those from deeper 

locations. The average time-of-travel for all particles was 

3,749 days, but travel times were shorter than the average for 

75.4 percent of the particles because a minority of particles 

originating from deep aquifer locations had very long 

simulated particle flowpaths. For particles that reached the 

water table during the 42-year simulation period, the median 

travel time was 984 days and the median horizontal distance 

between starting and ending particle locations was 13,194 ft. 

The model predicts that any changes to nutrient loading rates 

at these recharge locations should manifest as changes in 

measured nutrient concentrations in 984 days (median) or 

3,749 days (average) at the 89 model cells containing the wells 

of interest, with significant outliers for pathlines that have 

much longer and much shorter travel times. Individual particle 

pathlines ranged in duration from 22 days to the entire 42-year 

simulation period.

All particle pathlines and times-of-travel were calculated 

based on the results of a MODFLOW-2000 model of the 

Yakima River Basin that makes a series of assumptions 

about the complex natural groundwater-flow system. The 

1,000 ft grid cell size assumes homogeneous properties within 
References Cited 29

each model cell, and the seasonal variations in recharge 

and pumping were reduced to monthly mean values held 

constant for the duration of each monthly model stress 

period. Some geologic structures such as anticlines and faults 

were simulated in the model as flow barriers even though 

the location and continuity of these structures at depth is 

unknown. Hydraulic properties used in the model such as 

hydraulic conductivity were based on previously published 

ranges for the relevant geologic units, and calibrated using 

groundwater head and streamflow measurements that did 

not account for variable observation weighting based on 

measurement errors. Some model properties, such as porosity, 

were given default values in the absence of any reliable 

measurements. The 48 identified geologic units in the model 

domain were simulated using 24 computational layers and 

an averaged model cell property weighting scheme that 

effectively blurs the hydraulic property contrasts that are 

assumed to exist at well-defined natural geologic contacts. 

MODPATH uses a semi-analytic particle-tracking algorithm 

based on advective transport and that neglects effects of 

hydrodynamic dispersion, suggesting that the simulated 

zones of contribution are likely to be more accurate for the 

water contributed to a cell of interest than the contaminants 

measured in that water. A more complex analysis, such as an 

analysis using the MODFLOW-based contaminant transport 

code MT3D, could be combined with specific concentration 

loading maps of various potential sources to simulate 

advective and hydrodynamic-dispersive contaminant transport 

and thus to more accurately identify contaminant capture 

areas for specific model cells containing groundwater wells 

as part of a future study addressing these issues. The results 

presented in this report should be interpreted in the context of 

the recognized effect of model errors on groundwater flow and 

the derived particle-tracking results, and the known limitations 

inherent in the particle-tracking approach.

References Cited

Alley, W.M., 2003, Terminology used in studies of the 

source of water to wells under steady-state conditions: 

U.S. Geological Survey Office of Ground Water Technical 

Memorandum 2003.02, accessed July 22, 2015, at https://

water.usgs.gov/admin/memo/GW/gw03.02.html.

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 

00-342, 89 p.

Bauer, H.H., and Vaccaro, J.J., 1987, Documentation of a deep 

percolation model for estimating ground-water recharge: 

U.S. Geological Survey Open-File Report 86-536, 180 p.

Doherty, J.E., 2010, PEST, Model-independent parameter 

estimation—User manual (5th ed., with additions): 

Brisbane, Australia, Watermark Numerical Computing.

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.

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.

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, 81 p., 15 pls.

Harbaugh, A.W., Banta, E.R., Hill, M.C., and McDonald, 

M.G., 2000, MODFLOW-2000, the U.S. Geological 

Survey modular ground-water model—User guide to 

modularization concepts and the ground-water flow process: 

U.S. Geological Survey Open-File Report 2000-92, 121 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.

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., 7 pls.

Leavesley, G.H., Lichty, R.W., Troutman, B.M., and Saindon, 

L.G., 1983, Precipitation-runoff modeling system—

User’s manual: U.S. Geological Survey Water-Resources 

Investigations Report 83-4238, 207 p.

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.

Niswonger, R.G., and Prudic, D.E., 2005, Documentation 

of the Streamflow-Routing (SFR2) Package to include 

unsaturated flow beneath streams—A modification to SFR1: 

U.S. Geological Survey Techniques and Methods, book 6, 

chap. A13, 47 p.
30 Particle Tracking for Selected Groundwater Wells in the Lower Yakima River Basin, Washington

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.

Pollock, D.W., 1994, User’s guide for MODPATH/

MODPATH-PLOT, Version 3—A particle tracking postprocessing package for MODFLOW, the U.S. Geological 

Survey finite-difference ground-water flow model: U.S. 

Geological Survey Open-File Report 94-464, 236 p.

Pollock, D.W., 2012, User guide for MODPATH, Version 

6—A particle-tracking model for MODFLOW: U.S. 

Geological Survey Techniques and Methods, book 6, 

chap. A41, 58 p.

Prudic, D.E., Konikow, L.F., and Banta, E.R., 2004, A new 

streamflow-routing (SFR1) package to simulate streamaquifer interaction with MODFLOW-2000: U.S. Geological 

Survey Open-File Report 2004-1042, 96 p.

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, Washington, Pacific Northwest National 

Laboratory [variously paged].

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. Environmental Protection Agency, 2012, Relation 

between nitrate in water wells and potential sources in the 

lower Yakima Valley, Washington: U. S. Environmental 

Protection Agency, EPA-910-R-12-003.

U.S. Geological Survey, 2014, Hydrography—The National 

Hydrography Dataset and the Watershed Boundary Dataset: 

Web site, accessed July 20, 2015, at http://nhd.usgs.gov/

index.html. 

Vaccaro, J.J., 2007, A deep percolation model for estimating 

ground-water recharge—Documentation of modules for the 

modular modeling system of the U.S. Geological Survey: 

U.S. Geological Survey Scientific Investigations Report 

2006-5318, 30 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.

Vaccaro, J.J., and Olsen, T.D., 2007, Estimates of groundwater recharge to the Yakima River Basin aquifer system, 

Washington, for predevelopment and current land-use and 

land-cover conditions: U.S. Geological Survey Scientific 

Investigations Report 2007–5007, 30 p.

Vaccaro, J.J., and Sumioka, S.S., 2006, Estimates of groundwater pumpage from the Yakima River Basin aquifer 

system, Washington, 1960–2000: U.S. Geological Survey 

Scientific Investigations Report 2006–5205, 56 p.

Van Winkle, Walton, 1914, Quality of the surface waters of 

Washington: U.S. Geological Survey Water-Supply Paper 

339, 105 p.

Washington State Department of Ecology, 2010, Lower 

Yakima Valley groundwater quality preliminary assessment 

and recommendations document: Washington State 

Department of Ecology Publication No. 10-10-009.

Washington State Department of Ecology, 2015, Well logs—

Washington State well log viewer: Web site, accessed July 

16, 2015, at https://fortress.wa.gov/ecy/waterresources/map/

WCLSWebMap/default.aspx.

Zheng, C., Hill, M.C., and Hsieh, P.A., 2001, 

MODFLOW-2000, the U.S. Geological Survey modular 

ground-water model—User guide to the LMT6 package, 

the linkage with MT3DMS for multi-species mass transport 

modeling: U.S. Geological Survey Open-File Report 01-82, 

43 p.
Appendix A 31

Appendix A. Summarized Particle-Tracking Information

Well depths for each sample location where the screened 

interval depth of the well was not available from a well log 

were estimated based on the recorded average depth of all 

similar wells within a 1-mi radius of the sampling location. 

Model locations and estimated depths (table A1) were used 

to determine the model layer to which wells and thus particle 

starting locations were assigned. Particles were distributed 

throughout each model cell identified by these depths, which 

results in particles starting from a range of depths equal to the 

thickness of the assigned layer at that model location prior to 

initializing the backtracking process. 

Nitrate concentrations were greater than 10 mg/L in 

121 domestic wells; these wells were located within 89 

unique model cells. Particle pathlines for each of these 89 

model cells simulate the historical travel of groundwater in 

those cells at the end of the simulation period, and are shown 

in appendix table A1in approximate order form north to 

south. The summarized particle-tracking values represent the 

average value of all particles that were assigned to each of the 

89 starting cell locations, but each of those individual particles 

had a unique pathline and time-of-travel not reported here. 

The coordinate starting locations labeled “Starting 

X” and “Starting Y” are the locations of the center of the 

associated model cell, listed in Washington State Plane South, 

North American Datum of 1983. In this context, starting 

locations do not represent the actual starting locations on the 

cell faces but represent the model cells containing wells from 

which particle backtracking was calculated, and final locations 

are the places from which those particles originated before 

flowing towards the wells.

The averaged endpoint distances are the straight line 

distance between the particle starting and ending locations, 

and do not represent cumulative pathline length for particles 

that travel in varying directions or circuitous routes. The 

average final location model layer of backtracked particles is 

reported as 1.0 for any cell from which all tracked particles 

reached land surface within the 42-year simulation period. 

Values greater than 1.0 for any cell’s average final layer of 

particles indicate that at least some of those particles did not 

reach land surface. Similarly, an average pathline duration of 

15,341 days (42 years) suggests that all particles from that 

cell were still being actively tracked in the model domain 

at the end of the 42-year simulation period. In some cases, 

groundwater sampling locations with associated particle 

pathlines that ended the backtracking process in a deeper 

model layer than the one from which they started suggest 

upward groundwater flow at those locations. Such wells 

typically are acquiring older groundwater from deeper aquifer 

layers with longer flow paths, which may have associated 

recharge capture areas that are relatively far away and thus 

indicate distant nitrate loads. Particle pathlines that parallel 

a model no-flow boundary, such as those along the Horse 

Heaven Hills, may not indicate true groundwater-flow 

directions, which may include horizontal groundwater flow 

beneath or through the bounding structure. 
32 Particle Tracking for Selected Groundwater Wells in the Lower Yakima River Basin, Washington

Table A1. Model locations and estimated depth in feet used to determine the model layer to which wells, and thus particle starting 

locations, were assigned.

[Coordinate start location: “Starting X” and “Starting Y” are the locations of the center of the associated model cell, listed in Washington State Plane 

South, North American Datum of 1983]

Well 

location ID

Assigned 

depth 

(feet)

Coordinate start location Model starting location Average 

final 

location 

(layer)

Average 

endpoint 

distances 

(feet)

Average 

pathline 

duration 

(days)

Starting X 

(feet)

Starting Y 

(feet) Column Row Layer

1 100 245500 176500 246 424 2 1.0 17,256 855.6

2 280 245500 166500 246 434 5 1.0 37,638 5,276.7

3 153 320500 165500 321 435 3 1.0 7,938 1,399.4

4 152 320500 164500 321 436 3 1.0 8,204 1,516.4

5 147 323500 159500 324 441 3 1.0 10,071 705.9

6 140 326500 153500 327 447 3 1.0 11,191 564.0

7 143 326500 152500 327 448 3 1.0 13,684 515.4

8 120 352500 151500 353 449 3 1.0 6,197 1,137.9

9 105 352500 150500 353 450 3 1.0 6,655 1,401.0

10 120 354500 150500 355 450 3 1.0 7,161 1,452.2

11 105 355500 150500 356 450 3 1.0 9,209 2,254.6

12 185 373500 150500 374 450 4 2.6 18,555 10,868.9

13 123 329500 148500 330 452 3 1.0 18,006 2,891.8

14 145 349500 148500 350 452 3 1.0 7,536 853.9

15 185 345500 147500 346 453 4 1.0 6,754 1,075.9

16 154 347500 147500 348 453 3 1.0 7,211 695.0

17 157 348500 147500 349 453 3 1.0 7,308 706.8

18 145 349500 147500 350 453 3 1.0 7,902 743.3

19 115 354500 147500 355 453 3 1.0 8,585 870.2

20 125 356500 147500 357 453 3 1.0 12,171 1,486.6

21 105 386500 147500 387 453 3 1.3 41,040 6,411.3

22 140 348500 145500 349 455 3 1.0 8,371 766.9

23 181 368500 145500 369 455 4 1.0 8,737 3,025.5

24 162 372500 145500 373 455 4 1.0 8,408 1,436.2

25 210 386500 145500 387 455 4 1.8 41,813 8,265.5

26 181 370500 144500 371 456 4 1.0 13,496 3,100.5

27 175 379500 144500 380 456 4 1.3 44,217 9,472.5

28 135 385500 144500 386 456 3 1.1 36,766 4,703.0

29 212 386500 144500 387 456 4 1.1 45,411 6,086.2

30 200 387500 144500 388 456 4 2.0 43,189 8,798.3

31 315 387500 143500 388 457 6 4.8 19,548 15,213.7

32 126 345500 142500 346 458 3 1.0 13,885 1,147.9

33 140 349500 142500 350 458 3 1.0 10,488 1,039.2

34 145 354500 142500 355 458 3 1.0 11,355 627.8

35 145 370500 142500 371 458 3 1.0 8,401 929.6

36 157 373500 142500 374 458 3 2.5 25,648 8,518.1

37 182 381500 142500 382 458 4 1.5 36,521 8,080.9

38 185 381500 141500 382 459 4 1.1 28,949 5,858.0

39 185 387500 141500 388 459 4 2.3 44,271 7,792.5

40 175 382500 140500 383 460 4 1.1 34,367 4,851.5

41 190 387500 140500 388 460 4 1.0 48,707 1,689.1

42 140 351500 139500 352 461 3 1.0 10,784 620.4

43 137 354500 139500 355 461 3 1.3 17,189 3,403.4

44 145 363500 139500 364 461 3 1.0 10,160 669.4

45 115 346500 138500 347 462 3 1.0 11,143 569.3

46 185 385500 137500 386 463 4 1.0 52,845 1,735.8

47 98 344500 135500 345 465 2 1.0 6,036 80.9

48 141 363500 135500 364 465 3 1.0 15,978 907.8

49 140 374500 135500 375 465 3 1.0 14,060 1,517.2
Appendix A 33

Well 

location ID

Assigned 

depth 

(feet)

Coordinate start location Model starting location Average 

final 

location 

(layer)

Average 

endpoint 

distances 

(feet)

Average 

pathline 

duration 

(days)

Starting X 

(feet)

Starting Y 

(feet) Column Row Layer

50 141 365500 134500 366 466 3 1.0 16,688 886.3

51 145 373500 132500 374 468 3 1.0 18,141 435.3

52 145 384500 131500 385 469 3 1.0 22,302 422.9

53 160 384500 130500 385 470 4 1.2 51,483 2,270.0

54 175 384500 130500 385 470 3 1.0 19,644 471.7

55 143 344500 129500 345 471 3 1.8 33,392 7,750.0

56 145 376500 129500 377 471 3 1.0 17,521 475.0

57 163 398500 129500 399 471 4 1.2 33,009 4,496.7

58 142 402500 129500 403 471 3 1.0 20,924 2,907.2

59 132 396500 126500 397 474 3 1.0 20,150 735.5

60 150 401500 126500 402 474 3 1.0 22,679 1,700.2

61 142 403500 126500 404 474 3 1.0 21,108 1,369.0

62 115 406500 126500 407 474 3 1.0 26,205 3,162.1

63 250 412500 123500 413 477 5 1.0 10,883 9,019.2

64 118 387500 121500 388 479 3 1.1 46,079 2,354.1

65 100 393500 112500 394 488 3 1.0 24,055 293.6

66 100 397500 111500 398 489 3 1.0 33,941 851.1

67 110 385500 108500 386 492 3 1.0 37,167 849.6

68 110 385500 107500 386 493 3 1.0 35,269 468.1

69 122 396500 102500 397 498 3 1.0 29,424 2,358.6

70 142 392500 92500 393 508 3 2.4 47,572 7,680.7

71 142 392500 91500 393 509 3 4.6 64,334 14,373.4

72 139 382500 86500 383 514 3 1.0 9,787 759.3

73 110 385500 83500 386 517 3 2.3 37,180 5,026.8

74 130 386500 83500 387 517 3 2.7 19,671 5,304.6

75 110 389500 83500 390 517 3 4.3 17,225 5,594.3

76 100 399500 83500 400 517 3 5.0 52,788 14,720.4

77 110 388500 82500 389 518 3 4.0 9,000 5,470.4

78 90 395500 81500 396 519 2 1.0 23,625 1,874.2

79 139 405500 81500 406 519 3 8.8 18,723 15,341.0

80 130 386500 80500 387 520 3 2.7 18,735 5,378.4

81 90 394500 80500 395 520 2 1.0 9,050 614.6

82 90 395500 80500 396 520 2 1.0 21,185 1,573.5

83 90 396500 80500 397 520 2 1.0 20,771 1,107.2

84 130 386500 79500 387 521 3 2.7 27,998 6,039.7

85 90 394500 79500 395 521 2 1.0 13,317 1,199.0

86 125 368500 76500 369 524 3 4.6 86,411 12,180.2

87 125 385500 76500 386 524 3 3.0 51,891 8,616.3

88 240 378500 75500 379 525 5 9.1 20,601 15,341.0

89 125 385500 75500 386 525 3 3.2 40,115 7,655.5

Table A1. Model locations and estimated depth in feet used to determine the model layer to which wells, and thus particle starting 

locations, were assigned.—Continued

[Coordinate start location: “Starting X” and “Starting Y” are the locations of the center of the associated model cell, listed in Washington State Plane 

South, North American Datum of 1983]

Publishing support provided by the U.S. Geological Survey

Science 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
Bachmann—Particle Tracking for Selected Groundwater Wells in the Lower Yakima River Basin, Washington—Scientific Investigations Report 2015–5149

ISSN 2328-0328 (online)

http://dx.doi.org/10.3133/sir20155149