Cover photographs:

Top: Steelhead Falls on the Deschutes River near Crooked River Ranch, Oregon. 

Middle: Crooked River Canyon at Crooked River Ranch, Oregon. 

Bottom: North and Middle Sister with a wheel-line irrigation system in the foreground, 

near Sisters, Oregon. (Photographs by Rodney R. Caldwell, U.S. Geological Survey.) 
U.S. Department of the Interior

U.S. Geological Survey

Simulation of Regional Ground-Water Flow 

in the Upper Deschutes Basin, Oregon

By MARSHALL W. GANNETT and KENNETH E. LITE JR.

Water-Resources Investigations Report 03-4195 

Prepared in cooperation with Oregon Water Resources Department; 

Cities of Bend, Redmond, and Sisters; Deschutes and Jefferson Counties; 

The Confederated Tribes of the Warm Springs Reservation of Oregon; 

and U.S. Environmental Protection Agency

Portland, Oregon 

2004 
U. S. DEPARTMENT OF THE INTERIOR

GALE A. NORTON, Secretary

U.S. GEOLOGICAL SURVEY

CHARLES G. GROAT, Director

The use of trade, product, or firm names in this publication 

is for descriptive purposes only and does not imply 

endorsement by the U.S. Government.

For additional information: Copies of this report may be purchased from:

District Chief U.S. Geological Survey 

U.S. Geological Survey Branch of Information Services 

10615 S.E. Cherry Blossom Dr. Box 25286 

Portland, OR 97216-3103 Denver Federal Center 

E-mail: info-or@usgs.gov Denver, CO 80225-0286

Internet: http:// or.water.usgs.gov Telephone: 1-888-ASK-USGS

Suggested citation: Gannett, M.W., and Lite, K.E., Jr., 2004, Simulation of regional ground-water flow in the upper 

Deschutes Basin, Oregon: U.S. Geological Survey Water-Resources Investigations Report 03–4195, 84 p.
iii

CONTENTS

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

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

Background and Study Objectives....................................................................................................................................2

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

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

Acknowledgements...........................................................................................................................................................3

Ground-Water Hydrology ...........................................................................................................................................................3

Geologic Framework.........................................................................................................................................................3

Ground-Water Flow ..........................................................................................................................................................8

Hydraulic Characteristics................................................................................................................................................11

Horizontal Hydraulic Conductivity.......................................................................................................................11

Vertical Hydraulic Conductivity...........................................................................................................................13

Storage Coefficients..............................................................................................................................................13

Model Description.....................................................................................................................................................................14

Governing Equations and Model Code...........................................................................................................................14

Discretization ..................................................................................................................................................................14

Spatial Discretization ............................................................................................................................................15

Temporal Discretization........................................................................................................................................15

Boundary Conditions ......................................................................................................................................................15

Geographic Boundaries.........................................................................................................................................17

Hydrologic-Process Boundaries............................................................................................................................18

Model Parameters ...........................................................................................................................................................21

Horizontal Hydraulic Conductivity.......................................................................................................................21

Vertical Hydraulic Conductivity...........................................................................................................................23

Storage Coefficients..............................................................................................................................................23

Model Calibration .....................................................................................................................................................................26

Steady-State Calibration .................................................................................................................................................26

Steady-State Calibration Data...............................................................................................................................26

Steady-State Calibration Procedure ......................................................................................................................27

Steady-State Model Parameters ............................................................................................................................29

Parameter Sensitivity ............................................................................................................................................29

Final Parameter Values .........................................................................................................................................30

Steady-State Calibration Model Fit.................................................................................................................................32

Statistical Measures of Model Fit and Parameter Uncertainty .............................................................................32

Comparison of Simulated and Observed Steady-State Heads ..............................................................................35

Comparison of Simulated and Observed Steady-State Ground-Water Flow to and from Streams ......................44

Steady-State Model Water Budget........................................................................................................................48

Transient Calibration.......................................................................................................................................................48

Transient Calibration Data ....................................................................................................................................48

Transient Calibration Procedure ...........................................................................................................................49

Comparison of Simulated and Observed Transient Heads ...................................................................................50

Madras/Gateway Area.................................................................................................................................52

Lower Desert Area ......................................................................................................................................52

Black Butte/Upper Metolius River Area.....................................................................................................53

Peninsula Area ............................................................................................................................................53

Lower Squaw Creek/Deep Canyon Area ....................................................................................................54

Redmond Area ............................................................................................................................................54

Sisters/Tumalo Area....................................................................................................................................55

Powell Butte Area .......................................................................................................................................55

Lower Tumalo Creek/Bend Area................................................................................................................56
iv

Benham Falls Area......................................................................................................................................56

Central La Pine Subbasin............................................................................................................................57

Southern La Pine Subbasin .........................................................................................................................58

Summary of Simulated Water-Level Fluctuations......................................................................................58

Comparison of Simulated and Observed Transient Ground-Water Discharge to Streams...................................58

Odell Creek .................................................................................................................................................59

Streams above and below Crane Prairie Reservoir.....................................................................................59

Fall River.....................................................................................................................................................61

Metolius River Drainage.............................................................................................................................61

Deschutes River—Bend to Culver..............................................................................................................62

Lower Crooked River..................................................................................................................................63

Summary of Transient Ground-Water Discharge to Streams .....................................................................65

Transient Model Water Budget.............................................................................................................................65

Model Limitations.................................................................................................................................................65

Example Simulations ................................................................................................................................................................67

Simulating the Effects of Pumping .................................................................................................................................68

Simulated Pumping in the Bend Area...................................................................................................................68

Simulated Pumping in the Redmond Area............................................................................................................71

Simulated Pumping between the Deschutes and Crooked Rivers near Osborne Canyon ....................................74

Simulated Effects of Canal Lining..................................................................................................................................76

Summary and Conclusions........................................................................................................................................................79

References Cited .......................................................................................................................................................................83

FIGURES

1.– 5. Maps showing:

1. Location of the upper Deschutes Basin, Oregon, and major geographic and cultural features ..............................4

2. Lines of equal mean annual precipitation in the upper Deschutes Basin and graphs showing 

mean monthly precipitation at selected precipitation stations ................................................................................5

3. Generalized geology and hydrogeologic units of the upper Deschutes Basin........................................................6

4. Generalized lines of equal head, locations of field-inventoried wells, and land ownership in the 

upper Deschutes Basin ............................................................................................................................................9

5. Model grid, boundary conditions, locations of stream cells, and regions where evapotranspiration 

is active .................................................................................................................................................................16

6. – 7. Graphs showing:

6. Mean monthly recharge from precipitation, canal leakage, and on-farm loss for 1978 to 1997 ..........................19

7. Recharge from precipitation, canal leakage, and on-farm loss for each semiannual stress period 

from 1978 to 1997.................................................................................................................................................19

8. Map showing the distribution of irrigation pumpage and public water-supply pumpage on the 

model grid ...................................................................................................................................................................22

9. Graph showing the rates of ground-water pumpage for irrigation and public-supply use for each 

semiannual stress period from 1978 to 1997 ..............................................................................................................23

10. Maps showing hydraulic-conductivity zonation in all model layers ..........................................................................24

11. Map showing stream fluxes used for steady-state calibration ....................................................................................28

12.–15. Graphs showing:

12. Composite scaled sensitivities of steady-state parameters....................................................................................30

13. Weighted residuals as a function of weighted simulated values...........................................................................33

14. Normal probability plot of ordered weighted residuals ........................................................................................33

15. The final values, 95-percent confidence intervals, and ranges of expected values for optimized 

steady-state hydraulic-conductivity parameters....................................................................................................34
v

16. Maps showing lines of equal measured head and lines of equal simulated head for all model layers .......................36

17.–19. Graphs showing:

17. Measured or estimated steady-state stream fluxes, 95-percent confidence intervals on measured 

or estimated stream fluxes, and simulated steady-state stream fluxes..................................................................46

18. Simulated stream fluxes as a function of measured or estimated stream fluxes...................................................47

19. Measured or estimated steady-state stream fluxes, 95-percent confidence intervals on measured 

or estimated stream fluxes, and simulated steady-state stream fluxes for aggregated stream reaches.................47

20. Map showing subareas used for discussion of water-level fluctuations .....................................................................51

21.–41. Graphs showing: 

21. Simulated and observed water-level fluctuations in the Madras/Gateway area ...................................................52

22. Simulated and observed water-level fluctuations in the Black Butte/upper Metolius River area ........................53

23. Simulated and observed water-level fluctuations in the Peninsula area ...............................................................53

24. Simulated and observed water-level fluctuations in the lower Squaw Creek/Deep Canyon area ........................54

25. Simulated and observed water-level fluctuations in the Redmond area ...............................................................54

26. Simulated and observed water-level fluctuations in the Sisters/Tumalo area ......................................................55

27. Simulated and observed water-level fluctuations in the lower Tumalo Creek/Bend area ....................................56

28. Simulated and observed water-level fluctuations in the Benham Falls area.........................................................56

29. Simulated and observed water-level fluctuations in the central La Pine subbasin ...............................................57

30. Simulated and observed water-level fluctuations in the southern La Pine subbasin ............................................58

31. Daily and monthly mean flows of Odell Creek, and simulated ground-water discharge, 

1977 to 1997..........................................................................................................................................................59

32. Monthly mean flow of Cultus Creek, and simulated ground-water discharge, 1977 to 1997 ..............................60

33. Daily and monthly mean flows of the Deschutes River above Snow Creek, and simulated 

ground-water discharge, 1977 to 1997..................................................................................................................60

34. Simulated ground-water discharge to the Deschutes River between Crane Prairie Reservoir 

and Davis Creek, and estimated long-term ground-water discharge, 1977 to 1997 .............................................60

35. Monthly mean flow of Fall River and simulated ground-water discharge, 1977 to 1997 ....................................61

36. Monthly mean and October mean flow of the Metolius River, October mean flow of the 

Metolius River minus the October mean flows of Whitewater River and Jefferson Creek, 

and simulated ground-water discharge, 1977 to 1997 ..........................................................................................62

37. Measured monthly mean inflow to the Deschutes River between Bend and Culver, 1977 to 1997, 

and simulated inflow for all time steps .................................................................................................................63

38.Monthly mean flows of the Crooked River below Opal Springs and at Terrebonne, and the 

difference between the two gages, 1967 to 1973 ..................................................................................................64

39. July mean flows of the Crooked River below Opal Springs and at Terrebonne, and the difference 

between the two gages, 1968 to 1973 ...................................................................................................................64

40. July mean flows of the Crooked River below Opal Springs and simulated ground-water discharge 

to the lower Crooked River, 1977 to 1997............................................................................................................64

41. Variations in the transient model water budget by stress period ..........................................................................67

42. Maps showing simulated drawdown in model layer 2 and net reduction in ground-water discharge 

to streams after 2 and 10 years of pumping 10 ft3/s from model layer 2 in the Bend area.........................................69

43. Graph showing variations with time in the amounts of ground water coming from storage and from 

diminished streamflow in response to pumping 10 ft3/s from model layer 2 in the Bend area..................................70

44. Maps showing simulated drawdown in model layer 7 and net reduction in ground-water discharge 

to streams after 2 and 10 years of pumping 10 ft3/s from model layer 7 in the Bend area.........................................70 

45. Graph showing variations with time in the amounts of ground water coming from storage and from 

diminished streamflow in response to pumping 10 ft3/s from model layer 7 in the Bend area..................................71
vi

46. Maps showing simulated drawdown and net reduction in ground-water discharge to streams after 

2 and 10 years of pumping 10 ft3/s from model layer 2 in the Redmond area ...........................................................72

47. Maps showing simulated drawdown and net reduction in ground-water discharge to streams after 

2 and 10 years of pumping 10 ft3/s from model layer 7 in the Redmond area ...........................................................73 

48. Graphs showing variations with time in the amounts of ground water coming from storage and from 

diminished streamflow in response to pumping 10 ft3/s from model layers 2 and 7 in the Redmond area ...............73

49. Maps showing simulated drawdown in model layer 2 and net reduction in ground-water discharge 

to streams after 2 and 10 years of pumping 10 ft3/s from model layer 2 between the Deschutes 

and Crooked Rivers near Osborne Canyon.................................................................................................................74

50. Graph showing variations with time in the amounts of ground water coming from storage and from 

diminished streamflow in response to pumping 10 ft3/s from model layer 2 between the Deschutes 

and Crooked Rivers near Osborne Canyon.................................................................................................................75

51. Maps showing simulated drawdown in model layer 7 and net reduction in ground-water discharge 

to streams after 2 and 10 years of pumping 10 ft3/s from model layer 7 between the Deschutes 

and Crooked Rivers near Osborne Canyon.................................................................................................................75

52. Graph showing variations with time in the amounts of ground water coming from storage and from 

diminished streamflow in response to pumping 10 ft3/s from model layer 7 between the Deschutes 

and Crooked Rivers near Osborne Canyon.................................................................................................................76

53. Maps showing simulated head decline in model layer 1 and diminished streamflow due to a 50 percent 

reduction in leakage of main irrigation canals in the vicinity of Bend after 2 and 10 years ......................................77

54. Graph showing variations with time in the amounts of ground water coming from storage and from 

diminished streamflow in response to a 50 percent reduction in leakage of main irrigation canals in 

the vicinity of Bend.....................................................................................................................................................77

55. Maps showing simulated head decline in model layer 1 and diminished streamflow due to a 50 percent 

reduction in leakage of main irrigation canals north of Redmond after 2 and 10 years .............................................78

56. Graph showing variations with time in the amounts of ground water coming from storage and from 

diminished streamflow in response to a 50 percent reduction in leakage of main irrigation canals 

north of Redmond .......................................................................................................................................................79

TABLES

1. Stream groupings used for streambed conductance parameters .................................................................................29

2. Final hydraulic-conductivity values............................................................................................................................31

3. Measured or estimated steady-state stream gains and losses due to ground-water exchange, and 

simulated equivalents..................................................................................................................................................45

4. Measured or estimated gains and losses to aggregated stream reaches, and simulated equivalents...........................47

5. Steady-state model simulated water budget................................................................................................................48

6. Final storage coefficient values...................................................................................................................................49

7. Transient model water budget for each stress period..................................................................................................66
vii

CONVERSION FACTORS AND VERTICAL DATUM 

Elevation: In this report, elevation is referenced to the National Geodetic Vertical Datum of 1929 

(NGVD29). 

Specific conductance is given in microsiemens per centimeter (μS/cm) at 25°C.

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

°C = (°F − 32) / 1.8

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

°F = (1.8 × °C) + 32

MAPPING SOURCES: 

Base map modified from U.S. Geological Survey 1:500,000 State base map, 

1982, with digital data from U.S. Bureau of the Census, TIGER/Line (R), 1990, 

and U.S. Geological Survey Digital Line Graphs published at 1:100,000. 

Projection: Universal Transverse Mercator projection, Zone 10, 

1927 North American Datum.

Multiply By To obtain

inch (in) 25.40 millimeter (mm)

inches per year (in/yr) 25.40 millimeters per year (mm/yr)

foot (ft) 0.3048 meter (m)

feet per second (ft/s) 0.3048 meters per second (m/s)

feet per day (ft/d) 3.528 × 10-6 meters per second (m/s)

feet per year (ft/yr) 0.3048 meters per year (m/yr)

square feet (ft2) 0.0929 square meters (m2)

acre-feet per year (acre-ft/yr) 1,233.482 cubic meters per year (m3/yr)

mile (mi) 1.609 kilometer (km)

square mile (mi2) 2.590 square kilometer (km2)

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

gallon per minute (gal/min) 6.308 × 10-5 cubic meters per second (m3/s)
viii

LOCATION SYSTEM

The system used for locating wells, springs, and surface-water sites in this report is based on the 

rectangular system for subdivision of public land. The State of Oregon is divided into townships of 

36 square miles numbered according to their location relative to the east-west Willamette baseline and 

a north-south Willamette meridian. The position of a township is given by its north-south “Township” 

position relative to the baseline and its east-west “Range” position relative to the meridian. Each 

township is divided into 36 one-square-mile (640-acre) sections numbered from 1 to 36. For example, 

a well designated as 18S/11E-29AAC is located in Township 18 south, Range 11 east, section 29. 

The letters following the section number correspond to the location within the section; the first letter 

(A) identifies the quarter section (160 acres); the second letter (A) identifies the quarter-quarter section 

(40 acres); and the third letter (C) identifies the quarter-quarter-quarter section (10 acres). Therefore, 

well 29AAC is located in the SW quarter of the NE quarter of the NE quarter of section 29. When more 

than one designated well occurs in the quarter-quarter-quarter section, a serial number is appended.

Well- and spring-location system.

Each well is assigned a unique 8-digit identification number known as the log-id number. The first 

two digits of the log-id number indicate the county code from the Federal Information Processing 

Standards (FIPS) code file for the county in which the well exists. The FIPS codes for the counties 

in the study area are as follows: 13, Crook County; 17, Deschutes County; 31, Jefferson County; 

and 35, Klamath County. The last 6 digits of the number correspond to the State of Oregon well-log 

number (a unique number assigned by the Oregon Water Resources Department to the report filed 

by the well driller).
1

Simulation of Regional Ground-Water Flow 

in the Upper Deschutes Basin, Oregon

By Marshall W. Gannett and Kenneth E. Lite Jr. 

Abstract

This report describes a numerical model that 

simulates regional ground-water flow in the upper 

Deschutes Basin of central Oregon. Ground water 

and surface water are intimately connected in the 

upper Deschutes Basin and most of the flow of 

the Deschutes River is supplied by ground water. 

Because of this connection, ground-water pumping and reduction of artificial recharge by lining 

leaking irrigation canals can reduce the amount 

of ground water discharging to streams and, 

consequently, streamflow. The model described in 

this report is intended to help water-management 

agencies and the public evaluate how the regional 

ground-water system and streamflow will respond 

to ground-water pumping, canal lining, drought, 

and other stresses.

Ground-water flow is simulated in the model 

by the finite-difference method using MODFLOW 

and MODFLOWP. The finite-difference grid 

consists of 8 layers, 127 rows, and 87 columns. All 

major streams and most principal tributaries in 

the upper Deschutes Basin are included. Groundwater recharge from precipitation was estimated 

using a daily water-balance approach. Artificial 

recharge from leaking irrigation canals and onfarm losses was estimated from diversion and 

delivery records, seepage studies, and crop data. 

Ground-water pumpage for irrigation and public 

water supplies, and evapotranspiration are also 

included in the model. 

The model was calibrated to mean annual 

(1993–95) steady-state conditions using parameter-estimation techniques employing nonlinear 

regression. Fourteen hydraulic-conductivity 

parameters and two vertical conductance parameters were determined using nonlinear regression. 

Final parameter values are all within expected 

ranges. The general shape and slope of the simulated water-table surface and overall hydraulichead distribution match the geometry determined 

from field measurements. The fitted standard 

deviation for hydraulic head is about 76 feet. The 

general magnitude and distribution of groundwater discharge to streams is also well simulated 

throughout the model. Ground-water discharge to 

streams in the area of the confluence of the Deschutes, Crooked, and Metolius Rivers is closely 

matched. 

The model was also calibrated to transient 

conditions from 1978 to 1997 using traditional 

trial-and-error methods. Climatic cycles during this 

period provided an excellent regional hydrologic 

signal for calibration. Climate-driven water-level 

fluctuations are simulated with reasonable accuracy over most of the model area. The timing and 

magnitude of simulated water-level fluctuations 

caused by annual pulses of recharge from precipitation match those observed reasonably well, 

given the limitations of the time discretization in 

the model. Water-level fluctuations caused by 

annual canal leakage are simulated very well over 

most of the area where such fluctuations occur. 

The transient model also simulates the volumetric 

distribution and temporal variations in groundwater discharge reasonably well. The match 

between simulated and measured volume of and 

variations in ground-water discharge is, however, 

somewhat dependent on geographic scale. 

The rates of and variations in ground-water 

discharge are matched best at regional scales. 

Example simulations were made to demonstrate the utility of the model for evaluating the 

effects of ground-water pumping or canal lining. 

Pumping simulations show that pumped water 

comes largely from aquifer storage when pumping 

begins, but as the water table stabilizes, the 

pumping increasingly diminishes the discharge 

to streams and, hence, streamflow. The time it 

takes for pumping to affect streamflow varies 

spatially depending, in general, on the location 

of pumping relative to the discharge areas. Canallining simulations show similar effects. 
2

INTRODUCTION

Background and Study Objectives

The upper Deschutes Basin is presently one of 

the fastest growing population centers in Oregon. The 

number of people in Deschutes County, the most 

populous county in the basin, quadrupled between 

1970 and 2001. Approximately 160,000 people lived 

in the upper Deschutes Basin as of 2001. Growth in 

the region is expected to continue, and residents and 

government agencies are concerned about water supplies for the growing population and the consequences 

of increased development for existing water users. 

Surface-water resources in the area have been closed 

to additional appropriation for many years. Therefore, 

virtually all new development in the region must rely 

on ground water as a source of water. Prior to this 

study, very little quantitative information was available on the ground-water hydrology of the basin. 

This lack of information made ground-water resourcemanagement decisions difficult and was generally a 

cause for concern.

To fill this information void, the U.S. Geological 

Survey (USGS) conducted a cooperative study with 

the Oregon Water Resources Department (OWRD); 

the Cities of Bend, Redmond, and Sisters; Deschutes 

and Jefferson Counties; The Confederated Tribes 

of the Warm Springs Reservation of Oregon; and 

the U.S. Environmental Protection Agency. The 

objectives of this study were to provide a quantitative 

assessment of the regional ground-water system and 

provide the understanding and analytical tools for 

State and local government agencies, hydrologists, 

and local residents to make informed resourcemanagement decisions. This report is one in a series 

that presents the results of the upper Deschutes Basin 

ground-water study. The conceptual understanding of 

regional ground-water flow on which the numerical 

model presented herein is based was described by 

Gannett and others (2001). The geologic framework of 

the regional ground-water flow system was described 

by Lite and Gannett (2002). Much of the hydrologic 

data on which the conceptual understanding and simulation work is based can be found in Caldwell and 

Truini (1997).

Purpose and Scope

This report describes the mathematical simulation of regional ground-water flow in the upper 

Deschutes Basin in central Oregon. The report 

includes a description of the ground-water hydrology 

of the upper Deschutes Basin and how the hydrologic 

system was represented in the numerical model. The 

hydrologic data used for model calibration and the 

calibration procedures are also described. A discussion 

of model reliability and simulation of some hypothetical ground-water management scenarios are 

also included.

Numerical models presently provide the best 

method for synthesizing the data and the conceptual 

understanding for analyzing and predicting the behavior of a complex natural system. Numerical models 

allow the testing of conceptual models by placing 

them in a physically based mathematical framework 

constrained by data. Errors in conceptual models 

are often made obvious during numerical model 

development. If a numerical model is properly calibrated, meaning it accurately simulates the observed 

behavior of the system, then it can be used to predict 

the behavior of the system under new conditions. The 

uncertainty of the predictions, however, can be large 

if the new conditions differ markedly from the range 

of conditions over which the model was calibrated. 

The purpose of the simulation work in the upper 

Deschutes Basin was to test the reasonableness of the 

conceptual understanding of regional ground-water 

flow, provide insights into the behavior of the system, 

and evaluate the behavior of the ground-water system 

under possible future conditions. Of particular interest 

in the upper Deschutes Basin are changes in the location and amount of ground-water discharge to streams 

caused by increases in the rate of ground-water pumping and reduction of artificial ground-water recharge 

due to sealing of leaking irrigation canals. 

The simulation work described in this report 

was constrained by the limited geographic distribution 

of data and the scarcity of time-series data on such 

quantities as head and ground-water discharge. Most 

of the area of the upper Deschutes Basin is public land 

with very few wells to provide subsurface hydrologic 

information. Therefore, there are large geographic 

regions throughout which hydraulic head data and subsurface geologic data are scarce. Moreover, in the areas 

where there are wells, most penetrate only the upper 

part of the thick saturated section. Because of this, there 

is very little information on the deep parts of the flow 

system anywhere in the basin. This scarcity of data coupled with the geologic and hydrologic complexity of the 

area create inherent uncertainty in simulation results. 
3

Regardless, the simulation work described in this 

report has helped to refine the conceptual understanding of the regional flow system and provides important 

insights as to how the system will likely respond to 

changing conditions.

Study Area Description

The upper Deschutes Basin study area encompasses approximately 4,500 mi2 (square miles) of the 

Deschutes River drainage basin in central Oregon 

(fig. 1). The area is drained by the Deschutes River 

and its major tributaries: the Little Deschutes River, 

Tumalo and Squaw Creeks, and the Metolius River 

from the west, and the Crooked River from the 

east. Land-surface elevation ranges from less than 

1,300 feet near Gateway in the northern part of the 

study area to more than 10,000 feet in the Cascade 

Range. The study area includes the major population 

centers in the basin, where ground-water development 

is most intense and resource-management questions 

are most urgent. The major communities include 

Bend, Redmond, Sisters, Madras, Prineville, and 

La Pine. 

The study area boundaries were chosen to 

coincide as much as possible with natural hydrologic 

boundaries across which ground-water flow can be 

reasonably estimated or assumed to be negligible. The 

study area is bounded on the north by Jefferson Creek, 

the Metolius River, the Deschutes River, and Trout 

Creek; on the east by the generalized contact between 

the Deschutes Formation and the older, much less 

permeable John Day Formation; on the south by the 

drainage divide between the Deschutes Basin and 

the Fort Rock and Klamath Basins; and on the west 

by the Cascade Range crest.

There are approximately 164,000 acres 

(256 mi2) of irrigated agricultural land in the study 

area. The largest source of irrigation water is the 

Deschutes River. Most water is diverted from the 

Deschutes River near Bend and distributed to areas 

to the north through several hundred miles of canals. 

Smaller amounts of irrigation water are diverted 

from Tumalo and Squaw Creeks, the Crooked River, 

and Ochoco Creek.

The climate in the Deschutes Basin is controlled 

primarily by air masses that move eastward from 

the Pacific Ocean, across western Oregon, and into 

central Oregon. The climate is moderate, with cool, 

wet winters and warm, dry summers. Orographic 

processes result in large amounts of precipitation in 

the Cascade Range in the western part of the basin, 

with precipitation locally exceeding 200 in/yr (inches 

per year), mostly as snow during the winter (Taylor, 

1993). Precipitation rates diminish rapidly toward 

the east to less than 10 in/yr in the central part of the 

basin (fig. 2). Climate in the Deschutes Basin exhibits 

year to-year and longer-term variability. This variability generally parallels regional trends in the Pacific 

Northwest that have been correlated with large-scale 

ocean-atmosphere climate-variability patterns in the 

Pacific Basin such as the El Niño/Southern Oscillation 

(Redmond and Koch, 1991) and the Pacific Decadal 

Oscillation (Mantua and others, 1997). 

Acknowledgements

Numerous people contributed substantially to 

the modeling efforts in the upper Deschutes Basin 

described in this report. The authors gratefully 

acknowledge the following colleagues for their 

mentorship, insightful discussions, technical assistance, and general sharing of wisdom and ideas: 

David S. Morgan, Mary C. Hill, Leonard L. Orzol, 

Daniel T. Snyder, Henry M. Johnson, and Claire 

Tiedeman.

GROUND-WATER HYDROLOGY

Geologic Framework

 The stratigraphic and structural framework of 

the upper Deschutes Basin influences many aspects of 

ground-water recharge, flow, and discharge. Episodic 

volcanic activity in the region over the past several 

million years has resulted in a variety of volcanic, 

volcaniclastic, and volcanically derived sedimentary 

deposits (fig. 3). Volcanic activity and subsequent 

redeposition and alteration of some deposits have 

produced strata with a wide range of permeability.

Most ground water in the upper Deschutes Basin 

originates as recharge from precipitation on Quaternary deposits of the Cascade Range and Newberry 

Volcano (fig. 3). These deposits are highly permeable, 

and the fractured character of the lava flows facilitates 

rapid infiltration of precipitation and snowmelt and 

movement of ground water to lower elevations. Additional recharge occurs from leakage along sections 

of unlined irrigation canals near Bend constructed on 

lava flows from Newberry Volcano. Deposits of the 

Cascade Range and Newberry Volcano overlay or abut 

the Deschutes Formation and age-equivalent deposits.
4

Figure 1. Location of the upper Deschutes Basin, Oregon, and major geographic and cultural features. 
5

Figure 2. Lines of equal mean annual precipitation in the upper Deschutes Basin and graphs showing mean monthly precipitation at 

selected precipitation stations. 
6

Figure 3. Generalized geology and hydrogeologic units of the upper Deschutes Basin. 
7

Figure 3. Generalized geology and hydrogeologic units of the upper Deschutes Basin—Continued. 
8

Several saturated Quaternary-age sedimentary 

deposits occur within the study area (fig. 3). Pleistocene to Holocene sedimentary deposits resulted 

from several processes, including mass wasting of 

upland deposits, alluvial deposition by high-energy 

streams, deposition in low-energy streams and lakes, 

and glaciation. These sedimentary deposits are locally 

a few hundred feet thick. The La Pine subbasin occupies a complex graben structure where up to 1,000 feet 

of sediment has accumulated. 

The Deschutes Formation and age-equivalent 

deposits (fig. 3) are generally highly permeable and 

are the most widely used ground-water bearing units 

in the study area. The Deschutes Formation consists 

of a variety of volcanic and sedimentary deposits ranging in age from late Miocene (~7.5 million years) to 

Pliocene (~4.0 million years). Smith (1986) described 

distinct depositional environments and volcanic facies 

within the Deschutes Formation. The differentiated 

facies have shown to be useful hydrogeologic subdivisions. These include sedimentary deposits of an 

arc-adjacent alluvial plain, the ancestral Deschutes 

River, and an inactive basin margin, as well as 

proximal volcanic deposits that are chiefly lava flows. 

Ancestral Deschutes River deposits and some units 

within an arc-adjacent alluvial plain are among the 

highest yielding units within the Deschutes Formation, 

with some wells producing thousands of gallons per 

minute. Opal Springs basalt, Pelton basalt, and a 

rhyodacite dome complex near Steelhead Falls are 

particularly productive subunits within the Deschutes 

Formation and provide tens to hundreds of cubic feet 

per second of ground-water discharge to the Deschutes 

and Crooked Rivers, upstream of Round Butte Dam. 

The Prineville Basalt (fig. 3), which locally 

underlies the Deschutes Formation, is the oldest stratigraphic unit that bears substantial water in the upper 

Deschutes Basin. The Prineville Basalt is a series of 

middle Miocene lava flows that have a maximum 

thickness of nearly 700 feet near Bowman Dam 

(Smith, 1986). The unit occurs only locally near 

the eastern and northern boundaries of the upper 

Deschutes Basin study area. Wells constructed into 

water-bearing zones within the Prineville Basalt typically yield amounts that are inadequate for irrigation 

but are suitable for domestic and stock uses.

The oldest rock unit in the upper Deschutes 

Basin study area, the John Day Formation (fig. 3), is a 

sequence of upper Eocene to lower Miocene volcanic 

and sedimentary rocks. Weathering and alteration 

of the rocks has resulted in very low permeability; 

consequently, the unit forms the hydrologic basement 

for the regional ground-water flow system throughout 

much of the area and the lateral flow boundaries to the 

east and north. 

Hydrothermal alteration and secondary mineralization at depth below the Cascade Range and 

Newberry Volcano has drastically reduced the permeability of the material at depth in those regions, 

effectively restricting most ground water to the strata 

above the altered rocks. The top of the hydrothermally 

altered region is considered the base of the regional 

ground-water system beneath the Cascade Range and 

Newberry Volcano. These boundaries are discussed 

in more detail in subsequent sections.

Structural features influence ground-water flow 

within the upper Deschutes Basin mainly by juxtaposing materials with contrasting permeability. Several 

depositional centers have formed along the base of 

fault-line scarps or in grabens within the study area, 

and the infilling sedimentary deposits have permeability that differs from that of the surrounding rocks. 

Geologic units in the Deschutes Basin were 

divided into seven hydrogeologic units by Lite and 

Gannett (2002) (fig. 3). In some instances, the units 

correspond to existing stratigraphic divisions. In other 

instances, hydrogeologic units correspond to different 

facies within a single stratigraphic unit or formation. 

The hydrogeologic units include Quaternary sedimentary deposits, volcanic deposits of the Quaternary 

Cascade Range and Newberry Volcano, four zones 

corresponding to depositional facies within the late 

Tertiary and Quaternary Deschutes Formation and ageequivalent strata, and pre-Deschutes-age strata including the Prineville Basalt and John Day Formation.

Ground-Water Flow

The ground-water hydrology of the upper 

Deschutes Basin is described in detail by Gannett 

and others (2001), who built on the earlier works of 

Russell (1905), Stearns (1931), and Sceva (1960, 1968). 

Regional ground-water flow in the upper Deschutes 

Basin is primarily controlled by the distribution of 

recharge, the geology of the area, and the location 

and elevation of streams. Ground water flows from 

the principal recharge areas in the Cascade Range and 

Newberry Volcano toward discharge areas along the 

margin of the Cascade Range and near the confluence 

of the Deschutes, Crooked, and Metolius Rivers (fig. 4).
9

Figure 4. Generalized lines of equal head, locations of field-inventoried wells, and land ownership in the upper Deschutes Basin. 

(Generalized lines of equal head represent average conditions and are based on static water-level measurements in wells made 

during different periods and on spring and stream elevations.) 
10

At the regional scale, the distribution of recharge 

mimics that of precipitation. The annual precipitation 

rate shows considerable geographic variation throughout the upper Deschutes Basin (fig. 2). The Cascade 

Range, which constitutes the western boundary of the 

basin, locally receives in excess of 200 in/yr, mostly 

as snow. The central part of the study area, in contrast, 

typically receives less than 10 in/yr. The young 

Quaternary volcanic deposits and thin soils in the 

Cascade Range allow rapid infiltration of much of the 

rain and snowmelt, making the Cascades the locus of 

ground-water recharge for the basin. The 1993 to 1995 

average annual rate of recharge from precipitation 

basinwide is about 3,500 ft3/s (cubic feet per second) 

(Gannett and others, 2001). Precipitation provides 

relatively little ground-water recharge in the lowelevation areas in the central part of the basin; 

however, leaking irrigation canals (fig. 1) are locally 

a significant source of recharge. It is estimated that 

46 percent of the water diverted for irrigation is lost 

through canal leakage (Gannett and others, 2001). The 

average annual rate of leakage from irrigation canals 

during 1994 was estimated to be 490 ft3/s. Part of 

the ground water recharged in the Cascade Range 

discharges to spring-fed streams at lower elevations in 

the range and along margins of adjacent lowlands. 

The remaining ground water continues in the subsurface toward the central part of the basin, where 

most of it discharges to the Deschutes, Crooked, and 

Metolius Rivers in the vicinity of their confluence.

The principal ground-water fed streams along 

the topographic margin of the Cascade Range include 

the upper Deschutes River and its tributaries above 

Wickiup Reservoir, Fall and Spring Rivers farther 

downstream, and the upper Metolius River and its 

tributaries. Combined ground-water discharge along 

the margin of the Cascade Range is estimated to average about 2,600 ft3

/s, which is roughly one-half the 

total ground-water discharge of the upper Deschutes 

Basin (Gannett and others, 2001).

Stream gains and losses along the Deschutes and 

Little Deschutes Rivers in the La Pine subbasin east 

of the Cascade Range are small, indicating relatively 

little net exchange between ground water and surface 

water. North of Sunriver, however, the northward 

slope of the water table is larger than the slope of the 

land surface, so depths to ground water generally 

increase northward toward Bend (fig. 4). In this same 

reach, however, water is diverted from the Deschutes 

River into unlined irrigation canals. The combination 

of leakage from canals, and perhaps from the river 

itself, supports local shallow, possibly perched, saturated zones that leak water back to the river. Data from 

several stream gages operated along the Deschutes 

River between Benham Falls and Bend from 1944 to 

1953 (Oregon Water Resources Department, 1965) 

show that this reach contains both losing and gaining 

segments. Between Benham Falls and Lava Island, 

the Deschutes River consistently loses about 84 ft3/s 

(Gannett and others, 2001). The reach between Lava 

Island and the mill pond in Bend gains about 46 ft3

/s. 

The river then loses about 50 ft3/s between the mill 

pond and the gage below Bend. The net loss between 

Benham Falls and Bend indicated by the 1944 to 1953 

data is 88 ft3/s. This is consistent with the 89 ft3

/s loss 

indicated by long-term record (1945–95) available 

from gages at Benham Falls and below Bend.

The Deschutes and Crooked Rivers have incised 

canyons in the northern part of the study area. The 

canyons become increasingly deep northward toward 

Lake Billy Chinook, reaching depths of several hundred feet below the surrounding terrain. About 10 to 

15 miles above their confluence, the canyons of the 

Deschutes and Crooked Rivers are of sufficient depth 

to intersect the regional water table, and both streams 

gain flow from ground-water discharge. Streamflow 

measurements show that the Deschutes River and 

lower Squaw Creek combined gain about 400 ft3/s 

from ground-water discharge in this area prior to 

entering Lake Billy Chinook, and the lower Crooked 

River gains about 1,100 ft3/s before entering the lake. 

Ground-water discharge to Lake Billy Chinook is 

roughly 420 ft3/s. The total ground-water discharge 

in the confluence area is approximately 2,300 ft3/s. 

This ground-water discharge, along with the 1,200 

to 1,400 ft3/s flow of the Metolius River (which is 

predominantly ground-water discharge during the 

dry season), makes up virtually all the flow of the 

Deschutes River at Madras during the summer and 

early fall.

Geologic factors are the primary cause of the 

large ground-water discharge in the confluence area. 

The permeable Neogene deposits, through which 

nearly all regional ground water flows, become 

increasingly thin northward as the low-permeability 

John Day Formation nears the surface. The John 

Day Formation crops out in the Deschutes River 

canyon about 10 miles north of Lake Billy Chinook 

near Pelton Dam, marking the northern extent of 

the permeable regional aquifer system. Most of the 
11

regional ground water in the upper basin discharges 

to the Deschutes and Crooked Rivers south of this 

location. There is no appreciable ground-water discharge directly to the Deschutes River downstream 

of this point, and the small gains in streamflow that 

do occur result primarily from tributary inflow.

Ground-water head elevations and ground-water 

discharge fluctuate in the Deschutes Basin primarily 

due to decadal climate cycles (Gannett and others, 

2001). Ground-water levels in the Cascade Range 

fluctuate in excess of 20 feet in response to climate 

cycles. Ground-water discharge varies similarly in 

response. For example, August mean ground-water 

discharge to the Deschutes River between Bend and 

Culver varied over 100 ft3

/s between 1962 and 1997 

due to climate cycles. Ground-water discharge to the 

Metolius River, based on October mean flows, varied 

over 400 ft3/s from 1962 to 1997. Combined, climatedriven ground-water discharge fluctuations could 

account for variations in late-season monthly mean 

flows of the Deschutes River at Madras on the order 

of 1,000 ft3/s (Gannett and others, 2001). Head elevations also fluctuate in response to local stresses such as 

canal operation and stream-stage variations. 

Hydraulic Characteristics

The hydraulic characteristics considered in this 

section include horizontal and vertical hydraulic conductivity and storage coefficients. Sources of information on the hydraulic characteristics of subsurface 

materials in the Deschutes Basin include water-well 

reports (well logs), aquifer tests, and the literature. 

Information from each of these sources is discussed in 

the following paragraphs. Both horizontal and vertical 

hydraulic conductivity are considered. In the following discussion, the term “hydraulic conductivity” 

with no modifier refers to horizontal hydraulic conductivity. Hydraulic conductivity in the vertical direction is specifically referred to as vertical hydraulic 

conductivity.

Horizontal Hydraulic Conductivity

There is no evidence that horizontal hydraulic 

conductivity varies with direction at scales of concern in this study. Therefore, horizontal isotropy is 

assumed. Well logs filed with the State by drillers 

typically contain specific-capacity tests that can provide information on hydraulic-conductivity values. 

The methodologies for determining hydraulic-conductivity values from specific-capacity tests contain 

numerous assumptions that are seldom met, so the 

values are, at best, rough estimates. Regardless, 

specific-capacity tests are often used to determine 

patterns of spatial distribution of hydraulic conductivity for models.

Specific-capacity tests from well logs were of 

only limited use in the upper Deschutes Basin for 

two principal reasons. The first problem is that the 

distribution of wells is uneven both horizontally and 

vertically, with the majority of wells clustered in the 

shallowest parts of the flow system in the center part 

of the study area (fig. 4). The second problem stems 

from variability in the data. The specific capacity of 

a single well can vary with the pumping rate due to 

well inefficiency. Moreover, apparent specific capacities of different wells in the same aquifer can vary due 

to geologic heterogeneity, different well completion 

techniques, and pump-test methods. This means that 

there can be large variations in specific-capacity data 

from a single geologic unit. These variations are not 

randomly distributed spatially because controlling 

factors, such as well types and test methods, are not 

randomly distributed but may vary systematically 

with, for example, land-use patterns or working 

areas of particular drillers. The uneven distribution 

and large variability in the specific-capacity data 

precluded their use for creating a detailed map of 

hydraulic-conductivity distribution. These data were 

useful, however, in evaluating large-scale patterns in 

hydraulic conductivity.

Gannett and others (2001) presented an analysis 

of specific-capacity data from well logs in which only 

domestic wells were considered in order to minimize 

the variance that would have been introduced by 

including large-yield irrigation and municipal wells. 

Large-scale patterns, on the order of 10 miles, appear 

to be consistent with geologic factors occurring at 

equivalent scales. Specific-capacity data show that 

within the Deschutes Formation, hydraulic conductivities are relatively low along the eastern margin of 

the basin in areas roughly correlative with the inactive 

margin facies of Smith (1986) (fig. 3). Data also show 

that, in general, hydraulic conductivities are lower 

in the La Pine subbasin aquifers than in the Deschutes 

Formation to the north.

Aquifer tests are another source of hydraulicconductivity estimates. Four tests conducted as part 

of this study, as well as several preexisting tests, were 
12

analyzed to estimate hydraulic conductivity in parts 

of the basin. In addition, several tests conducted by 

consultants were also available. The aquifer tests are 

summarized in Gannett and others (2001). Many of 

these tests are discussed in the context of geologic units 

in Lite and Gannett (2002). The analyses of most tests 

were complicated by very high transmissivity values 

that resulted in very small drawdowns. Drawdowns 

in observation wells were commonly too small to 

resolve, even with pumping rates on the order of 100 

to 1,000 gal/min (gallons per minute). 

The small number of aquifer tests precludes 

developing a detailed spatial distribution of hydraulic 

conductivity. However, the tests provide good estimates of hydraulic characteristics of certain rock types 

at specific locations, provide insight into the likely 

conditions in similar geologic settings in the basin, and 

constrain the overall range of probable values. Tests in 

volcanic deposits dominated by basaltic lava and scoria 

of the Deschutes Formation and age-equivalent units 

yielded hydraulic-conductivity estimates of 14 ft/d 

(feet per day) to 2,300 ft/d (1.6 × 10-4 to 2.7 × 10-2 ft/s 

[feet per second]). Hydraulic-conductivity estimates 

from tests in sedimentary strata ranged from 16 to 

150 ft/d (1.9 × 10-4 to 1.7 × 10-3 ft/s). Aquifer testing 

of the City of Redmond Well Number 3, which produces from both volcanic and fluvial strata, yielded 

estimates ranging from 1,500 to 2,300 ft/d (1.7 × 10-2

to 2.7 × 10-2 ft/s). Although the small number of tests 

precluded any detailed assessment of the spatial distribution of hydraulic conductivity, aquifer-test results 

are consistent with the regional pattern in estimates 

derived from specific capacity, described previously, 

and with expectations based on the distribution of geologic facies within the basin.

The published literature is the third source of 

information on hydraulic characteristics used in this 

study. Literature values include those presented in general references as typical for particular rock types and 

those specific to the Deschutes Basin. Values from 

general references are usually presented as broad 

ranges, often spanning several orders of magnitude, for 

particular rock types. The ranges are usually too broad 

to be of much use, other than indicating the range of 

possible values. For example, Freeze and Cherry 

(1979, p. 29) show a range of 5 orders of magnitude 

for “permeable basalt”: 4 × 10-2 to 6 × 103 ft/d 

(5 × 10-7 to 7 × 10-2 ft/s). If unfractured basalt is 

included, the lower limit is 10-8 ft/d (10-13 ft/s),

extending the overall range to 11 orders of magnitude. 

Hydraulic-conductivity values given for sedimentary 

deposits ranging from silt to gravel span 9 orders 

of magnitude, from approximately 2.6 × 10-4 to 

2.6 × 105 ft/d (3 × 10-9 to 3 ft/s). 

Outside of this study, there are few published 

hydraulic-conductivity estimates for the upper Deschutes Basin, and these are average values given to 

represent regions within the basin. In work covering all 

of eastern Oregon, Gonthier (1985) estimated hydraulic conductivity for two generalized geologic units 

in the upper Deschutes Basin (basin-fill and alluvial 

deposits, and volcanic and sedimentary deposits) using 

the literature and specific-capacity data. Ranges of 

estimated hydraulic conductivity given are 25 to 

150 ft/d (2.9 × 10-4 to 1.7 × 10-3 ft/s) for the basin-fill 

and alluvial deposits and 10 to 500 ft/d (1.2 × 10-4

to 5.8 × 10-3 ft/s) for the volcanic and sedimentary 

deposits. Bolke and Laenen (1989) estimated hydraulic 

conductivity from specific-capacity data from the 

Warm Springs Indian Reservation. They estimate that 

the hydraulic conductivity of alluvium ranges from 

100 to 1,000 ft/d (10-3 to 10-2 ft/s) and give a range 

of 0.1 to 0.01 ft/d (10-7 to 10-8 ft/s) for John Day 

Formation strata. Lite and Gannett (2002) described 

the regional geology and volcanic-facies distribution 

in the upper Deschutes Basin and placed the overall 

large-scale hydraulic-conductivity distribution in a 

geologic context. 

There are very few direct measurements of 

hydraulic conductivity in the Cascade Range, and most 

knowledge of the large-scale permeability structure 

there is inferred from geothermal-gradient measurements and heat-flow studies. In simulating groundwater flow and heat transport in the Cascade Range, 

Ingebritsen and others (1992) estimated the permeability of rocks younger than 2.3 Ma (mega annum, 

or million years before present) to be about 10-13 ft2

(10-14 m2), which is equivalent to a hydraulic conductivity of about 0.018 ft/d (2.1 × 10-7 ft/s) assuming a 

water temperature of 5°C. The permeability of rocks 

with ages between 4 and 8 Ma was estimated to be 

5.4 × 10-15 ft2 (5.0 × 10-16 m2) which is equivalent 

to a hydraulic conductivity of about 9.1 × 10-4 ft/d 

(1.1 × 10-8 ft/s). Higher near-surface permeability, 

on the order of 10-14 to 10-12 m2 (.018 to 1.8 ft/d), 

was required in their simulation to match groundwater recharge estimates. Higher near-surface permeabilities are also suggested by well-test data. A 

specific-capacity test of a well near Mount Bachelor 

(18S/09E-20BDA) yielded a hydraulic-conductivity 

estimate of 9 ft/d (1.0 × 10-4 ft/s).
13

Blackwell and Priest (1996) suggested that the 

heat-transport pattern in the Cascade Range is predominantly conductive, indicating that the permeability in all but the youngest rocks below a depth of 

about 330 to 990 feet (100 to 300 meters) is orders 

of magnitude lower than the estimates of Ingebritsen 

and others (1992). They also suggested that groundwater flow at velocities sufficient to affect heat flow is 

restricted to local regions, except in the 0 to 1,600 feet 

(0 to 500 meters) depth range in the high Cascade 

Range. The volume of water moving through lowpermeability strata at depths greater than 1,600 feet 

in the Cascade Range is sufficiently small to be considered negligible compared to the overall groundwater budget, and these low-permeability strata are 

considered in this study to be the base of the regional 

ground-water flow system. The younger, permeable 

near-surface strata in the Cascade Range are included 

in the regional ground-water flow model. 

Mathematical modeling of ground-water discharge to spring-fed streams in the Cascade Range 

by Manga (1996, 1997) yielded permeability values 

for near-surface rocks less than about 2.0 Ma of about 

10-10 ft2 (10-11 m2), which equates to a hydraulic 

conductivity of about 18 ft/d (2 × 10-4 ft/s) assuming a 

water temperature of 5°C. This estimate is an order of 

magnitude larger than the upper value of Ingebritsen 

and others (1992) for near-surface rocks, where most 

ground-water flow occurs. The permeability estimates 

of Manga (1996, 1997) and Ingebritsen and others 

(1992) are considered to be a reasonable range of 

values for the younger, near-surface strata in the 

Cascade Range.

Well- and aquifer-test data are insufficient to 

allow continuous mapping of the distribution of 

hydraulic conductivity in the upper Deschutes Basin. 

Lite and Gannett (2002) delineated seven hydrogeologic units based on stratigraphy, facies distribution, 

and large-scale rock properties. These seven units 

provide a context in which to apply the aquifer-test 

data described above and provide the basis for the 

hydraulic-conductivity distribution in the model. 

Vertical Hydraulic Conductivity

There were no direct measurements or estimates 

of vertical hydraulic conductivity available for the 

Deschutes Basin prior to this study. The only vertical 

hydraulic-conductivity estimates available for volcanic terranes in the Pacific Northwest were derived 

from modeling studies in the Columbia River Basalt 

Group and related strata. 

Vertical and horizontal hydraulic-conductivity 

values are commonly related. The ratio of horizontal 

to vertical hydraulic conductivity, Kxy / Kz, is termed 

the vertical anisotropy. Values of vertical anisotropy 

vary widely with geologic settings, although vertical 

hydraulic conductivity is almost always less than 

horizontal. Anderson and Woessner (1992) state that 

vertical anisotropy ratios commonly range from 1 to 

1,000. Studies in the Columbia River Basalt Group, 

however, suggest that vertical anisotropy ratios may 

be considerably larger than 1,000 in some volcanic 

terranes. Hansen and others (1994) found that vertical 

anisotropy ratios of 500 to 3,300 gave the best overall 

results in simulations of regional ground-water flow 

in the Columbia Plateau in Washington and northern 

Oregon. Packard and others (1996) found that vertical 

hydraulic-conductivity values ranging from 4.92 × 10-4

to 1.08 × 10-2 ft/d (5.7 × 10-9 to 1.25 × 10-7 ft/s) gave 

best simulation results, and that vertical anisotropy 

ratios ranged from 43 to 2,500. They related vertical 

anisotropy ratios to structural settings in their model 

area, and attributed the variations to differing interbed 

lithologies. Tanaka and others (1974) found that their 

best vertical head distributions were achieved using a 

vertical hydraulic-conductivity value for Columbia 

River basalt of 2 × 10-5 ft/d (2.3 × 10-10 ft/s). Although 

it is unclear how this figure was applied in their model 

parameterization, when compared to the range of horizontal transmissivity values they used it suggests that 

vertical anisotropy ratios may be as high as 104 to 106. 

Mac Nish and Barker (1976) suggest that horizontal 

hydraulic conductivity in the Columbia River basalt 

may be as large as 432,000 ft/d (5 ft/s) and that 

vertical hydraulic conductivity may be as small as 

4.32 × 10-3 ft/d (5 × 10-8 ft/s), and that vertical 

anisotropy ratios may approach 109. No supporting 

information or references, however, are provided 

for these figures. Very large vertical hydraulic-head 

gradients occur in places in the upper Deschutes Basin, 

suggesting that large vertical anisotropy ratios occur at 

least locally. 

Storage Coefficients

Measurements of storage coefficients of aquifers 

in the upper Deschutes Basin are scarce. This is not 

only because there are few aquifer tests, but the pumping responses in observation wells necessary to calculate storage coefficients are commonly too small to be 
14

useful due to the large transmissivity values of the 

tested aquifers. Storage coefficients are available from 

two aquifer tests conducted as part of this study. The 

City of Madras Well 2, which produces from a confined sand and gravel aquifer, yielded a value of 10-4. 

The test of the City of Redmond Well 3, which produces from both lava and interbedded sand and gravel, 

yielded values ranging from 5 × 10-2 to 6 × 10-2. 

A third test conducted by a private consultant on the 

City of Redmond Well 4, which produces from sand 

and gravel, yielded a storage coefficient of 0.1. 

MODEL DESCRIPTION

Ground-water flow in the upper Deschutes Basin 

was simulated using the U.S. Geological Survey modular three-dimensional finite-difference ground-water 

flow model (MODFLOW) developed by McDonald 

and Harbaugh (1988). MODFLOW is a versatile 

finite-difference modeling program used to construct 

numerical flow models of specific areas. A MODFLOW model consists primarily of a set of input files 

that contain information on the physical properties of 

the modeled system such as the geometry, boundary 

conditions, internal properties (such as the distribution 

of hydraulic conductivity and storage coefficients), 

and sources and sinks such as ground-water recharge, 

streams, and pumping wells. Once these files are 

created, the model program is run to solve a set 

of equations that describe the distribution of head 

at discrete points within the system and the flow in 

response to that head distribution. 

Both steady-state and transient versions of the 

upper Deschutes Basin model were constructed. The 

steady-state model simulates the ground-water flow 

system as it would exist in equilibrium with a given 

set of conditions. The transient model simulates 

the ground-water flow system as it would change in 

response to varying conditions. 

Governing Equations and Model Code

The movement of ground water through porous 

media is described by the following partial differential 

equation, which is based on Darcy’s law and the conservation of mass (McDonald and Harbaugh, 1988):

 (1) ∂x

∂ Kxx∂x

∂h

⎝ ⎠ ⎛ ⎞

∂y

∂ Kyy∂y

∂h

⎝ ⎠ ⎛ ⎞ +

∂z

∂ Kzz∂z

∂h

⎝ ⎠ + ⎛ ⎞ – W Ss

∂t

∂h =

where 

Kxx, Kyy, and Kzz are values of hydraulic conductivity in the x, y and z directions along 

Cartesian coordinate axes, which are assumed 

to align with principal directions of hydraulic 

conductivity (LT-1),

h is hydraulic head (L),

W is a volumetric flux per unit volume and 

represents sinks and/or sources (T-1), 

Ss is the specific storage of the porous material 

(L-1), and 

t is time (T). 

Equation (1) describes the distribution of 

hydraulic head and flow throughout a continuous 

region. Derivations of equation (1) can be found 

in Freeze and Cherry (1979) and Anderson and 

Woessner (1992). As mentioned previously, there is 

no evidence of horizontal anisotropy in the upper 

Deschutes Basin, therefore Kxx and Kyy are considered 

to be equal at any given location and Kxx and Kyy are 

replaced in this discussion by the single term Kxy to 

describe horizontal hydraulic conductivity.

Equation (1) is continuous in space and time, and 

generally cannot be solved analytically for practical 

applications involving complex systems. In practice, 

numerical methods are employed in which the continuous system described by equation (1) is replaced by 

a set of spatially and temporally discrete points in a 

process known as discretization. Equation (1) is then 

replaced by a set of simultaneous algebraic equations 

that describe the distribution of hydraulic head at each 

point, and flow through the system in response to this 

head distribution. These simultaneous equations are 

set up in matrix form and then solved. A variety of 

techniques are available to solve the set of simultaneous equations. In this application, the preconditioned conjugate-gradient method of Hill (1990) was 

used. Discussions of the numerical technique used in 

this study, the finite-difference method, can be found 

in McDonald and Harbaugh (1988) and in Anderson 

and Woessner (1992).

Discretization

The numerical method used to approximate 

equation (1) requires that the modeled domain be 

divided into discrete volumes, called cells. The threedimensional array of cells is known as the model grid. 

The center of each cell defines the point for which 
15

hydraulic head is determined. The head is taken to represent the average head within the cell. For transient 

models, time must also be divided into discrete intervals known as time steps. Heads and flows are calculated at the end of each time step. 

Spatial Discretization

The regional ground-water flow system of the 

upper Deschutes Basin was represented as an array 

of cells arranged in 127 east-west trending rows, 87 

north-south trending columns, and 8 layers. The lateral 

dimensions of the cells range from 2,000 feet to 

10,000 feet on a side (fig. 5). The smaller cells coincide with areas where the most hydrologic data are 

available; generally, the more densely populated parts 

of the study area, where wells are numerous. The area 

of smallest cells also corresponds to the main area of 

ground-water discharge to streams, which is an area 

of particular interest to resource managers. The larger 

cells coincide with areas where hydrologic data are 

scarce; generally, unpopulated areas, where few wells 

have been drilled. Because the finite-difference model 

grid must be rectangular, certain cells may represent 

areas outside the modeled area (which can have any 

shape). Such cells are considered inactive and the 

ground-water flow equation is not solved for them.

Eight layers were used in order to represent 

vertical changes in the geology and allow simulation 

of vertical head gradients. Each layer has a uniform 

thickness throughout the model. The upper five layers 

are 100 feet thick; underlying layers, in order of 

increasing depth, are 200, 300, and 800 feet thick. The 

top of the uppermost layer is coincident with the water 

table. As with the horizontal dimension, the vertical 

layer thicknesses are smallest where data are most 

plentiful, and largest where data are sparse. Near the 

lateral and lower boundaries of the model, where 

the aquifer system may be thin, the thicknesses of the 

lowermost layer may be less in places than indicated 

above. For numerical stability, all layers were simulated as though they were confined, meaning that 

saturated thickness, and consequently transmissivity, 

were held constant with respect to head. In reality, it 

is possible for saturated thickness and transmissivity 

to vary in the uppermost parts of the flow system as 

head changes. However, water-level monitoring data 

indicate that head fluctuations in the basin are small 

compared to model layer thicknesses, so the error 

introduced by modeling all layers as confined is considered insignificant.

Temporal Discretization

For transient models, numerical methods require 

that time also be divided into discrete intervals. The 

temporal divisions are determined by the timing of the 

stresses, such as recharge or pumping. Time is divided 

into stress periods, during which the stresses are 

assumed to remain constant. For example, the irrigation season, during which irrigation pumping occurs, 

may define a stress period. Stress periods are further 

subdivided into time steps as described below. The 

ground-water flow equations in the model are solved 

for each time step. For this study, ground-water flow 

in the upper Deschutes Basin was simulated for the 

20-year period extending from 1978 to 1997. Each 

year was divided into two stress periods that were 

determined by the timing of recharge and irrigation 

pumping. The delineation of the stress periods is 

described in a later section. Each of the 40 stress 

periods was divided into 5 time steps. The 20-year 

simulation period, therefore, included 200 time steps.

Boundary Conditions

The ground-water flow system in the upper 

Deschutes Basin is controlled, to a large extent, 

by the hydrologic and geologic boundaries of the 

system. Boundary conditions define the geographic 

extent of the flow system as well as the movement of 

ground water into and out of the system, such as flow 

to or from streams.

Three types of boundary conditions were used 

to define the ground-water flow system in the upper 

Deschutes Basin: no-flow boundaries, specified-flux 

boundaries, and head-dependent flux boundaries. 

Geologic or hydrologic barriers to ground-water flow 

were simulated using no-flow boundaries. The contact 

between the permeable ground-water flow system 

and nearly impermeable bedrock is an example of a 

no-flow boundary. Known or estimated hydrologic 

fluxes, such as recharge and well discharge, are represented using specified-flux boundaries. A headdependent flux boundary is one across which ground 

water moves at a rate proportional to the hydraulichead gradient between the boundary and the groundwater system. Streams are usually represented as headdependent flux boundaries because the movement of 

ground water to or from a stream is proportional to 

the difference between the head in the aquifer and the 

stage of the stream.
16

Figure 5. Model grid, boundary conditions, locations of stream cells, and regions where evapotranspiration is active. 

(ET, evapotranspiration)
17

For the purposes of this discussion, boundaries 

in the upper Deschutes Basin regional ground-water 

flow model are divided into two broad categories. 

The first group includes boundaries related to the 

geographic extent of the regional aquifer system, such 

as contacts with effectively impermeable rock units on 

the margin of the basin. The second group includes 

boundaries within the model extent related to hydrologic processes or features, such as stream systems. 

Geographic Boundaries

The geographic boundaries of the upper 

Deschutes Basin regional ground-water model were 

chosen to correspond as closely as possible with 

natural hydrologic boundaries across which groundwater flow can be assumed negligible, such as groundwater divides, or can be reasonably estimated. Major 

topographic divides are often considered no-flow 

boundaries because topographic divides are typically 

coincident with ground-water divides. Ground water 

on either side of a ground-water divide flows away 

from the divide and not across it, so the divide itself 

acts as a no-flow boundary. Topographic divides often 

coincide with ground-water divides because upland 

areas commonly have larger amounts of precipitation 

and recharge than surrounding areas, so the watertable surface develops a coincident high-elevation 

region from which ground-water flow diverges. 

Major streams can also act as no-flow boundaries, particularly if they are regional ground-water 

discharge points. In these circumstances, ground water 

flows toward and discharges to the stream from both 

sides but does not flow across the vertical plane 

extending downward from the stream. Streams may 

also coincide with no-flow boundaries when they 

flow parallel to ground-water flow, basically following ground-water flow streamlines. In situations such 

as this, ground water would not flow either toward or 

away from the stream. In some circumstances, however, the hydraulic-head distribution at depth beneath 

streams may be such that there is a component of 

regional ground-water flow at depth across the stream. 

Therefore, streams should be considered no-flow 

boundaries only when there is hydraulic-head information supporting their function as regional sinks, or 

when their orientation is parallel to ground-water flow.

The lateral boundaries of the upper Deschutes 

Basin flow model are generally represented as no-flow 

boundaries, with the exception of two areas where 

they are represented as specified-flux boundaries 

(fig. 5). The western boundary of the model coincides 

with the topographic crest of the Cascade Range. 

The distribution of precipitation along this western 

margin generally mimics the topography, with the 

highest precipitation occurring at the crest of the 

range. In the northern part of the western margin 

(the Metolius River subbasin), however, the area of 

maximum precipitation occurs west of the topographic 

divide (see Taylor, 1993). This suggests that there may 

be ground-water flow eastward beneath the crest of the 

Cascade Range, a phenomenon also indicated by the 

hydrologic budget of the Metolius River subbasin. 

Average ground-water discharge to the Metolius River 

in the study area above the gage near Grandview 

(station number 14091500, which is just upstream 

of Lake Billy Chinook) is approximately 1,300 ft3/s, 

whereas the mean annual recharge from precipitation 

in the Metolius River subbasin above this point in 

the study area is approximately 500 ft3/s (Gannett 

and others, 2001). The difference, 800 ft3/s, must 

come from subsurface flow from an adjacent basin. 

Hydraulic-head gradients shown by Gannett and 

others (2001) indicate that the most plausible sources 

for this additional water are the upper Santiam and 

the North Santiam River Basins to the west. To simulate this probable interbasin ground-water flow and 

balance the hydrologic budget, a constant 800 ft3/s 

flux is specified at the western boundary of the model 

in the Metolius River subbasin.

The northern boundary to the model is defined 

by the east-west segment of the Metolius River and the 

generally east-west trending Jefferson Creek. These 

streams generally parallel regional ground-water flow 

directions, and ground-water flow into or out of the 

model area from the adjacent area is considered to 

be negligible. This part of the model boundary is 

represented as a no-flow boundary. 

North of the confluence of the Metolius and 

Deschutes Rivers, the western boundary of the model 

is defined by the Deschutes River. Hydraulic-head 

measurements and streamflow measurements by 

Bolke and Laenen (1989) indicate that the Deschutes 

River is a regional ground-water sink in this area. 

North of Pelton Dam, the river has cut down to the 

nearly impermeable John Day Formation. Regional 

ground-water flow across the Deschutes River north 

of the Metolius River confluence is unlikely, and it is 

therefore modeled as a no-flow boundary.
18

The eastern boundary of the model to south 

of Prineville Reservoir generally coincides with the 

contact between the permeable Deschutes Formation 

and older rocks of the John Day Formation which are, 

in general, nearly impermeable. The part of the boundary coinciding with the contact with John Day Formation is modeled as a no-flow boundary. 

From south of Prineville Reservoir, through 

Millican and the China Hat area, the eastern model 

boundary does not coincide with either a topographic 

divide or a geologic contact. The region east of this 

area was excluded from the model because of the 

lack of subsurface hydrologic information, very low 

recharge, and distance from the areas of primary 

concern. Head data, however, indicate there is some 

flow across this boundary into the modeled area from 

the southeast. Therefore, this area was represented as a 

specified-flux boundary (fig. 5). The boundary flux 

was estimated using a variety of methods. The part of 

the Deschutes Basin east of this boundary is very dry 

(10 to 15 in/yr average annual precipitation) and has a 

poorly developed drainage system with no perennial 

streams. The divide between this part of the Deschutes 

Basin and the Fort Rock and Christmas Lake Basins 

to the south is poorly defined and interbasin flow is 

likely. Miller (1986) states that flow to the Deschutes 

Basin from the Fort Rock Basin “probably exceeds 

10,000 acre-ft/yr” (acre feet per year) which equals 

about 14 ft3

/s. Estimates based on the Darcy equation 

using measured head gradients, and estimated hydraulic conductivity and aquifer thickness suggest that the 

flux into the modeled area could possibly be as high 

as 100 ft3/s. Additional estimates were derived using a 

water-budget approach. The probable area contributing to the boundary flux was defined using hydraulichead maps from the Deschutes Basin and the Fort 

Rock Basin from Miller (1986). Flux rates were 

calculated using a range of recharge values from 

Newcomb (1953), Miller (1986), and McFarland and 

Ryals (1991). Using a contributing area of 648 mi2

and recharge estimates ranging from 0.5 to 3.0 in/yr, 

the boundary flux could range from 25 to 145 ft3/s. 

The boundary flux in this region was ultimately set at 

a constant 48 ft3

/s (the figure using 1.0 in/yr recharge 

in the contributing area) because it provided the 

best fit to the head data during preliminary model 

calibration.

The southern model boundary coincides entirely 

with the drainage divides between the Deschutes 

Basin and the Fort Rock and Klamath Basins. It is 

represented as a no-flow boundary. 

The lower boundary of the model corresponds to 

the contact of the permeable rocks of the Deschutes 

Formation, age-equivalent units, and younger strata, 

with underlying low-permeability rock units. The 

underlying low-permeability rocks correspond to either 

the John Day Formation or equivalent rocks, or to 

younger rocks in which the permeability has been 

nearly eliminated by hydrothermal alteration or secondary mineralization. The subsurface permeability 

structure is discussed in detail in Lite and Gannett 

(2002). This lower boundary of the model is simulated 

as a no-flow boundary.

Hydrologic-Process Boundaries

The second broad category of boundaries in the 

model are those related to hydrologic processes including recharge, ground-water flow to and from streams, 

evapotranspiration from the water table, and pumpage.

Ground-water recharge constitutes a specifiedflux boundary condition. Recharge to the ground-water 

system from infiltration of precipitation, canal leakage, 

and deep percolation of applied irrigation water (onfarm loss) is simulated as a specified flux to the uppermost layer of the model. The specified flux varies from 

cell to cell. Recharge also varies with time, and is the 

largest factor influencing the transient behavior of 

the regional ground-water flow system. Therefore, the 

specified recharge flux to each uppermost cell is varied 

with time during the transient simulation. 

The methods used to estimate recharge from all 

sources are described in detail in Gannett and others 

(2001). Recharge from precipitation was estimated 

using the Deep Percolation Model (DPM) of Bauer 

and Vaccaro (1987), which uses a variety of data to 

calculate a daily water balance. This technique yields 

daily estimates of recharge, but practical considerations require that these be aggregated into longer time 

periods. A graph of average monthly recharge from 

precipitation for the basin (fig. 6) shows a period of 

high recharge from November to April, and a period 

of low recharge from May to October. 

Recharge from canal leakage was estimated by 

Gannett and others (2001) using diversion and delivery 

records from irrigation districts, canal-seepage 

measurements from the Oregon Water Resources 

Department and Bureau of Reclamation, and information on canal geology and geometry. Recharge from 

on-farm loss was estimated using information on irrigation deliveries, prevalent irrigation techniques, croptype distribution, and evapotranspiration estimates. 
19

Recharge from canal leakage occurs primarily during 

the irrigation season from April to October, with 

the highest rates from May to September (fig. 6). 

On-farm losses were determined by comparing consumptive use figures for specific crop types with 

information on water application rates and irrigation 

method efficiencies (Gannett and others, 2001). The 

timing of recharge from on-farm loss is similar to that 

of canal leakage, and is essentially the inverse of the 

timing of recharge from precipitation.

Figure 6. Mean monthly recharge from precipitation, canal 

leakage, and on-farm loss for 1978 to 1997. 

The temporal distribution of recharge was 

used to determine appropriate stress periods for the 

transient simulation. The timing of recharge from 

precipitation and canal leakage define two distinct 

periods during each year: a period of high recharge 

from precipitation and low recharge from canal leakage from November to April, and a period of lower 

recharge from precipitation and high recharge from 

canal leakage from May to October. Pumping stresses 

in the basin also match the semiannual stress periods. 

Municipal pumping is lowest and irrigation pumping 

is effectively zero during the November to April 

period. Municipal pumping increases and irrigation 

pumping occurs during the May to October period. 

In this report, we define a model water year that runs 

from November 1 to October 31. Recharge rates from 

precipitation, canal leakage, and on-farm loss are 

shown for each semiannual stress period from model 

water years 1978 to 1997 in figure 7. 

Figure 7. Recharge from precipitation, canal leakage, and on-farm loss for each semiannual stress period from 1978 to 1997. 
20

The movement of ground water to and from 

streams is another important boundary condition. 

The flux of water between the ground-water system 

and streams is, in part, dependent on the hydraulic 

head in the ground-water system and is simulated as a 

head-dependent flux boundary. Cells in the model that 

correspond to the locations of major streams (fig. 5) 

are mathematically represented in a manner that 

allows ground water to move between the aquifer and 

stream with a direction and magnitude that depends on 

the head relations. Certain stream reaches were simulated in model layers 2, 3, or 4 where the stream elevation was substantially below that of the surrounding 

terrain, for example in an incised canyon (fig. 5). 

Streams were simulated using the streamflow-routing 

(STR1) module for MODFLOW by Prudic (1989). 

The formula for calculating the flux of water 

between the ground-water system and streams is based 

on Darcy’s law and is given by Prudic (1989) as 

(2)

where 

Ql is leakage to or from the aquifer through 

 the streambed, (L3/T), 

Hs is head in the stream, (L), 

Ha is head in the aquifer, (L), and 

CSTR is conductance of the streambed (L2/T), 

 which is the hydraulic conductivity of 

 the streambed multiplied by the area of 

 the streambed divided by its thickness. 

The rate of ground-water flow to or from streams 

is proportional to the difference between the head 

in the stream and the head in the aquifer, and to the 

streambed conductance. If the head in the aquifer is 

greater than the head in the stream, then Ql

 is negative, 

indicating flow from the aquifer to the stream; if the 

head in the aquifer is less than the head in the stream, 

then the flow is from the stream to the aquifer. 

Cells in the model representing stream-aquifer 

relations (stream cells) are shown in figure 5. These 

include all major perennial streams and most perennial 

tributaries. Some stream cells actually represent 

multiple tributaries where the scale of the cell is large 

compared to the spacing of the streams. Parts of the 

Crooked River that flow on low-permeability bedrock 

are very poorly connected to the regional groundwater system and are not represented in the model. 

The movement of ground water to or from streams 

is a function of the head in the stream, Hs, also known 

as stream stage. Stream stage is the elevation of the 

water surface in the stream, and a single value must 

be chosen to represent the stream across the entire 

cell. Stream stages used in this simulation were determined by picking an approximate median stream 

elevation within the cell using USGS 1:24,000-scale 

topographic maps. 

The streambed conductance, CSTR, controls 

the rate at which water moves to or from a stream 

in response to a given head gradient. Direct measurements of streambed conductance are rare and 

these terms are usually derived empirically during 

model calibration. As generally conceptualized, 

the streambed-conductance term is used to simulate 

a low-permeability streambed material that impedes 

the movement of water to or from the stream 

(McDonald and Harbaugh, 1988). This is a reasonable representation of losing streams in the upper 

Deschutes Basin, where sediment on the bottom of a 

stream and in fractures below the stream impede water 

leakage from the stream. It is not, however, a reasonable representation of many of the gaining reaches 

where water flows to the stream from springs emerging from open fractures at or above the elevation of the 

stream. In this case, there is no low-permeability layer 

impeding the exchange of water and the rate is controlled primarily by the permeability of the surrounding rock. Streambed-conductance terms were initially 

calculated using a streambed thickness of 10 feet, 

stream widths of 20 to 50 feet (depending on the 

stream size), and a vertical hydraulic conductivity of 

0.5 to 5.0 ft/d. The hydraulic-conductivity value was 

then adjusted during model calibration.

Evapotranspiration can remove water from 

both the unsaturated and saturated zones. Evapotranspiration from the unsaturated zone is calculated 

and accounted for by the DPM. Evapotranspiration 

from the saturated zone occurs where certain plants 

(phreatophytes) transpire water directly from the water 

table. This occurs generally in areas where the water 

table is relatively close to land surface and within the 

rooting depth of plants. Evapotranspiration (ET) from 

the water table is represented by a head-dependent 

flux boundary condition in the model. 

Ql CSTR Hs Ha = ( ) –
21

Evapotranspiration from the water table is simulated using the evapotranspiration (EVT) module 

for MODFLOW (McDonald and Harbaugh, 1988). 

Evapotranspiration is calculated assuming that when 

the head in the aquifer (the water-table elevation) 

is at or above a specified elevation known as the ET 

surface (generally taken to be land-surface elevation), 

evapotranspiration from the water table will occur 

at a specified maximum rate. When the water-table 

elevation falls below a specified depth, known as the 

extinction depth (generally corresponding to the maximum rooting depth of plants), evapotranspiration from 

the water table stops. When the water-table elevation 

is between the ET surface elevation and the extinction 

depth, evapotranspiration varies linearly from the 

maximum rate to zero depending on the water-table 

elevation. 

The maximum evapotranspiration rate used in 

this simulation was based on information from the 

DPM. The DPM uses information on vegetation, soils,

elevation, and climate to calculate the total potential 

evapotranspiration of an area as well as the actual 

evapotranspiration from the unsaturated zone. The difference between these two figures is the portion of the 

potential evapotranspiration not satisfied by evapotranspiration from the unsaturated zone. This residual 

evapotranspiration is the amount that could potentially 

be satisfied by evapotranspiration from the water table 

and is the figure used for the maximum evapotranspiration rate in the EVT module. The maximum evapotranspiration rate for the steady-state model was set 

at 2 ft/yr (feet per year). For the transient model, the 

maximum evapotranspiration rate varied with stress 

periods ranging from between zero and 1 foot during 

the winter, and up to 4 feet during the summer. The 

extinction depth used in the simulation to represent the 

maximum rooting depth of plants was 5 feet, which 

is within the range of 2 to 5 feet used for various landcover types in the DPM (Boyd, 1996). 

Evapotranspiration from the water table as simulated by the EVT module only occurs where the water 

table is very close to land surface. This occurs in the 

upper Deschutes Basin in only a few settings. Very 

shallow water-table conditions occur over a wide area 

in the La Pine subbasin, and in the drainage of Indian 

Ford Creek. The water table can also be close to 

land surface in riparian zones along certain streams. 

Because the evapotranspiration process can lead to 

numerical instability, it is applied only to regions of 

the model where evapotranspiration from the water 

table is thought to represent a significant part of the 

hydrologic budget. Therefore, the evapotranspiration 

from the water table is only simulated in the Black 

Butte–Sisters area, alluvial areas northwest of Bend, 

and in the La Pine subbasin (fig. 5).

The final hydrologic-process boundary condition 

in this discussion is ground-water pumpage, which 

is simulated as a specified flux. Only irrigation and 

public-supply pumping is simulated in the model. 

Pumpage by individual domestic wells was not 

simulated because most of the pumpage is returned 

to the ground-water system through on-site sewagedisposal systems and the amount of consumptive use 

is negligible—less than 2 ft3/s. Specified fluxes are 

removed from cells corresponding to the geographic 

locations and production depths of wells. Pumpage 

from wells with open intervals that span multiple 

model layers was proportioned among the layers based 

on the percentage of the open interval in each layer. 

Pumping rates applied during the simulation vary with 

time. Techniques used to determine the locations and 

rates of pumping are described in Gannett and others 

(2001). The pumping rates for irrigation wells were 

determined using information from the OWRD waterrights files, satellite imagery, and data from State 

and County agricultural agencies. Pumping rates 

for public-supply wells were determined using data 

provided by the municipalities and water companies. 

The geographic distribution of irrigation and publicsupply pumpage in the model is shown in figure 8. 

Rates of pumping for each semiannual stress period 

are shown in figure 9. 

Model Parameters

Horizontal Hydraulic Conductivity

The scarcity of information on the subsurface 

geology in the saturated zone and of aquifer-test data 

precludes mapping hydraulic conductivity on a cellby-cell basis. A practical alternative is to represent 

hydraulic conductivity in a set of discrete subregions 

or zones within which conditions are considered 

uniform. The zones defined for the upper Deschutes 

Basin represent the hydraulic-conductivity distribution 

at a scale that can be supported with available geologic and hydrologic data. Lite and Gannett (2002) 

delineated seven hydrogeologic units in the upper 

Deschutes Basin that provide a general framework 

for the spatial distribution of hydraulic conductivity. 
22

Figure 8. Distribution of irrigation pumpage and public water-supply pumpage on the model grid. 
23

Figure 9. Rates of ground-water pumpage for irrigation and public-supply use for each semiannual stress period from 1978 to 1997. 

These hydrogeologic units are based on stratigraphy 

and facies distribution within the basin as determined 

from surficial geologic maps, well logs, analysis of 

well cuttings, and geophysical data. The 7 hydrogeologic units were ultimately further divided into 

17 zones based on several criteria, including the geographic density of hydraulic-head measurements, 

hydraulic gradients, recharge distribution, and analysis 

of residuals and sensitivity during calibration. Because 

the final zonation scheme was defined in part on 

hydraulic-head distribution, it is more detailed in the 

populated parts of the basin, where wells are plentiful. 

Hydraulic-conductivity zonation varies with depth 

in the model in places where geologic and hydrologic 

information are sufficient to support the level of detail. 

Model layer 1 includes a zone (6) to represent shallow 

alluvial deposits. Modifications to the hydraulic-conductivity zonation made during the model-calibration 

process are discussed in a later section. The final 

hydraulic-conductivity zonation for the model is 

shown in figure 10. 

Vertical Hydraulic Conductivity

The spatial distribution of vertical hydraulic 

conductivity follows the same zonation as horizontal 

hydraulic conductivity. Vertical hydraulic conductivity is uniform within each zone, but the value of 

vertical hydraulic conductivity varies between zones. 

In MODFLOW, the vertical hydraulic conductivities of two vertically adjacent cells are used 

to compute a vertical leakance term. The vertical 

leakance term (VCONT) is given by McDonald and 

Harbaugh (1988) as 

(3)

where 

∆zu is the thickness of the upper layer (L), 

∆zl is the thickness of the lower layer (L), 

Kz

u is the vertical hydraulic conductivity 

of the upper cell (LT-1), and 

Kz

l is the vertical hydraulic conductivity 

of the lower cell (LT-1). 

Storage Coefficients

The spatial distribution of storage coefficients 

follows the hydraulic-conductivity zonation scheme. 

Storage coefficients are constant only within each 

layer in each zone and can vary between layers within 

zones. Although all layers were modeled as confined, 

meaning transmissivity remains constant with respect 

to head fluctuations, storage coefficients in the uppermost layers were set to values typical of unconfined or 

semiconfined aquifers. This allowed the upper layers 

to behave unconfined or semiconfined with respect to 

their storage characteristics.

VCONT 1

z ∆ u ( ) ⁄ 2

Kzu

------------------- z ∆ l ( ) ⁄ 2

Kzl

+ ------------------

= ----------------------------------------------
24

Figure 10. Hydraulic-conductivity zonation in all model layers. (K, hydraulic conductivity, in feet per day)
25

Figure 10. Hydraulic-conductivity zonation in all model layers—Continued. (K, hydraulic conductivity, in feet per day)
26

MODEL CALIBRATION

Model calibration is the process whereby model 

parameter structure and parameter values are adjusted 

and refined to provide the best match between measured and simulated values of hydraulic heads and 

flows. Traditionally, models are calibrated by a 

trial-and-error process in which model parameters are 

adjusted within reasonable limits from one simulation 

to the next to achieve the best model fit. Model fit is 

commonly evaluated by visual comparison of simulated and measured heads and flows or by comparing 

root mean square (RMS) errors of heads and flows 

between simulations.

Models can also be calibrated using inverse 

methods in which the optimal parameter values for 

a given parameter structure are determined using a 

mathematical technique such as nonlinear regression 

(Cooley and Naff, 1990; Hill, 1992, 1998). This 

technique is sometimes referred to as parameter 

estimation. The regression process finds the set of 

parameter values that minimizes the difference 

between simulated and measured quantities such as 

hydraulic heads and flows.

The upper Deschutes Basin ground-water model 

was calibrated in two steps: (1) steady-state calibration 

to average conditions during the period 1993 to 1995, 

and (2) calibration to transient conditions from 1978 to 

1997. Inverse methods were used for the steady-state 

calibration, and traditional trial-and-error methods 

were used for the transient calibration. Inverse methods were not used for transient calibration because the 

combination of the sparse transient data set and relatively complicated model parameterization (developed 

using the rich steady-state data set) was not conducive 

to stable regression. The trial-and-error method provided satisfactory results for the transient calibration 

and met the overall goals of the project. The steadystate and transient calibration processes are discussed 

in the following paragraphs. 

Steady-State Calibration

The upper Deschutes Basin regional steady-state 

ground-water model was calibrated to average conditions from 1993 to 1995. This period was chosen 

because much of the data collection for the project 

occurred during this time and because average conditions during this period were reasonably close to longterm average conditions. For example, the average rate 

of recharge from precipitation during the 1993–95 

period is estimated to be about 3,500 ft3/s, while the 

1962–97 average is estimated to be about 3,800 ft3

/s 

(fig. 7). There is no evidence to suggest that the 

regional ground-water system in the upper Deschutes 

Basin is not in long-term equilibrium with the natural 

climate cycles and human activity. For example, 

no long-term water-level declines due to pumping 

were observed in the data, and (with few exceptions) 

ground-water discharge measurements show trends 

that can be related only to climate.

Steady-State Calibration Data 

Head observations for calibration consisted 

of water-level measurements from 983 wells. Head 

observations were not evenly distributed throughout 

the model domain but were clustered geographically 

in the populated areas of the study area (fig. 4) and 

vertically toward the upper parts of the flow system. 

Just over one-half of the head observations are in the 

upper two layers of the model, and 91 percent are in 

the upper five layers. Head observations were related 

to specific cells (horizontally and vertically) based on 

the well location and on the depth of the well’s open 

interval. Head observations from wells in which 

open intervals span multiple layers were compared 

to simulated heads in the associated layers according 

to methods described by Hill (1992). Although the 

steady-state model was calibrated to average conditions during the 1993 to 1995 time period (which 

approximates average long-term conditions), head 

measurements made outside of this time period were 

used during steady-state model calibration. This was 

necessary because obtaining water-level measurements during the calibration period was not possible 

in approximately one-half of the wells because of 

accessibility problems. In these cases, water levels 

measured by the driller or the USGS at an earlier 

time were used. This was acceptable hydrologically 

because time-series water-level data showed that 

water-level fluctuations in most of the model area 

were on the order of a few feet, rarely exceeding 

10 feet, and within the range of the measurement error 

from topographic maps (Caldwell and Truini, 1997). 

In addition, the temporal fluctuation is small in 

relation to the overall range in heads and large hydrologic gradients throughout most the modeled area. 

Head data were scrutinized carefully, and anomalous 

values due to measurement or location errors, pumping effects, or the obvious influence of immediately 

adjacent canals were removed from the calibration 

data set.
27

Calibration of this model benefited from a particularly rich set of flow observations. Flow observations were based on gain and loss estimates described 

in detail by Gannett and others (2001). Estimates of 

the average ground-water discharge to, or leakage 

from, 60 separate stream reaches were used in steadystate model calibration (fig. 11). Estimates along 

major streams, which account for the largest fluxes in 

the model, were based on gain-loss studies (seepage 

runs) conducted during the calibration period. Repeat 

measurements in subsequent years and comparison 

with gaging-station data indicated that the fluxes 

used for major streams were reasonably representative 

of average conditions. Estimates for some of the 

smaller streams and tributaries were based on gage 

data from both active and historic gages with periods 

of record ranging from a few years to decades. Some 

flux estimates in the upper Metolius River drainage 

were based, in part, on miscellaneous streamflow 

measurements. All flow observations were weighted 

to reflect their reliability.

Steady-State Calibration Procedure

The steady-state calibration was done using 

MODFLOWP, a version of MODFLOW modified 

by Hill (1992) to estimate parameters using nonlinear 

regression. In MODFLOWP, nonlinear least-squares 

regression is used to find the set of parameter values 

that will minimize the weighted sum of squared errors 

(SOSE) objective function (Hill, 1998):

(4)

where 

b is a vector containing values for each of the 

parameters being estimated, 

ND is the number of observations (hydraulic head 

and flow measurements), 

yi is the ith observation being matched by the 

regression, 

y’i(b)is the simulated value corresponding to the ith 

observation, and 

wi is the weight assigned to the ith observation. 

The differences between the observed and simulated values [yi - y’i (b)] are referred to as residuals. 

Residuals are weighted by the term wi for two reasons. 

First, weighting allows meaningful comparisons 

of measurements with different units because the 

weighted values are dimensionless. Secondly, weighting can be used to reduce the influence of measurements with large errors or uncertainty and increase 

the influence of measurements with small errors. 

In MODFLOWP, the weight for a particular observation is formulated as the inverse of the variance of the 

measurement error. 

The observation weights used in this study were 

developed in a manner similar to that suggested by 

Hill (1998). For hydraulic-head measurements, error 

was dominated by the accuracy of topographic maps 

used to determine land-surface elevation (and subsequently head). Other sources of error were considered negligible compared to map accuracy limitations. 

Map accuracy was taken to equal plus or minus one 

contour interval, which was usually either 10 or 

20 feet. A standard deviation of measurement error 

was estimated for each measurement by assuming that 

the 95-percent confidence interval for the true head 

was equal to the head measurement plus or minus the 

contour interval of the topographic map used to determine the well elevation.

Flow observations (estimates of the flux of water 

moving between the ground-water and surface-water 

systems) were also weighted by estimating standard 

deviations of measurement errors. Unlike head 

observations, however, there was no single dominant 

source of error. Flow observations assumed to represent average conditions were estimated from a variety 

of data sources and were generally known with much 

less certainty than head observations. The weighting 

of flow observations was based on subjective estimates of 95-percent confidence intervals developed 

considering the data and methods used to derive the 

flow observation. For the steady-state calibration, 

the 95-percent confidence interval associated with a 

particular observation was based on a subjective evaluation of how likely the observed value represented 

average 1993–95 conditions. In situations where the 

value was based on repeated measurements during 

the 1993 to 1995 period, confidence intervals were 

as small as plus or minus 10 percent. Confidence 

intervals were larger where values were inferred 

from gage data. Confidence intervals were largest, 

in some instances exceeding the flow estimate, where 

values were estimated from miscellaneous streamflow 

measurements or gage data outside the 1993 to 1995 

time period.

S b( ) wi yi y'

i [ ] – ( ) b 2

i = 1

ND

= ∑
28

Figure 11. Stream fluxes used for steady-state calibration. (Numbers correspond to observation numbers in table 3, p. 45–46.) 
29

Nonlinear regression was done using the 

MODFLOWP code of Hill (1992), which minimizes 

the objective function, equation (4), by nonlinear 

regression using a modified Gauss-Newton optimization method. The general basis of the technique is 

described by Cooley and Naff (1990) and the code 

and instructions are found in Hill (1992). Methods 

and guidelines for model calibration using inverse 

methods are presented by Hill (1998). 

Steady-State Model Parameters

During steady-state model calibration, the parameter structure allowed adjustment of 32 parameters. 

The boundary fluxes along the Cascade Range crest 

in the Metolius River drainage and east of Newberry 

Volcano, pumpage rates, and evapotranspiration 

parameters were not adjusted during calibration. These 

model parameters were not set up for optimization and 

not adjusted during calibration for a variety of reasons. 

First, they were derived independently and considered 

reasonably reliable. Boundary flux and evapotranspiration parameters were potential sources of problems 

with the regression due to correlation with other parameters and nonlinearity. The parameters adjusted 

during calibration include hydraulic conductivity in 

17 separate zones (fig. 10), streambed conductance 

values, recharge, modelwide vertical anisotropy, and 

additional vertical leakance parameters in 6 of the 

hydraulic-conductivity zones. The six vertical leakance parameters were set up for zones 3, 7, 8, 10, 11, 

and 16, where large vertical gradients existed and/or 

data were considered sufficient to support refinement 

of vertical hydraulic-conductivity estimates. The vertical leakance parameter allows adjustment of the vertical leakance (VCONT) between layers in subregions 

of the model in addition to the modelwide vertical 

anisotropy. This is done by adding another conductance term in the denominator of equation (3).

Numerous hydraulic-conductivity zonation 

schemes were tried during the calibration process. 

Zonation schemes were modified to reduce the sum of 

squared errors after optimization, reduce the number 

and size of geographic areas dominated by positive or 

negative residuals (model bias), minimize parameter 

correlation, and minimize the number of insensitive 

parameters. During this process, care was taken to 

remain faithful to the geologic data and overall geologic understanding.

Streams in the model were organized into 

seven groups, based on their geologic setting (table 1). 

Streambed conductance (CSTR) values within each of 

these groups could be adjusted up or down in unison.

Parameter Sensitivity

The ability to estimate values for particular 

parameters using nonlinear regression is affected by 

parameter sensitivity, which reflects the amount of 

information about a particular parameter provided by 

the available observation data. Changing values of 

parameters with low sensitivity has little effect on the 

sum of squared errors. Changing values of parameters 

with high sensitivities can have a relatively large 

effect on the sum of squared errors. 

MODFLOWP calculates the scaled sensitivity of 

a simulated value to a single parameter (Hill 1998) as 

(5)

where 

y’i is the simulated value associated with the ith 

observation, 

bj is the jth estimated parameter,

is the sensitivity of the simulated value 

associated with the ith observation with 

respect to the jth parameter and is evaluated 

at b (the set of parameter values at which the 

sensitivity is being evaluated), and 

wii is the weight of the ith observation. 

Table 1. Stream groupings used for streambed conductance 

parameters

Parameter 

number Streams included

KST1 Drainages in and on the margin of the 

Cascade Range south of Mt. Bachelor

KST2 Fall and Spring Rivers

KST3 Low gradient sections of the Deschutes and

Little Deschutes Rivers and tributaries 

in the La Pine subbasin

KST4 Deschutes River, Sunriver to Bend

KST5 Upper Tumalo and upper Squaw Creeks

KST6 Deschutes, Crooked, and Metolius Rivers and 

tributaries in the ground-water discharge area 

near and below Lake Billy Chinook 

KST7 Upper Metolius River and tributaries 

ssij bj ∂

∂y'

i

⎝ ⎠ ⎜ ⎟ ⎛ ⎞

bj

wii

1 2⁄ =

bj ∂

∂y'

i
30

The combined scaled sensitivity of all simulated 

equivalents of the observations to an individual 

parameter is provided by the composite scaled sensitivity (Hill, 1998). The composite scaled sensitivity is 

calculated for each parameter by MODFLOWP (Hill, 

1998) as

(6)

where 

ND is the number of observations, and 

ssij is the scaled sensitivity described above 

 (evaluated at b). 

Composite scaled sensitivities are dimensionless and their magnitudes are a measure of the total 

amount of information on a particular parameter provided by all the observations. Generally speaking, 

the regression method employed by MODFLOWP 

has difficulty estimating values for parameters with 

composite scaled sensitivities that are considerably 

lower than those for other parameters in the model. 

Models commonly include some parameters that have 

low composite scaled sensitivities and cannot be estimated using regression methods. These parameters 

are usually set to reasonable values and are not 

included in the regression process. 

cssj ssij ( )2

i = 1

ND

∑

b

⁄ ( ) ND

1 2⁄

=

The composite scaled sensitivities of all parameters in the upper Deschutes Basin model configured 

for optimization are shown graphically in figure 12. 

Parameters in figure 12 are grouped into those that 

were optimized and those that were not optimized. 

It can be seen that the nonoptimized parameters generally have low sensitivity. Low sensitivity parameters 

often destabilize the regression process, and can result 

in parameter value estimates that are physically unreasonable. Sensitivity is not the only factor, however, 

influencing the ability of the regression process to 

optimize certain parameters. It can be seen in figure 12 

that some parameters with low sensitivity were successfully optimized, and others with larger sensitivities (for example, recharge) were not. 

Figure 12. Composite scaled sensitivities of steady-state parameters. (Parameter types are as follows: K, horizontal hydraulic conductivity; 

KV, vertical leakance; ANIV, modelwide vertical anisotropy; RCH, recharge; KST, streambed conductance.) 

Final Parameter Values 

The final values for horizontal and vertical 

hydraulic-conductivity parameters in the steady-state 

calibration, both optimized and set, are shown in 

table 2. Parameter values not determined using 

regression were set using trial and error to values 

that were within reasonable limits prior to final optimization of other parameters. Generally, only small 

adjustments were made. Hydraulic conductivity in 

zone 6 (K6 on fig. 12), which represents Quaternary 

alluvial deposits along the margin of the Cascade 

Range from the Metolius River to Bend and in the 

La Pine subbasin, was set to 194 ft/d (2.25 × 10-3 ft/s). 
31

Table 2. Final hydraulic-conductivity values 

[Kxy, horizontal hydraulic conductivity; Kz, vertical hydraulic conductivity; 

parameters in bold type optimized]

Zone 

no.

Kxy

ft/s

Kxy

ft/day

Kz

ft/s 

Kz

ft/day Kxy / Kz

1 a2.41×10-4

a Parameter optimized.

a2.09×10+1 2.41×10-7 2.09×10-2 1,000

2 a6.48×10-5 a5.60 6.48×10-8 5.60×10-3 1,000

3 a3.54×10-4 a3.06×10+1 3.03×10-7 2.62×10-2 1,170

4 a8.99×10-3 a7.77×10+2 8.99×10-6 7.77×10-1 1,000

5 a2.32×10-4 a2.00×10+1 2.32×10-7 2.00×10-2 1,000

6 2.25×10-3 1.94×10+2 2.25×10-6 1.94×10-1 1,000

7 a7.00×10-4 a6.05×10+1 a1.66×10-8 a1.43×10-3 a42,200

8 a1.17×10-4 a1.01×10+1 a3.37×10-9 a2.91×10-4 a34,754

9 a3.60×10-4 a3.11×10+1 3.60×10-7 3.11×10-2 1,000

10 a1.14×10-2 a9.84×10+2 9.74×10-6 8.42×10-1 1,169

11 a2.26×10-4 a1.95×10+1 1.19×10-7 1.03×10-2 1,904

12 6.60×10-6 5.70×10-1 6.60×10-9 5.70×10-4 1,000

13 a1.69×10-4 a1.46×10+1 1.69×10-7 1.46×10-2 1,000

14 a1.07×10-3 a9.25×10+1 1.07×10-6 9.25×10-2 1,000

15 a2.79×10-3 a2.41×10+2 2.79×10-6 2.41×10-1 1,000

16 a3.68×10-3 a3.18×10+2 2.20×10-6 1.90×10-1 1,674

17 1.00×10-7 8.64×10-3 1.00×10-10 8.64×10-6 1,000

Hydraulic conductivity in zone 12 (K12 on fig. 12), 

which represents the central part of Newberry 

Volcano, was set to 0.57 ft/d (6.6 × 10-6 ft/s). This is 

an order of magnitude lower than the final value used 

for the central Cascade Range, but near the middle of 

the range given by Ingebritsen and others (1992) 

for near-surface rocks in the Cascade Range. The 

hydraulic conductivity of zone 17 (K17 on fig. 12), 

which represents a low-permeability barrier along 

the Green Ridge fault zone, was set to 8.6 × 10-3 ft/d 

(1.0 × 10-7 ft/s). 

The modelwide vertical anisotropy, the ratio of 

horizontal to vertical hydraulic conductivity (ANIV on 

fig. 12), was set to 1,000. The KV (vertical leakance) 

parameters act to further limit vertical conductance of 

certain zones. KV parameter values for zones 7 and 

8 were determined by regression, and for zones 3, 10, 

11, and 16 were set by trial and error.

Recharge rates were kept at the values determined by the DPM except in the southern part of the 

Cascade Range and the La Pine subbasin, where rates 

were reduced by 25 percent (a total of about 400 ft3/s) 

to better fit heads and flows. Recharge rates and 

hydraulic-conductivity values are often correlated 

parameters in many ground-water models, meaning 

that the observed head distribution can be simulated 

by a variety of values as long as the ratio of recharge 

and hydraulic conductivity remains the same. Such 

correlation was not generally an issue in the upper 

Deschutes Basin model because recharge rates were 

independently constrained by ground-water discharge 

measurements.

Streambed conductances were based on an 

assumed streambed thickness of 10 feet and streambed

hydraulic-conductivity values ranging from 0.1 to 

50 ft/d (with one reach set at 250 ft/d). Nearly all 

reaches had hydraulic conductivities less than 10 ft/d.

Streambed conductances were adjusted manually to 

achieve a good fit between simulated and observed 

stream gains and losses. Values were most commonly 

adjusted downward in losing reaches to better simulate 

the impediment of stream leakage by streambed materials and sediment-filled fractures. Stream reaches 

with large hydraulic-conductivity values were typically spring areas where ground water discharges 

from coarse sediments and open fractures in lava 

along streambanks and canyon walls.

All final parameter values (table 2) fall within or 

close to expected ranges, with the possible exception 

of the vertical hydraulic conductivity of zones 7 and 8. 

The regression process tended to set the KV parameter 

values for these two zones such that vertical anisotropies were on the order of 104, which is an order of 

magnitude higher than in other zones. Although 

anomalous, these values are not unreasonable. First, 

this area (particularly zone 7) is characterized by 

very high vertical head gradients, which suggests low 

vertical hydraulic conductivity. Anisotropies of 104

were necessary to simulate the observed gradients. 

Although high, these anisotropies are similar to or 

less than those found in some modeling studies in 

the Columbia River Basalt Group (Tanaka and others, 

1974; Mac Nish and Barker, 1976). It is unclear what 

geologic conditions may cause high vertical anisotropies in these zones; however, sparse subsurface data 

suggest this region of the model contains a higher 

proportion of layered lava sequences than elsewhere. 

The final value of hydraulic conductivity for zone 6 

(194 ft/d) is about twice the value of 88 ft/d estimated 

for the La Pine subbasin from aquifer tests (Century 

West Engineering Corporation, 1982). 
32

Steady-State Calibration Model Fit 

A number of statistical and graphical analyses of 

model fit, as described by Hill (1998), are presented in 

this section. Comparisons of simulated and measured 

heads and flows are included as well. Statistical and 

graphical techniques to evaluate model fit generally 

involve the analysis of residuals. A residual is the 

difference between a measured quantity, for example 

a head elevation, and its simulated equivalent [see 

equation (4)]. Of particular interest is the overall magnitude of the residuals and the degree to which they 

are random, independent, and normally distributed. 

Some of these characteristics are important in identifying model error and in evaluating confidence intervals 

on estimated parameter values. A number of statistical 

and graphical analyses of residuals are presented in 

this section along with more traditional comparisons 

of simulated and observed heads and flows using maps 

and graphs.

Statistical Measures of Model Fit 

and Parameter Uncertainty

The calculated error variance is an overall 

measure of the magnitude of residuals. The calculated 

error variance is simply the sum of squared errors 

from equation (4) divided by the number of observations minus the number of parameters (see equation 14 

of Hill, 1998). The square root of the calculated error 

variance is known as the standard error of regression.

If the model fit is consistent with the data accuracy, 

and the weighting appropriately reflects the data 

accuracy, both the calculated error variance and the 

standard error of regression should equal 1.0 (Hill, 

1998). Values larger than 1.0 suggest model error. 

In practice, the calculated error variance is almost 

always greater than 1.0. The calculated error variance 

and standard error of regression for the steady-state 

calibration described herein are 55.47 and 7.45, 

respectively. Sources of error are discussed later in 

the report. 

A more intuitive measure of overall model fit is 

the product of the standard error of regression and the 

statistics used to determine weights (usually standard 

deviations or coefficients of variation). These values 

are referred to by Hill and others (1998) as fitted standard deviations or fitted coefficients of variation. The 

weights used for head observations in the steady-state 

model calibration are based on the assumed error in 

topographic maps used to determine well elevations. 

The standard deviation of measurement error for 

most wells was 10.2 feet. Multiplying this figure by 

the standard error of regression yields a fitted standard 

deviation for heads of 76 feet.

It is desirable for weighted residuals to be 

independent, random, and normally distributed. 

Weighted residuals that are not so may be an indication of model error. Moreover, calculation of certain 

statistics, such as confidence intervals on parameter 

values, requires residuals to have these properties. 

Independence and randomness of weighted residuals 

can be assessed by plotting them against weighted 

simulated values (Cooley and Naff, 1990, p. 170; Hill, 

1994, p. 3). Residuals from the steady-state calibration 

are plotted in figure 13. If weighted residuals are 

random, they should be equally distributed above and 

below the horizontal zero axis. In addition, weighted 

residuals should show no trends as a function of 

weighted simulated value; in other words, there should 

be no distinct slope or curve to the data. Weighted 

residuals plotted in figure 13 generally meet both these 

criteria with the exception of a cluster of points on the 

right-hand side of the graph with weighted residuals 

that are all negative. These correspond to head observations in the La Pine subbasin, primarily in zone 5. 

The predominantly negative residuals indicated that 

simulated heads are systematically too high in this 

area. Various parameterization schemes were tried 

and none satisfactorily resolved the problem, but 

usually resulted in problems matching stream fluxes 

or heads elsewhere. Efforts to devise a solution were 

hampered by the lack of deep subsurface geologic 

information, particularly in the area just north of 

the La Pine subbasin (zone 3). The high simulated 

heads in the La Pine subbasin do not significantly 

affect overall simulation results, but may be problematic for detailed analysis in that area. 

Certain inferential statistics, such as confidence 

intervals on parameter estimates, require that residuals 

be normally distributed. The degree to which weighted 

residuals are normally distributed can be assessed by a 

normal probability plot (Draper and Smith, 1998). If 

ordered residuals are normally distributed, they will 

fall approximately on a straight line on such a plot. 

Figure 14 is a normal probability plot of head and flow 

residuals for the steady-state model. The plot shows 

both head and flow residuals are generally linear but 

that there is slight curvature, particularly at the tails of 

the plot. There are, however, very few outliers.
33

Figure 13. Weighted residuals as a function of weighted simulated values.

Figure 14. Normal probability plot of ordered weighted residuals.

An associated statistic for testing the normality 

of the distribution of residuals is the correlation 

coefficient between the ordered weighted residuals 

and the order statistics from a normal probability 

distribution function, denoted as R2

N (Hill, 1992, 

1998). Values of R2

N close to 1.0 indicate that 

weighted residuals are independent and normally 

distributed. If R2

N is significantly less than 1.0, the 

weighted residuals are not likely independent and 

normally distributed. Critical values of R2

N below 

which residuals are not likely independent and 

normally distributed (at 0.05 and 0.1 significance 

levels) are given by Hill (1998, appendix D). 

The critical values are a function of the number of 

observations, with the values increasing with the 

number of observations. The largest number of 

observations for which Hill (1998) gives critical 

values of R2

N is 200, and the steady-state model 

has 1,043 observations. The critical R2

N for 

200 observations is 0.987 so the value for 1,043 

observations would be slightly higher but still less 

than 1.0. The R2

N value for the steady-state model 

calibration is 0.977, slightly below the critical value 

for 200 observations at the 0.05 significance level. 
34

By this measure, the weighted residuals are not likely 

to be truly independent and normally distributed. 

Numerical models of complex ground-water 

systems are usually nonlinear. The degree of nonlinearity can be evaluated by the modified Beale’s 

measure (Cooley and Naff, 1990). The modified 

Beale’s measure is calculated from MODFLOWP 

output using the program BEALEP (Hill, 1994). 

The program also calculates critical values for evaluating the modified Beale’s measure. These critical 

values indicate that the upper Deschutes Basin model 

would be considered effectively linear if the modified 

Beale’s measure is less than 0.05 and would be considered highly nonlinear if the modified Beale’s measure 

is greater than 0.6. The modified Beale’s measure for 

the calibrated steady-state model is 0.3, indicating 

that it is moderately nonlinear. Principal causes of 

nonlinearity are head-dependent boundaries such as 

streams and evapotranspiration processes. 

In order to provide information on parameter 

uncertainty, MODFLOWP calculates 95-percent 

linear confidence intervals for parameter estimates 

(fig. 15). Confidence intervals represent the precision 

with which the parameters are estimated given the 

available data. Parameter confidence intervals from 

a model can be compared to the range of reasonable 

parameter values, determined independently of the 

model using aquifer tests, specific-capacity tests, 

and the literature, to assess the reasonableness of the 

estimated values. The expected ranges for values in 

this model are described in the preceding discussion 

of hydraulic conductivity and by Gannett and others 

(2001). The relatively small range of the confidence 

intervals on figure 15 indicates that the available 

data provide sufficient information to narrowly constrain the parameter estimates. The closeness of the 

estimated parameter values and confidence intervals 

to the expected values indicates that the estimates 

are consistent with independent hydrogeologic information. The confidence intervals for the hydraulicconductivity parameters for zones 4 and 5 are larger 

than for other hydraulic-conductivity parameters likely 

because of sparse or poorly distributed observations 

in those zones. These larger confidence intervals are 

consistent with the low composite scales sensitivities 

for these parameters on figure 12. 

Figure 15. The final values, 95-percent confidence intervals, and ranges of expected values for optimized steady-state hydraulic-conductivity 

parameters. 
35

The method used to calculate the 95-percent 

confidence intervals assumes that the model is linear 

in the vicinity of the optimized parameter values 

and that residuals are independent and normally 

distributed. As previously discussed, the requirements 

of model linearity and independent normally distributed residuals are not strictly met; however, measures 

indicate that the model is not highly nonlinear and that 

residuals are reasonably close to being independent 

and normally distributed. Because the assumptions 

are not strictly met, the computed linear confidence 

intervals shown in figure 15 must be considered 

approximate indicators of parameter uncertainty and 

used with caution. 

Parameter correlation can be a problem in 

ground-water models. Two parameters can be highly 

correlated when available information is insufficient 

to determine unique values for each. Coordinated 

changes in correlated parameter values will result in 

the same heads and flows. MODFLOWP calculates 

a matrix of parameter correlation coefficients using 

the methods described by Hill (1998). Parameter 

correlation coefficients range from -1.0 to +1.0. Correlation coefficients with absolute values close to 1.0 

indicate a high degree of correlation. In such cases the 

effects of parameter correlation need to be evaluated. 

Correlation coefficients with absolute values less than 

about 0.95 are considered acceptable (Hill, 1998). All 

correlation coefficients for the steady-state calibration 

were well within the acceptable range, with the largest 

being 0.86. 

Comparison of Simulated and Observed 

Steady-State Heads

Model calibration also can be assessed qualitatively by visually comparing simulated hydraulic 

heads and flows with measured values. Figure 16 

shows the simulated head distribution, representing 

mean heads for the 1993–95 period, contoured by 

computer (using the program of Harbaugh, 1990) 

as well as contours drawn by hand from field measurements for each model layer. Hand-drawn head 

contours are based primarily on field measurements 

from wells and from elevations of major springs and 

gaining streams. At certain locations in some layers, 

head contours were constrained by field measurements 

in overlying and underlying layers. The variation in 

hydraulic head with depth was an issue when hand 

contouring measurements from wells. In areas with 

high vertical gradients, the measured head at any 

particular location can vary depending on the vertical 

location within the layer. Moreover, the open intervals 

of some wells span multiple layers, and the measured 

head is a function of the heads in all the water-bearing 

zones penetrated. In general, data for each layer 

(including the multilayer measurements) were fairly 

consistent, and hand-drawn contours are well constrained where there are abundant data points. The 

placement of contours is much less certain where they 

are constrained by only one or two data points or 

where data are widely spaced. In figure 16, contours 

based on field measurements are drawn only where 

sufficient data were available. The hand-contoured 

area and number of contours present diminished with 

depth as data become sparse. Hand-drawn contours are 

solid where well constrained and dashed where data 

are sparse or ambiguous. Hand-drawn head contours 

closely follow or overlie simulated head contours 

where the simulated contours are entirely consistent 

with field data.

Contours of the simulated head distribution are 

generally consistent with contours hand drawn based 

on field measurements. The direction and magnitude 

(slope) of simulated head gradients match measured 

regional gradients well. There are some areas, 

however, where simulated head contours consistently 

differ from those drawn from field measurements. 

Simulated heads are generally higher than measured 

heads in the La Pine subbasin, as described earlier in 

the discussion of figure 13. This is most noticeable 

in layer 1, but occurs in lower model layers as well. 

Another area where simulated heads appear to be consistently different from measured heads is in layers 2 

through 5 north of Cline Buttes between Redmond 

and Sisters. The simulated 2,600 foot contour in these 

layers generally forms an east-west arc gently convex 

toward the south. Field measurements suggest that 

the 2,600 foot contour is more curved toward the 

south, becoming v-shaped with the apex close to Cline 

Buttes.

In areas where head data are insufficient to draw 

contours, a direct comparison between simulated 

and actual head distributions is not possible. In such 

areas, for example the Cascade Range and Newberry 

Volcano, simulated heads are generally consistent 

with what would be expected given the topography, 

precipitation distribution, and inference from streams 

and sparse well data. 
36

Figure 16. Lines of equal measured head and lines of equal simulated head for all model layers. 
37

Figure 16. Lines of equal measured head and lines of equal simulated head for all model layers—Continued. 
38

Figure 16. Lines of equal measured head and lines of equal simulated head for all model layers—Continued. 
39

Figure 16. Lines of equal measured head and lines of equal simulated head for all model layers—Continued. 
40

Figure 16. Lines of equal measured head and lines of equal simulated head for all model layers—Continued. 
41

Figure 16. Lines of equal measured head and lines of equal simulated head for all model layers—Continued.
42

Figure 16. Lines of equal measured head and lines of equal simulated head for all model layers—Continued. 
43

Figure 16. Lines of equal measured head and lines of equal simulated head for all model layers—Continued. 
44

Comparison of Simulated and Observed 

Steady-State Ground-Water Flow 

to and from Streams 

The reasonableness of simulated groundwater discharge to, or recharge from, streams can 

be assessed by comparing simulated values with estimates derived by independent methods. Estimates 

of the average rates and distribution of ground-water 

flow to and from the Deschutes, Crooked, and Metolius Rivers and most major tributaries in the upper 

Deschutes Basin were developed by Gannett and 

others (2001) based on a variety of data from seepage 

runs conducted by OWRD, gaging stations, and 

miscellaneous streamflow measurements. Groundwater flow estimates to or from 60 stream reaches 

were used in the steady-state model calibration. The 

locations of the 60 reaches are shown in figure 11. 

Estimated and simulated flows are given in table 3, 

and a comparison of estimated and simulated flows 

is shown in figure 17. These data are also shown as 

a scatter plot in figure 18. 

The pattern of stream gains and losses represented by the 60 steady-state flow estimates 

contains spatial variability that cannot be reasonably 

simulated with a regional-scale model. The distribution of stream gains and losses is controlled by 

the interplay of geology, topography, streambed 

characteristics, and the hydraulic-head distribution in 

the ground-water system. Data describing the spatial 

variability of these factors (particularly subsurface 

geology) is not available in most parts of the upper 

Deschutes Basin (for example the Cascade Range) 

at the scales comparable to the resolution of the gain 

and loss estimates. In addition, topography, a critical 

factor in the interaction between ground water and 

streams, is greatly generalized in the regional-scale 

model. 

The regional-scale model, therefore, does not 

include the level of complexity required to simulate 

the spatial variations in stream gains and losses 

represented by the available flow estimates. It is 

useful, therefore, to aggregate stream reaches and 

consider net flows of water moving between the 

ground-water and surface-water systems in subregional areas at scales that reflect model complexity 

and represent overall hydrologic settings. 

Comparisons of estimated and simulated 

ground-water flows to and from streams aggregated 

in 13 hydrologic settings (table 4) show that the 

steady-state model, in general, reasonably simulates 

the observed magnitude and direction of groundwater flow to and from streams at regional and subregional scales. Estimated and simulated values 

are shown graphically in figure 19. Simulated 

ground-water flows to streams are within a factor 

of 1.4 of estimated values in most major discharge 

areas.

There are two settings where simulated flows 

are off by a factor of nearly 20. The largest differences between simulated and estimated flows to 

and from streams occur along the Little Deschutes 

River downstream of Crescent Creek in the flat 

part of the La Pine subbasin, where the simulated 

ground-water discharge of 88 ft3

/s greatly exceeds 

the estimated value of 5 ft3/s. This is expected given 

that, as previously described, simulated heads are 

consistently higher than observed heads in that area. 

The excess discharge here is compensated for in the 

overall budget by simulated discharge to streams that 

is slightly lower than observed in the Cascade Range 

to the west, and simulated losses slightly larger than 

expected along the Deschutes River to the north. 

Flow residuals are also large in the area drained 

by tributaries to the upper Metolius River, where 

the simulated ground-water discharge of 13 ft3/s 

is only about 5 percent of the observed groundwater discharge of 260 ft3/s. Simulated discharge 

to the Metolius River main stem is reasonably close to 

the estimated value. 

The rate and distribution of ground-water 

discharge to major inflow areas along the main 

stems of the Deschutes, Crooked, and Metolius 

Rivers, including the area near their confluence, 

are, in general, close to the measured or estimated 

values. Critical stream reaches along which simulated 

ground-water discharge rates are close to estimated 

rates include the main stem of the Metolius River 

from the headwaters to Jefferson Creek, the Deschutes River between Lower Bridge and the gage 

near Culver, the Crooked River from Osborne 

Canyon to the gage below Opal Springs, and Lake 

Billy Chinook downstream to Lake Simtustus. A 

significant part of the total ground-water discharge 

in the upper Deschutes Basin occurs along these 

reaches, so the close fit between simulated and 

estimated rates of ground-water discharge is noteworthy.
45

Table 3. Measured or estimated steady-state stream gains and losses due to ground-water exchange, and simulated equivalents 

[Observation locations are shown by number on figure 11; RM, river mile; all values are cubic feet per second; Res., reservoir]

 Estimated 95% confidence 

interval on measured or 

estimated flux

Observation 

number Stream reach

Measured or 

estimated 

stream flux

Estimated 

2 sigma Lower Upper 

Simulated 

stream flux

1 Little Deschutes River—above RM 86 -19.7 10 -29.7 -9.7 0.0

2 Little Deschutes River—RM 86–80 -11.8 10 -21.8 -1.8 0.0

3 Little Deschutes River—RM 80–71 2.6 10 -7.4 12.6 -6.6

4 Little Deschutes River—RM 71–64.5 2.0 10 -8.0 12.0 -4.0

5 Big Marsh Creek -21.0 10 -31.0 -11.0 -47.5

6 Crescent Creek—above RM 18.5 -23.7 10 -33.7 -13.7 -12.0

7 Crescent Creek—RM 18.5–2.2 1.5 10 -8.5 11.5 -24.3

8 Little Deschutes River—RM 57.5–43 -2.3 10 -12.3 7.7 -32.1

9 Little Deschutes River—RM 43–33.5 -2.9 10 -12.9 7.1 -25.7

10 Little Deschutes River—RM 33.5–26.7 0.0 1 -1.0 1.0 -7.9

11 Little Deschutes River—RM 26.7–13.7 0.6 10 -9.4 10.6 -22.4

12 Paulina Creek—RM 13–5.2 6.1 1 5.1 7.1 12.7

13 Odell Lake -41.0 20 -61.0 -21.0 -23.6

14 Odell Creek—Lake outlet to gage -41.0 10 -51.0 -31.0 -21.1

15 Odell Creek and Davis Creek above RM 3 -191.0 25 -216.0 -166.0 -102.0

16 Davis Creek below RM 4 (west arm Wickiup Res.) -15.0 5 -20.0 -10.0 -7.9

17 Cultus River above Crane Prairie Res. -63.0 20 -83.0 -43.0 -7.9

18 Deschutes River above Crane Prairie Res. -151.0 20 -171.0 -131.0 -14.9

19 Deschutes River between Crane Prairie Res. and Wickiup Res. -229.0 50 -279.0 -179.0 -146.0

20 Wickiup Res. 0.0 50 -50.0 50.0 -20.1

21 Deschutes River—Wickiup Res. to RM 217.6 1.5 10 -8.5 11.5 -7.2

22 Deschutes River—RM 217.6–208.6 -6.0 10 -16.0 4.0 -6.4

23 Deschutes River—RM 208.6–199.7 -24.0 10 -34.0 -14.0 -46.2

24 Fall River above RM 5 -119.0 20 -139.0 -99.0 -140.0

25 Fall River—RM 5–0.4 7.0 5 2.0 12.0 -34.0

26 Deschutes River—RM 191.5–189 (includes Spring River) -271.0 50 -321.0 -221.0 -63.8

27 Deschutes River—RM 198–181.4 23.8 10 13.8 33.8 28.4

28 Deschutes River—RM 181.4–164.4 (Benham Falls to Bend) 89.0 10 79.0 99.0 157.0

29 Tumalo Creek above RM 11 (Columbia Southern diversion) -68.0 20 -88.0 -48.0 0.0

30 Tumalo Creek—RM 11–3.2 -1.0 5 -6.0 4.0 0.0

31 Tumalo Creek—RM 2.3–0.3 0.2 1 -0.8 1.2 0.0

32 Deschutes River - RM 164.4–145.3 (Bend to Cline Falls) -15.0 3 -18.0 -12.0 -10.4

33 Deschutes River—RM 138–134 (near Lower Bridge) -10.0 2 -12.0 -8.0 -22.5

34 Deschutes River—RM 134–124.9 (below Lower Bridge) -170.0 20 -190.0 -150.0 -210.0

35 Deschutes River—RM 124.9–123.3 -102.0 10 -112.0 -92.0 -31.1

36 Squaw Creek above RM 26.8 (gage) -65.0 10 -75.0 -55.0 -109.0

37 Indian Ford Creek above RM 10.7 (Black Butte Ranch springs) -6.1 3 -9.1 -3.1 0.0

38 Indian Ford Creek—RM 10.7–8.0 1.7 2 -0.3 3.7 0.0

39 Indian Ford Creek—RM 8.0–3.6 0.6 1 -0.4 1.6 0.0

40 Indian Ford Creek—RM 3.6–2.1 0.6 1 -0.4 1.6 0.0

41 Indian Ford Creek—RM 2.1–1.3 2.8 1 1.8 3.8 0.0

42 Indian Ford Creek—RM 1.3-–0.8 0.5 1 -0.5 1.5 0.0

43 Indian Ford Creek—RM 0.8 to mouth 0.0 1 -1.0 1.0 0.0

44 Squaw Creek—RM 19.0–16.5 -6.6 1 -7.6 -5.6 0.0
46

45 Squaw Creek—RM 1.7 to mouth (includes Alder Springs) -94.0 10 -104.0 -84.0 -124.0

46 Deschutes River—RM 123.0–120.0 (Culver gage) -25.0 5 -30.0 -20.0 -88.8

47 Crooked River—RM 20.4–13.8 -70.0 15 -85.0 -55.0 -204.0

48 Crooked River—RM 13.8–6.7 (Opal Springs gage) -1,006.0 200 -1,206.0 -806.0 -990.0

49 Lake Billy Chinook -420.0 100 -520.0 -320.0 -558.0

50 Deschutes River Round Butte Dam to Lake Simtustus -200.0 75 -275.0 -125.0 -22.4

51 Lake Simtustus -51.0 15 -66.0 -36.0 -14.4

52 Metolius River—RM 41.2 (headwaters) to 38.1 -257.0 50 -307.0 -207.0 -244.0

53 Suttle Lake -31.0 10 -41.0 -21.0 0.0

54 Lake Creek—Suttle Lake outlet to mouth -36.0 10 -46.0 -26.0 -5.0

55 First Creek—upper part of drainage -1.5 5 -6.5 3.5 0.0

56 Jack Creek—upper part of drainage -46.0 10 -56.0 -36.0 0.0

57 Canyon Creek—entire drainage above RM 0.5 -60.0 10 -70.0 -50.0 -6.2

58 Abbot Creek—entire drainage above RM 1.5 -12.0 5 -17.0 -7.0 -2.1

59 Candle Creek—entire drainage above RM 2.0 -73.0 25 -98.0 -48.0 -0.1

60 Metolius River—RM 38.1–28.8 (Jefferson Creek) -724.0 100 -824.0 -624.0 -512.0

Table 3. Measured or estimated steady-state stream gains and losses due to ground-water exchange, and simulated equivalents 

—Continued

[Observation locations are shown by number on figure 11; RM, river mile; all values are cubic feet per second; Res., reservoir]

 Estimated 95% confidence 

interval on measured or 

estimated flux

Observation 

number Stream reach

Measured or 

estimated 

stream flux

Estimated 

2 sigma Lower Upper 

Simulated 

stream flux

Figure 17. Measured or estimated steady-state stream fluxes, 95-percent confidence intervals on measured or estimated stream fluxes, 

and simulated steady-state stream fluxes. (The left-hand axis corresponds to observations left of the bold vertical line, and the 

right-hand axis corresponds to observations with larger values to the right of the bold vertical line.) 
47

Figure 18. Simulated stream fluxes as a function of measured or 

estimated stream fluxes. 

Table 4. Measured or estimated gains and losses to aggregated stream reaches, and simulated equivalents

[All values are cubic feet per second; RM, river mile]

 95% Confidence 

interval on measured or 

estimated flux

Stream reaches

Figure 19 

designation

Measured or 

estimated 

stream flux Lower Upper 

Simulated 

stream flux

Little Deschutes River and tributaries above Crescent Creek LDAC -70 -140 0 -94

Little Deschutes River below Crescent Creek LDBC -5 -36 26 -88

Odell and Davis Creek drainages ODC -288 -348 -228 -155

Wickiup Reservoir and tributaries from the north WR -443 -583 -303 -189

Deschutes River and tributaries from the west between 

Wickiup Reservoir and Sunriver

DWS -412 -517 -307 -298

Upper Tumalo and upper Squaw Creeks UTS -133 -163 -103 -109

Deschutes River between Sunriver and Bend DSB 113 93 133 185

Deschutes River between Bend and Lower Bridge DBLB -25 -30 -20 -33

Deschutes River between Lower Bridge and the Culver gage DLBC -391 -436 -346 -454

Crooked River between Osborne Canyon and the Opal Springs gage COO -1,076 -1,291 -861 -1,194

Lake Billy Chinook and downstream down to and including 

Lake Simtustus

LBC -671 -861 -481 -595

Metolius River main stem above RM 28.8 MM -981 -1,131 -831 -756

Metolius River tributaries above Jefferson Creek MT -260 -335 -185 -13

Figure 19. Measured or estimated steady-state stream fluxes, 

95-percent confidence intervals on measured or estimated 

stream fluxes, and simulated steady-state stream fluxes for 

aggregated stream reaches. (Stream-reach abbreviations 

explained in table 4.) 
48

Steady-State Model Water Budget

The steady-state model water budget is evaluated 

by looking at the instantaneous rates of the various 

budget components (table 5). The water balance in 

the model closes to within about 0.02 percent. The 

recharge, boundary flux, and pumpage components 

of the budget are specified. Stream leakage and 

evapotranspiration are head-dependent quantities and 

are calculated by the model. Total stream-leakage 

figures in table 5 do not match the total measured 

(or estimated) stream fluxes given in tables 3 and 4 

because not all stream cells are represented in tables 3 

and 4. Moreover, the gain and loss observations in 

tables 3 and 4 are net values integrated over many 

cells, whereas the budget values are the gross sums of 

all individual cells. The comparison of measured and 

simulated stream fluxes presented earlier indicates that 

simulated stream leakage is reasonable.

It is more difficult to evaluate the simulated 

evapotranspiration rate of 420 ft3/s. Gannett and 

others (2001) estimated that evapotranspiration from 

the water table in the La Pine subbasin was probably 

less than 80 ft3/s. The simulated value of 420 ft3/s 

includes evapotranspiration from a much larger 

area that includes the upper Indian Ford Creek and 

Metolius River drainages, where shallow ground water 

is common. Because heads are systematically too high 

in the La Pine subbasin, as previously described, it can 

be assumed that evapotranspiration is overestimated in 

that area.

Transient Calibration

The upper Deschutes Basin regional groundwater model was also calibrated to transient conditions for the 20-year period between 1978 and 1997. 

This period was chosen for several reasons. First, 

it was a period of rapid growth in population in the 

basin. It is also a period of sufficient length to include 

decadal-scale climate variations that are known to be 

the principal force driving head and discharge 

fluctuations. Lastly, this period includes time-series 

water-level measurements collected by the USGS 

in the late 1970s as part of an earlier project.

The initial conditions for the transient simulation 

were based on climate conditions for the 4 years prior 

to the calibration period. The first stress period for the 

transient calibration is the winter of 1978 (November 

1977 through April 1978). Conditions during the 4 

years prior to the beginning of the transient calibration 

period were slightly wetter than the average conditions 

during the 1993–95 water years used for steady-state 

calibration by a factor of about 1.15. The starting 

heads for the transient calibration were derived from 

a steady-state simulation using the 1993–95 average 

recharge rate multiplied by 1.15.

Transient Calibration Data

Time-series head data and ground-water discharge data reflecting the dynamic behavior of the 

regional ground-water system were used for transient 

calibration. During the calibration period, head and 

discharge exhibited decadal and seasonal variations in 

recharge from precipitation, and seasonal variations in 

recharge from canal operation. The driving boundary 

conditions have been described in previous sections.

There are approximately 103 wells with timeseries water-level data in the upper Deschutes Basin 

study area (Caldwell and Truini, 1997). The frequency 

of measurements varies from 2 hours to once or twice 

per year, with some wells having multiple-year gaps in 

records. The duration of record for wells varies from 

less than a year to over 50 years. Of available wells 

with time-series data, 64 were evaluated during transient calibration. A number of factors were considered 

while selecting wells for use in transient calibration. 

The overall goal was to have good geographic 

coverage, all hydrologic settings represented, and, 

where possible, consistent records throughout the 

calibration period. Wells were generally not used if 

they were monitored for only a short period of time or 

included too many nonstatic (pumping or recovering) 

water levels. Wells that appeared to be responding 

to influences that cannot be represented in a regionalscale model, such as large influences of nearby canals 

or pumping wells, were not included. 

Table 5. Steady-state model simulated water budget 

[All values in cubic feet per second; ---, not applicable]

Budget component In Out

Recharge (precipitation) 3,106 ---

Recharge (canal) 492 ---

Recharge (on-farm) 68 ---

Recharge (total) 3,666 ---

Subsurface inflow 847 ---

Stream leakage 403 4,464

Evapotranspiration --- 420

Pumpage --- 31

Total 4,916 4,915
49

Information on temporal variations in groundwater discharge to streams (stream fluxes) was derived 

solely from stream-gage data. Stream-gage data can be 

used to estimate temporal variations in ground-water 

discharge only where certain conditions exist. Where a 

gage is on a stream that is entirely spring fed, the gage 

provides a direct record of ground-water discharge. 

For some streams, late season flows can be considered 

to be largely or entirely ground-water discharge. In 

these situations, the gage record can be used to 

estimate year-to-year variations but are not useful for 

evaluating seasonal variations. Where two gages exist 

on a stream, the ground-water inflow to the intervening reach can be estimated by comparing the gages, 

provided tributary inflows and diversions can be 

quantified. All of these situations occur in the upper 

Deschutes Basin.

The maximum ET rate was varied each stress 

period according to residual ET calculated by the 

DPM, as explained in the preceding section on boundary conditions. Public water-supply pumping was varied each stress period, with values based on data from 

water providers and population records. Irrigation 

pumping for the transient calibration was based on 

the 1993–95 steady-state values but varied according 

to climate-driven variations in consumptive use. 

Pumping for individual irrigation wells was included 

only for stress periods postdating the issuance of the 

associated water right. 

Transient Calibration Procedure

Transient calibration was achieved using traditional trial-and-error methods. Regression methods 

were not employed because the model parameterization scheme, which was developed using the steadystate data set, was not well suited for allowing meaningful, stable, nonlinear regression using the transient 

data set. Transient data were much fewer than steadystate data and had a different spatial distribution. 

Trial-and-error methods for the transient calibration 

provided satisfactory results and met the overall 

project goals.

The same model parameter structure was used 

for the transient calibration as was used for the 

steady-state calibration. The only difference was 

the addition of storage coefficients. The storage 

coefficient parameterization followed the hydraulicconductivity zonation. After experimentation it was 

determined that satisfactory transient calibration 

could be achieved using all of the parameter values 

from the steady-state calibration and adjusting only 

storage coefficients. Storage coefficients were adjusted by layer in each of the hydraulic-conductivity 

zones. In general, best results were achieved when 

storage coefficients decreased with depth (table 6). 

Table 6. Final storage coefficient values

[na, zone not present in the layer]

Storage coefficient

Zone Layer 1 Layer 2 Layer 3 Layer 4 Layer 5 Layer 6 Layer 7 Layer 8

1 0.004 0.002 0.002 0.00004 0.00004 0.00004 0.00004 0.000004

2 0.008 0.002 0.002 0.00004 0.00004 0.00004 0.00004 0.000004

3 0.15 0.0015 0.0015 0.00015 0.00015 0.00015 0.00015 0.00015

4 0.01 na na na na na na na

5 na 0.01 0.01 0.01 0.01 0.01 0.01 0.001 a

6

a Metolius River drainage. 

0.3 na na na na na na na b

6

b All other areas. 

0.03 na na na na na na na

7 0.04 0.02 0.002 0.0006 0.0002 0.0002 0.00005 0.00005

8 0.01 0.0005 0.0005 0.0005 0.0005 0.0005 0.0001 0.0001

9 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.0001

10 0.05 0.05 0.05 0.05 0.05 0.01 0.0005 0.0005

11 0.02 0.004 0.004 0.004 0.002 0.002 0.0004 0.0004

12 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.0001

13 0.08 0.08 0.08 0.08 0.08 0.04 0.008 0.0008

14 0.05 0.05 0.05 0.01 0.01 0.005 0.001 0.0001

15 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.0001

16 0.002 0.001 0.001 0.001 0.001 0.001 0.001 0.0001

17 0.01 0.001 0.001 0.001 0.001 0.001 0.001 0.0001
50

Storage coefficients ranged from 0.3, for coarse, 

near-surface alluvial deposits, to 4 × 10-6, for materials 

greater than 1,000 feet below the top of the saturated 

zone in the Cascade Range (zones 1 and 2). Specific 

storage values ranged from approximately 10-3 to 10-8

and also decreased with depth. This change in storage 

coefficients with depth is to be expected in volcanic 

terranes where the system is unconfined or semiconfined near the surface and becomes increasingly 

confined with depth. Moreover, specific storage 

values can be expected to decrease with depth as 

well, as rocks become increasingly altered and have 

less porosity. The final storage coefficient values are 

within reasonable and expected ranges, and are consistent with aquifer-testing results described in preceding 

sections (also see Gannett and others, 2001, table 1). 

The spatial distribution of final storage coefficient 

values is consistent with geologic information, with 

larger values generally occurring in parts of the Deschutes Formation where volcaniclastic and alluvial 

deposits are prominent in the section. Lower values 

were generally related to regions dominated by finegrained sediments or thick lava sequences without 

coarse interbeds.

Comparison of Simulated and Observed Transient Heads

Time-series water-level data from 103 wells 

were available for the modeled area. As previously 

discussed, not all of these wells were suitable for use 

in model calibration because water levels in some 

wells responded to stresses that cannot be simulated in 

a regional-scale model. Water levels fluctuate in the 

upper Deschutes Basin in response to a variety of driving forces that work at different temporal and spatial 

scales. Long-term climate cycles are the dominant 

driving force, affecting the entire region and working 

on a decadal time scale. Annual recharge cycles also 

work over the entire region. Canal leakage and irrigation, also important driving forces, affect only parts of 

the model area and have an annual cycle that is out of 

phase with the timing of natural recharge. Which of 

these two seasonal processes dominates is largely a 

function of location. Other forces include stream-stage 

variations and ground-water pumping. These work at a 

variety of scales, generally small relative to other 

stresses. The model simulates the fluctuations caused 

by climate cycles, natural recharge from precipitation, 

and canal leakage and irrigation. Stream-stage 

variation is not explicitly represented in the model, 

so water-level fluctuations resulting from variations 

in stream stage are not simulated. Although groundwater pumping is included in the model, drawdown 

effects are averaged over entire model cells, so 

large fluctuations close to pumping wells are not 

simulated.

As was discussed in the previous section on 

steady-state calibration, observed water-level elevations are not matched exactly, particularly in regionalscale models with large head gradients. In the upper 

Deschutes Basin, the water-table elevation ranges over 

at least 4,500 feet. Elevations of head measurements 

from wells used in the model calibrations range from 

1,797 to 5,586 feet. A comparison of measured and 

simulated heads (fig. 16) shows that the general 

geometry and elevation of the water table is simulated 

well at the regional scale. At any particular well, 

however, simulated and measured heads will differ. 

Fitted error statistics (Hill, 1998) for head residuals 

in the steady-state calibration indicate that the fitted 

standard deviation for heads is about 76 feet. This 

means that about 68 percent of the simulated heads are 

within 76 feet of the measured values, and that approximately 95 percent are within 152 feet. This error must 

be kept in mind when evaluating the head fluctuations 

simulated by the transient model. 

Measured temporal head fluctuations in the 

upper Deschutes Basin during the calibration period 

range from less than 1 foot to as much as 30 feet. 

Most measured fluctuations are in the 5 to 10 feet 

range. This range of fluctuations is small compared 

to the steady-state head residuals. It is unlikely, and 

unexpected, that the transient calibration would 

improve on the fit to absolute head elevations. The 

principal goal in transient calibration is to match 

observed head fluctuations. Therefore, when assessing 

the transient calibration with respect to heads, comparing simulated and observed fluctuations is more 

important than comparing absolute head elevations. 

To aid in comparing water-level fluctuations, graphs 

in this section have separate vertical axes for measured 

and simulated water levels. The axes have the same 

vertical scale but may be offset; the amount of offset 

can be determined by comparing the vertical axes on 

the graph. 

For the purposes of discussing simulated waterlevel fluctuations, the modeled area is divided into 

12 subareas (fig. 20). The subareas are geographically 

distinct regions within which wells generally show 

similar water-level fluctuation patterns. Each subarea 

contains 2 to 10 wells with time-series water-level 

data. In the following paragraphs, the subareas are 

discussed in an approximate north to south order.
51

Figure 20. Subareas used for discussion of water-level fluctuations. 
52

(1) Madras/Gateway Area

The Madras/Gateway area represents the entire 

model area north and east of Round Butte Dam 

(fig. 20, area 1). There are six wells with time-series 

water-level data in this area, representing the top 

four layers in the model. The general climate-driven 

pattern of water levels in this area, exemplified by well 

09S/14E-20DDA (fig. 21), is a slight decline starting 

in the early 1980s followed by an abrupt rise starting 

in late 1996 or early 1997. Three wells in the area (typified by well 11S/12E-26AAC, fig. 21) do not follow 

this pattern but, contrary to climatic trends, exhibit a 

rising water level from the mid-1980s to the mid1990s. The cause for this anomalous rise is not known, 

but it may be related to long-term effects of the filling 

of Lake Billy Chinook (Gannett and others, 2001). 

Because there is no process or boundary condition 

within the model to cause the anomalous rise, such as 

the filling of Lake Billy Chinook, it was not possible 

to simulate the rise. Simulated and observed waterlevel fluctuations match reasonably in wells such as 

09S/14E-20DDA and 11S/13E-11DBB (fig. 21) 

where the fluctuations are driven by climatic factors. 

Figure 21. Simulated and observed water-level fluctuations 

in the Madras/Gateway area. 

(2) Lower Desert Area

The Lower Desert area includes the region 

just south of the Metolius River arm of Lake Billy 

Chinook (fig. 20, area 2). There are three wells in this 

area with time-series water-level data, all of which 

exhibit anomalous behavior related to Lake Billy 

Chinook (see Caldwell and Truini, 1997, p. 22). Two 

of the wells exhibit the anomalous long-term rise 

contrary to climatic conditions described in the last 

section. The third well, for which only about 3 years 

of data are available, shows anomalous fluctuations 

likely related to either pumping or stage changes in 

the lake. Because of the anomalous behavior of these 

wells, comparisons with simulated fluctuations is not 

considered meaningful.
53

(3) Black Butte/Upper Metolius River Area

The Black Butte/upper Metolius River area 

extends from the uppermost Metolius River to south of 

Black Butte (fig. 20, area 3). The two wells with timeseries data, 14S/09E-08ABA and 13S/09E-03AAC02 

(fig. 22), represent model layers 1, 3, and 4. Data from 

the wells start in 1993 and 1994, respectively, so only 

the latter part of the calibration period is represented. 

Water levels in the area were generally declining until 

early 1995, when a snowmelt event caused a slight 

rise in water level. Water levels rose sharply during 

the winters of 1996 and 1997. Simulated water 

levels show both the early decline and subsequent 

water-level rises. The timing and magnitude of the 

measured and simulated fluctuations match reasonably 

well in the Black Butte/upper Metolius River area. 

Well 14S/09E-08ABA is particularly important, 

because it is one of the few wells representative 

of conditions in the Cascade Range, the principal 

recharge area. 

Figure 22. Simulated and observed water-level fluctuations in 

the Black Butte/upper Metolius River area. 

(4) Peninsula Area

The Peninsula area includes the region between 

the Deschutes and Crooked Rivers north of Lower 

Bridge (fig. 20, area 4). Time-series water-level measurements are available for five wells in the area, but 

data only go back to 1993, so only the later part of the 

calibration period is covered. One well that has less 

than 2 years of data and an anomalous shallow water 

level is not considered further. Another well showed 

fluctuations dominated entirely by nearby canal leakage and so is also not considered useful for regional 

model calibration. The remaining three wells, which 

represent the top three model layers, showed slight 

declines (0.5 to 2.0 feet) between 1994 and 1996, 

followed by a general rise in water levels starting 

in late 1996 in response to wetter conditions (fig. 23). 

The wells also showed slight seasonal fluctuations 

of generally less than one-half a foot. Simulated head 

fluctuations match both the measured long-term 

decline and rise and also the seasonal variations. 

The simulated rise in water levels starting in 1996, 

however, is generally larger than that observed. The 

timing of fluctuations are simulated accurately in well 

13S/12E-21AAC, but less so in well 14S/12E-02CCC.

Figure 23. Simulated and observed water-level fluctuations 

in the Peninsula area. 
54

(5) Lower Squaw Creek/Deep Canyon Area 

The lower Squaw Creek/Deep Canyon area 

includes the region between the Deschutes River 

and Squaw Creek, extending south to about Cline 

Falls on the Deschutes River and to about Henkle 

Butte on Squaw Creek (fig. 20, area 5). Data were 

available for five wells in the area representing the 

top three model layers. All wells include measurements from 1993 to 1997, and three include some 

measurements from 1978 and 1979. There are no 

measurements in the area from the 1980s and early 

1990s. As depicted by wells 15S/12E-09BDB and 

14S/11E-18CAC (fig. 24), measurements show 

that water levels in the area declined approximately 

5 feet between the late 1970s and mid-1990s. The 

decline continued until 1996, when water levels 

stopped declining and started to rise in response 

to the prevailing wet conditions. Simulated waterlevel fluctuations also show the long-term decline 

followed by the rise in response to wetter conditions. 

Simulated fluctuations generally match observed 

fluctuations well in both timing and magnitude 

in the lower Squaw Creek/Deep Canyon area. 

Figure 24. Simulated and observed water-level fluctuations in 

the lower Squaw Creek/Deep Canyon area. 

(6) Redmond Area 

The Redmond area, as used herein, refers to the 

region between the Deschutes River and Dry River that 

extends from Terrebonne south to about Long Butte 

(fig. 20, area 6). There are 11 wells with water-level 

fluctuation data available in this area representing 

the top four model layers. One of the wells, 15S/13E18ADD (fig. 25) provides a continuous record of 

water-level trends throughout the calibration period. 

This well shows that the water table in the area rose 

slightly between the late 1970s and late 1980s, and 

then declined until late 1996 or 1997. Other wells 

in the area with shorter data sets, 15S/12E-14CDD 

and 14S/13E-28CBB01, show the end of the decline 

with more detail, and also show seasonal fluctuations. 

Simulated water levels match both the timing and 

magnitude of long-term water-level fluctuations in 

the Redmond area.

Figure 25. Simulated and observed water-level fluctuations in the 

Redmond area. (A “p” indicates measurements affected by pumping.) 
55

Simulated and observed water-table fluctuations also 

match at shorter time scales. The simulated water 

levels in well 14S/13E-28CBB01 (fig. 25) match both 

seasonal fluctuations, which occur locally in response 

to irrigation and canal operation, and the longer term 

climatic signal (the year-to-year variation) in the 

observed data. The seasonal fluctuations observed 

in well 15S/12E-14CDD are out of phase with the 

irrigation-caused fluctuations, and are likely due 

to fluctuations in river stage. Because river-stage 

variations are not included in the model, simulated 

water levels in well 15S/12E-14CDD do not show this 

seasonal fluctuation, but only the longer-term climatedriven fluctuation. This well is largely unaffected by 

canal operation and irrigation.

(7) Sisters/Tumalo Area

The Sisters/Tumalo area extends from the region 

around the City of Sisters and McKinney Butte southeast to Tumalo (fig. 20, area 7). There are 10 wells in 

this area with time-series data representing mostly the 

top three model layers. One well is open to layers 5 

and 6. Well 15S/10E-08ACD (fig. 26) has continuous 

record throughout the calibration period. Three of the 

wells have some data from the late 1970s and mid1990s, but no data from the 1980s and early 1990s. 

Wells with long-term data show that water levels in the 

area were declining during the late 1970s. Water levels 

rose from the early 1980s until about 1986, which 

marked the start of a decline that lasted until about 

1996. Water levels began to rise again in 1996 in 

response to the return to wet climatic conditions.

Simulated water levels match the timing of 

climate-driven changes in observed water levels reasonably well, usually within a year (fig. 26). However, 

the simulated water levels show a long-term downward 

trend that is steeper than actually observed. This mismatch is seen in a few other places in the model, but 

seems to be restricted to the upper few layers. The 

reasons for the mismatch early in the calibration period 

could be due to initial conditions. The initial conditions 

chosen gave the best overall model fit and were consistent with climatic conditions immediately prior to the 

calibration period. The initial conditions are, however, 

steady state and may not reflect the spatial variability 

of true initial conditions throughout the modeled 

area. As typified by wells 15S/10E-36AAD02 and 

15S/11E-34CCA, the match between simulated 

and observed water-level trends during the 1990s 

is generally quite good in the Sisters/Tumalo area. 

Figure 26. Simulated and observed water-level fluctuations 

in the Sisters/Tumalo area. 

(8) Powell Butte Area

The Powell Butte area includes the region 

around and north of the town of Powell Butte 

(fig. 20, area 8). There are two wells in this area with 

time-series data from 1994 into 1997. Measurements 

in the area show stable or slightly declining water 

levels with seasonal fluctuations of 1 to 2 feet. Neither 

of the wells appears to have responded to the wet 

winters of 1996 and 1997. Simulated water levels 

show larger than measured seasonal fluctuations and 

a rise in water levels in response to the wet winters. 

Because both of these wells are near the margin of 

the pre-Deschutes Formation rocks, and on the model 

boundary, their value in assessing model calibration is 

unclear.
56

(9) Lower Tumalo Creek/Bend Area 

The lower Tumalo Creek/Bend area refers 

to the vicinity of the City of Bend, Awbrey Butte, 

and lower Tumalo Creek (fig. 20, area 9). There are 

four wells with time-series data representing model 

layers 1, 2, 4, 5, and 7. No wells have continuous 

observations throughout the calibration period. 

There are measurements spanning the mid-1990s 

from all four wells, and some measurements from the 

late 1970s from two wells, 17S/11E-25CBA and 

17S/12E-23BBC (fig. 27). Water levels in the area 

were declining during the late 1970s. No data are 

available for the 1980s and early 1990s, but based 

on data from other areas, water levels probably rose 

during the first one-half of the 1980s and then started 

to decline until the mid-1990s. Simulated water 

levels match the general trends observed in the 

measurements. The decline in simulated water levels 

between the late 1970s and the early 1990s is steeper 

than measured in the shallow layers (well 17S/11E25CBA) and less steep than measured in the deep 

layers (well 17S/12E-23BBC). There are also slight 

short-term (year-to-year) variations in the measurements, probably due to pumping and/or nearby canal 

operation not present in the simulated water levels. 

Figure 27. Simulated and observed water-level fluctuations 

in the lower Tumalo Creek/Bend area. 

(10) Benham Falls Area

The Benham Falls area includes the region 

around Benham Falls and Lava Butte, and extends 

north to Lava Island (fig. 20, area 10). There are six 

wells with time-series data in this area representing 

the top four model layers. None of the wells has 

continuous record through the calibration period, 

all have measurements from the early to mid-1990s, 

and three also have measurements in the late 1970s. 

Measurements indicate that water levels in the area 

were stable or showed only very slight declines 

between the late 1970s and the 1990s (fig. 28). 

Figure 28. Simulated and observed water-level fluctuations 

in the Benham Falls area. 
57

This is in contrast to many other regions, where water 

levels have fluctuated to a much larger degree in 

response to climate cycles. The stable water levels 

likely result from the moderating influence of the 

Deschutes River, which leaks water across this region 

at a long-term rate that is comparatively constant 

relative to climate variations. Wells in this area do, 

however, exhibit seasonal fluctuations of up to 6 to 

10 feet (wells 19S/11E-24CDA and 19S/11E-16ACC), 

with fluctuations generally largest closest to the river. 

Gannett and others (2001) show that water-level 

fluctuations in well 19S-11E-16ACC result from 

variations in river stage caused primarily by releases 

from Wickiup Dam. 

The match between simulated and observed 

water-level fluctuations in the Benham Falls area 

is poor. In general, simulated water levels show 

climate-driven decline in water levels of about 5 to 

10 feet that is not observed in the data (fig. 28, wells 

19S/11E-24CDA and 18S/11E-36ACB). The model 

is probably not reflecting the moderating effect of 

leakage from the Deschutes River on long-term 

fluctuations. In addition, simulated water levels do 

not show the seasonal fluctuations in the measured 

water levels. This is expected because the streamstage variations that drive water-table fluctuations in 

this area are not present in the model. Simulated water 

levels do show the water-level rise caused by the onset 

of wet climatic conditions seen in the measured data. 

(11) Central La Pine Subbasin

There are five wells with time-series data in the 

central La Pine subbasin (fig. 20, area 11). Three wells 

represent model layer 1, one well is open to model 

layers 5 and 6, and one well is open to model layer 8. 

One well has measurements for the entire calibration 

period, and a second has measurements from the mid1980s to the mid-1990s. Long-term data from well 

21S/11E-19CCC (fig. 29) show water-level trends 

similar to those observed elsewhere in the basin. 

Water levels rose from the early to mid 1980s, and 

then started a decline that lasted until the mid-1990s. 

The return to wet weather in 1996 caused water levels 

to start rising again. The character of water-level 

fluctuations in the La Pine subbasin varies with depth. 

Seasonal fluctuations that are prominent in shallow 

wells, such as well 21S/11E-19CCC, are virtually 

lacking in deep wells, such as 22S/10E-10DDB. The 

magnitude of longer term climate-driven fluctuations 

appears to be subdued at depth as well. 

Simulated water-level fluctuations match those 

observed in most regards. The magnitude of the longterm, climate-driven trends are reasonably matched in 

both shallow and deep wells. In well 22S/10E-10DDB 

(fig. 29), open to model layer 8, the maximum and 

minimum simulated water levels appear to lag the 

observations by about 1 year. This lag does not 

appear to be present in well 22S/10E-14CCA, which 

is open to model layers 5 and 6. Both the timing and 

magnitude of long-term fluctuations are simulated at 

shallow depths (well 21S/11E-19CCC). The magnitude of annual fluctuations is also simulated well 

at shallow depths. The match between the timing 

of simulated and observed annual water-level fluctuations at shallow depths is variable. This is due to 

the limitations of the semiannual stress periods in the 

model. 

Figure 29. Simulated and observed water-level fluctuations 

in the central La Pine subbasin. 

Observed

Simulated

Open interval—

 layer 2 

Observed

Simulated

Open interval—

 layer 8

Observed

Simulated

Open interval—

 layers 5 and 6

22S/10E-14CCA

4170

4175

4180

4185

4190

4270

4275

4280

4285

4290

4180

4190

4200

4210

4280

4290

4300

4310

22S/10E-10DDB

21S/11E-19CCC

4170

4180

4190

4200

4210

4190

4200

4210

4220

4230

P

OBSERVED WATER-LEVEL ELEVATION, IN FEET

SIMULATED WATER-LEVEL ELEVATION, IN FEET

1975 1980 1985 1990 1995 2000

YEAR
58

(12) Southern La Pine Subbasin

There are two wells with time-series data in 

the southern part of the La Pine subbasin (fig. 20, 

area 12). One is open to model layer 2, and the second 

is open to model layers 2 through 5. Measurements are 

available for both wells only for the mid-1990s, 

and no data are available on long-term water-level 

trends in the area. It is likely that long-term trends 

mirror the climate-driven trends seen elsewhere 

in the basin. Both wells in the area show apparent 

climate-driven water-level rise on the order of 2 feet 

during 1996 and 1997 as typified by well 23S/09E36BBC (fig. 30). Neither well shows any annual 

fluctuation. 

Simulated and observed water-level fluctuations 

do not match well in this part of the model. Simulated 

water levels exhibit annual fluctuations not observed 

in the measurements. Moreover, simulated water 

levels show a greater year-to-year variation than 

observed. Data from elsewhere in the La Pine subbasin 

show that shallow wells exhibit pronounced annual 

water-level fluctuations and that fluctuations in deeper 

wells are much more subdued. The observed water 

levels in this area are similar to those in deep wells 

elsewhere in the La Pine subbasin, while the simulated 

water levels are similar to those measured in shallow 

wells elsewhere in the subbasin. In other words, the 

model is not properly simulating the variation in 

water-level fluctuations with depth in this area. This 

should not be a problem with the use of the model, 

since this area is on the periphery of the model and far 

from critical stream reaches and areas of likely water 

use.

Figure 30. Simulated and observed water-level fluctuations in the 

southern La Pine subbasin. 

Summary of Simulated Water-Level Fluctuations

Water-level fluctuations in the upper Deschutes 

Basin are primarily driven by climate cycles, canal 

leakage and irrigation, and variations in stream stage. 

Water-level fluctuations can also be caused by pumping. The upper Deschutes Basin ground-water model 

simulates fluctuations caused by climate cycles, 

annual recharge cycles, and canal leakage and 

irrigation. Water-level fluctuations caused by streamstage variations are not simulated.

Climate-driven water-level fluctuations are 

simulated with reasonable accuracy over most of the 

model area. The timing of the simulated water-level 

response to the change to wet conditions starting in 

1996 matches the observed response within months, 

in most cases, and generally within a year. The timing 

and magnitude of simulated water-level fluctuations 

caused by annual pulses of recharge from precipitation 

match those observed reasonably well, given the limitations of the time discretization in the model. Waterlevel fluctuations caused by annual canal leakage are 

simulated very well over most of the area where this 

occurs. Time discretization was not a problem in simulating canal-leakage effects because the stress periods 

in the model were designed to approximate the irrigation season.

There are some areas in the model where 

simulated long-term water-level declines exceed the 

observed decline (fig. 26). This can be the case even 

where the simulated and observed short- and longterm fluctuations match well. This appears to be 

an artifact of the starting heads and of starting with 

steady-state conditions.

Comparison of Simulated and Observed Transient 

Ground-Water Discharge to Streams

Streamflow data were used to estimate temporal 

variations in ground-water discharge to seven stream 

reaches in the modeled area representing major 

ground-water discharge areas. The ability to accurately simulate ground-water discharge variations 

to some reaches was limited by the scale of the 

model grid relative to the scale of the stream network. 

This was a particular problem with streams above 

Crane Prairie Reservoir, where model cells are up 

to 10,000 feet on a side and streams are often less 

than a mile apart. In such areas, it is more meaningful 

to look at ground-water discharge to regions of the 

model than to evaluate individual stream reaches. 

Observed

Simulated

Open interval—

 layers 2-5

23S/09E-36BBC

4280

4290

4300

4310

4320

4335

4345

4355

4365

4375

OBSERVED WATER-LEVEL ELEVATION, IN FEET

SIMULATED WATER-LEVEL ELEVATION, IN FEET

1975 1980 1985 1990 1995 2000

YEAR
59

The ability to simulate the exact timing of observed 

annual variations in ground-water discharge was 

limited by the time descretization in the model. The 

semiannual stress periods in the model were the same 

each year of the simulation: November to April and 

May to October. The actual timing of recharge from 

winter storms and snowmelt varies from year to year. 

Comparison of simulated and observed year-to-year 

variations is less affected by this timing. Simulated 

temporal variations in ground-water discharge are 

compared with estimated variations based on streamgage data in the following paragraphs. 

Odell Creek

Temporal variations in ground-water discharge 

to the Odell Creek drainage were estimated from 

streamflow gage data from Odell Creek (USGS station 

number 14055600). Streamflow at this point contains 

both ground-water and surface-water components, 

including flow from Odell Lake. A hydrograph of 

daily and monthly mean flows (fig. 31) shows many 

individual storm peaks, in addition to the annual 

spring runoff peaks. Baseflow, the ground-water 

component of flow inferred from low-flow periods on 

the hydrograph, ranged between 50 and 200 ft3/s during the calibration period. Simulated ground-water 

discharge to the Odell Creek drainage above the gage 

(fig. 31) closely approximates the baseflow inferred 

from the stream hydrograph in both magnitude and 

timing of variations.

Figure 31. Daily and monthly mean flows of Odell Creek, and simulated ground-water discharge, 1977 to 1997. 

0

100

200

300

400

500

600

700

800

STREAMFLOW OR DISCHARGE,

IN CUBIC FEET PER SECOND

Observed daily mean streamflow

Observed monthly mean streamflow

Simulated ground-water discharge

1975 1980 1985 1990 1995 1999

YEAR

Odell Creek

Streams above and below Crane Prairie Reservoir

Streams above Crane Prairie Reservoir are considered together because the scale of the stream network in this area is small compared to that of the 

model grid. Because of the scale discrepancy, not 

all streams above Crane Prairie Reservoir could be 

included in the model. Cultus Creek and the Deschutes 

River were represented in the model because they 

spanned the region drained by ground-water fed 

streams in the area. 

Measured monthly mean flows and simulated 

baseflow of Cultus Creek are shown in figure 32. The 

flow of Cultus Creek responds rapidly to snowmelt 

and then goes dry (or nearly so) each fall, indicating 

that it is dominated by surface runoff and shallow 

ground-water discharge. Simulated baseflow (fig. 32) 

exceeds actual baseflow (as inferred from monthly 

mean values for low-flow months), particularly during 

the fall. The amplitude of the simulated baseflow variation is significantly less than that measured. Patterns 

of long-term baseflow variation are similar. Cultus 

River, a ground-water dominated stream with a mean 

annual discharge of 63 ft3

/s, is roughly a mile from 

Cultus Creek. Because Cultus River is not explicitly 

represented in the model, Cultus Creek, as simulated, 

must represent at least part of the ground water that 

actually discharges to Cultus River. Therefore, it is 

reasonable that simulated baseflow to Cultus Creek 

should exceed that actually measured. 
60

Figure 32. Monthly mean flow of Cultus Creek, and simulated ground-water discharge, 1977 to 1997. 

0

50

100

150

Observed monthly mean streamflow

Simulated ground-water discharge

Cultus Creek

STREAMFLOW OR DISCHARGE,

IN CUBIC FEET PER SECOND

1975 1980 1985 1990 1995 1999

YEAR

Ground-water discharge in the vicinity of the 

uppermost Deschutes River can be evaluated from 

the streamflow gage below Snow Creek (USGS 

station number 14050000), just above Crane Prairie 

Reservoir. The hydrograph of daily and monthly mean 

flows (fig. 33) shows that the Deschutes River at this 

point is largely ground-water fed with little response to 

individual storms. The annual peaks caused by snowmelt are much larger than in other ground-water fed 

streams in the area, suggesting that the stream is fed 

in part by a relatively shallow ground-water system. 

Simulated ground-water discharge to the Deschutes 

River is lower than observed (fig. 33). Given that the 

Deschutes River must represent (in part) discharge to 

Cultus River as well, it is apparent that overall groundwater discharge in the area of the model above Crane 

Prairie reservoir is too low by roughly 50 to 100 ft3/s. 

This limits the usefulness of the model for evaluating 

ground-water discharge fluctuations in this area.

Figure 33. Daily and monthly mean flows of the Deschutes River above Snow Creek, and simulated ground-water discharge, 1977 to 1997.

0

100

200

300

400

500

Observed daily mean streamflow

Observed monthly mean streamflow

Simulated ground-water discharge

Deschutes River above Snow Creek

STREAMFLOW OR DISCHARGE,

IN CUBIC FEET PER SECOND

1975 1980 1985 1990 1995 1999

YEAR

Simulated ground-water discharge volumes 

below Crane Prairie Reservoir and above Davis Creek 

appear reasonable. Although there are no stream-gage 

data to estimate temporal variations in ground-water 

discharge to this reach of the Deschutes River, the 

simulated discharge compares favorably with the 

estimated average discharge from 1993 to 1995 used 

in the steady-state calibration (fig. 34).

Figure 34. Simulated ground-water discharge to the Deschutes River between Crane Prairie Reservoir and Davis Creek, and estimated 

long-term ground-water discharge, 1977 to 1997. 

0

100

200

300

400

Estimated long-term average ground-water discharge

Simulated ground-water discharge

Deschutes River between Crane Prairie Reservoir and Davis Creek

DISCHARGE,

IN CUBIC FEET PER SECOND

1975 1980 1985 1990 1995 1999

YEAR
61

Fall River

Fall River is entirely spring fed, so the gage 

record (USGS station 14057500) directly represents 

ground-water discharge. Simulated ground-water 

discharge to Fall River exceeds measured discharge by 

roughly one-third (about 50 ft3/s) (fig. 35). The magnitude and timing of discharge variations, however, is 

considered good given the constraints of model grid 

scale and stress-period definition. 

Figure 35. Monthly mean flow of Fall River and simulated ground-water discharge, 1977 to 1997. 

0

100

200

300

Fall River

STREAMFLOW OR DISCHARGE,

IN CUBIC FEET PER SECOND

1975 1980 1985 1990 1995 1999

YEAR

Observed monthly mean streamflow

Simulated ground-water discharge

Metolius River Drainage

A hydrograph of monthly mean discharge of the 

Metolius River just upstream of Lake Billy Chinook 

(USGS station number 14091500) (fig. 36) shows 

that seasonal variations are small relative to the total 

discharge, as is typical of ground-water dominated 

streams. The Metolius River accounts for a large 

part of the total ground-water discharge in the upper 

Deschutes Basin. Much of the discharge is to springs 

along the main stem above Jefferson Creek and to 

lower parts of tributaries. Temporal variations in 

ground-water discharge to the Metolius River are difficult to estimate precisely because the only long-term 

gaging station on the river is low in the drainage, just 

upstream of Lake Billy Chinook. The streamflow at 

this point includes flow from both ground-water and 

surface-water dominated tributaries, some of which 

include glacial meltwater. It is clear, however, that a 

large proportion of the total flow of the Metolius River 

comes from ground water, and the flow seldom drops 

below 1,200 ft3

/s. October mean flows can be taken 

to represent regional ground-water discharge to the 

Metolius River. To correct for possible glacial contributions, the flow of Jefferson Creek and the Whitewater River (gaged since 1984 and 1983, respectively) 

can be subtracted from the flow of the Metolius River.

Simulated ground-water discharge matches 

October mean flows (corrected for Whitewater River 

and Jefferson Creek) reasonably well, particularly 

with respect to decadal-scale fluctuations (fig. 36). 

The notable departure between simulated groundwater discharge and October mean flows occurs in 

1996 and 1997. The extremely wet winters of 1996 

and 1997 caused a substantial increase in flow from 

the Metolius River drainage. It is apparent that during 

the fall of both 1996 and 1997 the flow of the Metolius 

River was still receding from the previous winter.
62

Figure 36. Monthly mean and October mean flow of the Metolius River, October mean flow of the Metolius River minus the October mean 

flows of Whitewater River and Jefferson Creek, and simulated ground-water discharge, 1977 to 1997. 

0

1,000

500

1,500

2,000

2,500

3,000

Observed monthly mean streamflow

Observed October mean streamflow

Observed October mean flow minus Whitewater River and Jefferson Creek

Simulated ground-water discharge

STREAMFLOW OR DISCHARGE, IN CUBIC FEET PER SECOND

1975 1980 1985 1990 1995 1999

Metolius River

YEAR

Deschutes River—Bend to Culver

Temporal variations in ground-water discharge 

to the middle section of the Deschutes River can be 

evaluated by comparing data from gaging stations 

below Bend (USGS station number 14070500) and 

near Culver (USGS station number 14076500). During 

the summer months, tributary inflow to this reach 

is negligible, and the increase in flow of the river 

is attributable entirely to ground-water discharge 

(including discharge to the lower 2 miles of Squaw 

Creek). The amount of ground-water discharge can 

be determined by subtracting the flow at Bend from 

the flow at Culver during summer months. During 

the winter months, however, ungaged tributary inflow 

precludes the use of gage data to estimate groundwater discharge. Gage data, therefore, can be used to 

evaluate the year-to-year variations in ground-water 

discharge, but not seasonal variations.

The magnitude and timing of year-to-year 

variations in simulated fall ground-water discharge 

match observed variations in August mean groundwater discharge reasonably well (fig. 37). The total 

rate of simulated ground-water discharge exceeds 

observed discharge by roughly 100 to 150 ft3/s. 

The larger discrepancy early in the simulation period 

is likely an artifact of starting conditions. The difference late in the simulation period is consistent with 

the steady-state calibration. 
63

Figure 37. Measured monthly mean inflow to the Deschutes River between Bend and Culver, 1977 to 1997, and simulated inflow for all 

time steps. (Inflow includes both surface-water and ground-water components during winter months. Measured August mean inflow 

and simulated fall inflow is due almost entirely to ground-water discharge.) 

0

100

200

300

400

500

600

700

800

900

1,000

Observed monthly mean inflow (includes surface water and ground water)

Observed August mean inflow (primarily ground water)

Simulated inflow all time steps (includes surface water and ground water)

Simulated fall inflow (primarily ground water)

INFLOW, IN CUBIC FEET PER SECOND

1975 1980 1985 1990 1995 1999

Deschutes River between Bend and Culver

YEAR

Lower Crooked River

The lower Crooked River from Osborne Canyon 

to the mouth is one of the principal ground-water 

discharge areas in the upper Deschutes Basin, with 

ground-water discharge exceeding 1,000 ft3/s. 

Ground-water discharge to the lower Crooked River 

was estimated using data from the streamflow gage 

below Opal Springs. During the late summer and 

early fall, the flow below Opal Springs consists 

almost entirely of ground-water discharge. The late 

summer/early fall flow of the Crooked River below 

Opal Springs has been used to evaluate long-term 

changes in ground-water discharge to the river (Sceva, 

1968; Gannett and others, 2001). This can be done 

because long-term changes in ground-water discharge 

are large compared to variations in late-summer/early 

fall flow from upstream. Flow measurements below 

Opal Springs are of limited use in evaluating seasonal 

variations in ground-water discharge, however, because 

these variations are small compared to variations in 

flow from upstream. In order to evaluate short-term 

variations in ground-water discharge to the lower 

Crooked River, flow from upstream must be accounted 

for. Measurements of streamflow upstream of the 

ground-water discharge area are not available for the 

simulation period; however, insight can be gained 

from measurements available for a few years prior to 

the calibration period.

Streamflow just upstream of the ground-water 

inflow area (near Terrebonne) was measured from 

October 1968 to July 1973. These measurements can 

be subtracted from the flow at Opal Springs to evaluate 

short-term variations in ground-water discharge to the 

lower Crooked River. A comparison of monthly mean 

discharge between the two gages (fig. 38) shows that 

almost all of seasonal flow variation observed at the 

Opal Springs gage can be attributed to surface-water 

flow from upstream. Ground-water discharge, calculated as the difference between the two gages, shows 

much less variation. The short-term variations in the 

difference between the gages correspond exactly to 

the variations in streamflow. This suggests that the 

variations in the difference may not be due to variations 

in ground-water discharge, but are likely related to 

ungaged tributary inflow or systematic measurement 

error. A comparison of July mean flows at Terrebonne 

and Opal Springs (fig. 39) shows that even the yearto-year variations in discharge at Opal Springs are 

due predominantly to surface water upstream of the 

ground-water discharge area. In summary, variations in 

regional ground-water discharge to the lower Crooked 

River are small on a year-to-year basis and most pronounced at longer (decadal) time scales. This is to be 

expected given the scale of the ground-water system 

contributing to discharge at this point.

Simulated transient ground-water discharge to 

the lower Crooked River is consistent with inferences 

made from gage data (fig. 40). Simulated ground-water 

discharge shows small year-to-year variations and 

only slightly larger long-term (decadal) variation. The 

amount of simulated ground-water discharge is very 

close to the average late summer/early fall streamflow 

at Opal Springs during the calibration period.
64

Figure 38. Monthly mean flows of the Crooked River 

below Opal Springs and at Terrebonne, and the 

difference between the two gages, 1967 to 1973. 

0

500

1,000

1,500

2,000

2,500

3,000

3,500

4,000

4,500

1967 1968 1969 1970 1971 1972 1973

Monthly mean stream- 

 flow—Crooked River 

 at Terrebonne

Monthly mean stream- 

 flow—Crooked River 

 below Opal Springs

Difference—ground-

 water discharge

STREAMFLOW, IN CUBIC FEET PER SECOND

YEAR

Figure 39. July mean flows of the Crooked River 

below Opal Springs and at Terrebonne, and the 

difference between the two gages, 1968 to 1973. 

0

200

400

600

800

1,000

1,200

1,400

July mean streamflow— 

Crooked River at Terrebonne

July mean streamflow— 

Crooked River below Opal Springs

Difference—ground-water discharge

1967 1968 1969 1970 1971 1972 1973

STREAMFLOW, IN CUBIC FEET PER SECOND

YEAR

Figure 40. July mean flows of the Crooked River below Opal Springs and simulated ground-water discharge to the lower Crooked River, 

1977 to 1997. 

1,000

1,500

Observed July mean streamflow

Simulated ground-water discharge

STREAMFLOW OR DISCHARGE,

IN CUBIC FEET PER SECOND

Crooked River below Opal Springs

1975 1980 1985 1990 1995 1999

YEAR
65

Summary of Transient Ground-Water Discharge to Streams

The transient model simulates the volumetric 

distribution and temporal variations in ground-water 

discharge to streams reasonably well. The match 

between simulated and measured volume of and variations in ground-water discharge is somewhat dependent on geographic scale. Simulated and observed 

discharge fluxes and variations are most similar for 

streams that receive regional ground-water discharge. 

At smaller scales (as would be represented by firstorder tributaries), geologic heterogeneities, topographic complexity, and model discretization combine 

to reduce the accuracy of matches between observed 

and simulated discharge fluxes and variations. 

Nevertheless, simulated fluxes and variations are very 

close for certain small streams, such as Odell Creek. 

The fit between observed and simulated discharge fluctuations is also somewhat dependent on 

temporal scales. Overall, the fit between simulated and 

observed discharge fluctuations is best at decadal time 

scales, which reflect long-term climate cycles (such 

as droughts). The ability of the model to match the 

timing of annual fluctuations is somewhat limited 

by the time discretization, particularly the semiannual 

stress periods, as previously discussed. 

Transient Model Water Budget

The transient model water budget can be evaluated by observing the instantaneous rates of various 

budget components at the end of each stress period 

(table 7 and fig. 41). The most noticeable feature 

of the transient budget is the way the ground-water 

system attenuates the wide seasonal and interannual 

variations in recharge. Recharge varies interannually 

by a factor of nearly 5 while stream discharge varies 

by a factor of only about 0.2 (which is consistent 

with streamflow measurements). The large swings 

in recharge are moderated by storage. During the 

wet winters of 1982 to 1984, it can be seen that large 

amounts of ground water go into storage (as the water 

table rises), and the amount of ground water going into 

storage during the winter exceeds the amount coming 

out of storage the following summers. This situation is 

reversed during the dry winters of 1990 to 1992, when 

little water goes into storage and a much larger amount 

comes out of storage the following summers (as the 

water table drops). The variation in stream discharge 

over the calibration period is consistent with observed 

variations in discharge of the Deschutes River near 

Madras.

The transient model water budget shows a cumulative reduction in ground-water storage of about 

3.5 × 1010 ft3 (about 810,000 acre-ft) over the 1978 

to 1997 transient calibration period. This is equivalent 

to about 0.31 feet of water over the model area and 

is consistent with the water-level declines observed 

over most of the study area between the late 1970s 

and mid-1990s. 

Model Limitations

Numerical models of ground-water flow are 

only approximations of complex natural systems 

and, as such, have intrinsic error and uncertainty. 

Error stems largely from the fact that certain spatially 

variable properties, such as hydraulic conductivity 

and stream stage, must be represented as uniform 

values in discrete model cells. Simplification can also 

occur at larger scales within a model, as is the case 

here, where hydraulic conductivity is represented as 

uniform in zones composed of multiple model cells. 

Model uncertainty stems from random error in the 

field measurements used for model calibration, which 

is translated through model calibration to uncertainty 

in the calibrated parameter values.

Because of intrinsic error and uncertainty, the fit 

between simulated and observed hydraulic heads and 

fluxes, described in previous sections, is not perfect. 

For example, steady-state model head residuals have a 

fitted standard deviation of 76 feet, and a root-meansquare (RMS) error of 78 feet. The RMS error for 

heads in the transient calibration is 89.6 feet. This 

error affects the differences between simulated 

ground-water elevations and stream elevations and, 

consequently, the spatial distribution of simulated 

ground-water/surface-water exchanges. The comparison of simulated and measured stream gains and losses 

described in preceding sections, however, shows that 

the model fit is reasonably good, particularly at scales 

larger than a few to several miles. At smaller scales, 

the fit is not as good. The numerical model also has 

certain error with regard to transient phenomena. 

This can also be evaluated by assessing the fit between 

simulated and observed variations in hydraulic head 

and discharge to streams. In general, simulated and 

observed responses to regional stresses matched 

within a year or two. In the Redmond area, the model 

fit between simulated and observed responses to shortterm stresses, such as canal operation and irrigation, 

was good. 
66

Table 7. Transient model water budget for each stress period[All values in cubic feet per second; Nov, November; Apr, April; Oct, October; ET, evapotranspiration]

Stress 

period Season

Model 

water year

Recharge Storage

ET

Precipitation Canal On-farm Total In from Out to Change 

1 Nov–Apr 1978 7,605 104 13 7,722 124 3,627 3,503 48 2 May–Oct 1978 1,041 741 109 1,891 2,622 245 -2,377 622 3 Nov–Apr 1979 4,466 75 13 4,554 428 1,122 695 57 4 May–Oct 1979 1,279 790 109 2,179 2,304 242 -2,062 732 5 Nov–Apr 1980 5,093 104 13 5,210 289 1,622 1,334 42 6 May–Oct 1980 1,381 730 109 2,220 2,086 188 -1,898 624 7 Nov–Apr 1981 5,268 98 13 5,379 304 1,815 1,511 62 8 May–Oct 1981 1,500 675 109 2,284 2,007 157 -1,850 631 9 Nov–Apr 1982 10,814 88 13 10,915 57 6,362 6,305 48 

10 May–Oct 1982 2,821 757 109 3,687 1,362 522 -840 516 11 Nov–Apr 1983 8,983 86 13 9,083 202 4,408 4,206 37 12 May–Oct 1983 1,655 787 109 2,552 2,600 427 -2,172 674 13 Nov–Apr 1984 9,438 82 13 9,533 170 4,708 4,538 28 14 May–Oct 1984 914 813 109 1,837 3,411 376 -3,035 758 15 Nov–Apr 1985 4,849 117 13 4,978 266 878 612 60 16 May–Oct 1985 870 791 109 1,770 3,111 244 -2,867 790 17 Nov–Apr 1986 8,017 100 13 8,129 177 3,702 3,525 51 18 May–Oct 1986 874 758 109 1,740 3,222 227 -2,996 850 19 Nov–Apr 1987 4,203 130 13 4,346 493 736 243 73 20 May–Oct 1987 1,725 768 109 2,602 2,270 325 -1,946 807 21 Nov–Apr 1988 4,539 105 13 4,657 428 1,063 635 73 22 May–Oct 1988 626 697 109 1,433 3,121 148 -2,973 806 23 Nov–Apr 1989 7,480 101 13 7,594 256 3,667 3,411 32 24 May–Oct 1989 621 765 109 1,496 3,102 197 -2,905 723 25 Nov–Apr 1990 2,531 120 13 2,664 1,306 236 -1,070 84 26 May–Oct 1990 880 714 109 1,703 2,649 159 -2,489 774 27 Nov–Apr 1991 2,507 91 13 2,612 1,297 422 -875 57 28 May–Oct 1991 1,624 682 109 2,415 1,686 159 -1,527 606 29 Nov–Apr 1992 3,660 96 13 3,770 843 1,111 268 64 30 May–Oct 1992 910 591 109 1,611 2,378 105 -2,272 652 31 Nov–Apr 1993 7,160 54 13 7,227 284 3,832 3,547 15 32 May–Oct 1993 2,766 724 109 3,600 929 486 -443 525 33 Nov–Apr 1994 1,407 103 13 1,523 1,984 53 -1,931 158 34 May–Oct 1994 636 627 109 1,372 2,606 128 -2,478 677 35 Nov–Apr 1995 5,111 79 13 5,204 533 2,278 1,745 30 36 May–Oct 1995 1,841 658 109 2,608 1,329 171 -1,158 471 37 Nov–Apr 1996 7,982 81 13 8,076 270 4,418 4,149 60 38 May–Oct 1996 1,167 749 109 2,025 2,397 291 -2,106 638 39 Nov–Apr 1997 8,990 95 13 9,098 83 4,953 4,870 63 40 May–Oct 1997 1,665 725 109 2,500 2,019 333 -1,686 576 
67

Figure 41. Variations in the transient model water budget by stress period. 

-8,000

-6,000

-4,000

-2,000

0

2,000

4,000

6,000

8,000

10,000

12,000

0246 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40

STRESS PERIOD

CUBIC FEET PER SECOND

1978 1980 1982 1984 1986 1988 1990 1992 1994 1996

Storage Pumpage (X10)

Evapotranspiration Recharge

Stream Leakage In Stream Leakage Out

Model error and uncertainty propagate into 

model predictions. Therefore, predictions of the location and timing of streamflow and head changes based 

on simulations must be considered general indicators 

of how the ground-water system will respond to 

stresses and should not be considered exact. A good 

sense of the accuracy and reliability of model predictions can be obtained by evaluating the match between 

simulations and observations as is done in the preceding calibration sections of this report. 

EXAMPLE SIMULATIONS

The calibrated numerical model can be used to 

estimate the response of the regional ground-water 

system to new stresses, such as pumping wells or 

modifying recharge through canal lining. Both of these 

types of stresses involve changes in the hydrologic 

budget. In the case of pumping, ground water is 

removed from the system. In the case of canal lining, 

the recharge to the system is reduced. In both cases, 

the amount of pumping or reduced recharge must 

be offset by reductions in natural discharge from 

the ground-water system or increases in recharge 

elsewhere. Time is required for the ground-water 

system to compensate for these changes and reach 

a new equilibrium. During this time, changes in 

storage within the ground-water system occur as 

the water table changes to reflect the new equilibrium. 
68

When analyzing the results of simulations, important 

factors to evaluate include the location and rates of 

changes in stream leakage and spring discharge, 

changes in the water-table elevation, and changes in 

ground-water storage. 

Example simulations showing the effects of 

pumping and of canal lining are included in this 

section. When evaluating the results of these, or any 

simulations, it is important to remember that model 

predictions reflect the error and uncertainty in the 

model. Model predictions are best used as general 

guides to help understand how the ground-water 

system will respond to new stresses and should not be 

considered exact. A general sense of model accuracy 

and performance can be obtained from the comparison 

of observed and simulated heads and flows presented 

in the Model Calibration section of this report. 

Simulating the Effects of Pumping

Example simulations were run to demonstrate 

the effects of pumping wells in three areas. Two 

of the areas are the population centers of Bend and 

Redmond. The third area is an arbitrary location 

between the Crooked and Deschutes Rivers near 

Osborne Canyon used to demonstrate the difference in 

the type of response close to a major discharge area. 

At each of these areas, both shallow and deep wells 

are simulated.

In all simulations, wells are pumped at a 

continuous rate of 10 ft3/s for 50 years. The 50-year 

stress period is divided into 30 time steps. Successive 

time steps increase in length by a factor of 1.25, so 

they range in length from 5.7 days to 10 years. 

The simulations start with simulated steady-state 

conditions, and all stresses other than pumping are 

held constant. 

When pumping of a well begins, virtually all the 

water comes from storage as the hydraulic head in 

the area immediately around the well drops, forming a 

cone of depression. With time, the cone of depression 

expands and starts to stabilize as sufficient ground 

water to offset the pumping is intercepted, and the 

proportion of pumped water coming from storage 

decreases. As the proportion of pumped water from 

storage decreases, the proportion of water from 

streams or other sources increases. At equilibrium, 

no water comes from storage and all pumped water 

comes from other sources. In the upper Deschutes 

Basin, the principal source of water to pumped wells 

once equilibrium has been attained is diminished 

streamflow.

The cone of depression from a pumping well can 

affect streamflow in both losing and gaining reaches. 

Along losing stream reaches where the ground-water 

elevation is sufficiently close to the stream elevation, 

lowering the ground-water head increases the head 

gradient between the ground-water system and the 

stream, causing increased stream leakage. If the 

ground-water head is far below the stream, the 

ground-water system and the stream are effectively 

uncoupled, and changes in head have little or no effect 

on stream loss. Along gaining reaches, the amount 

of ground-water discharge to the stream can be diminished by lowering the head in the ground-water 

system and reducing the head gradient, and groundwater flow, toward the stream.

The only other possible sources of water to 

pumped wells in the upper Deschutes Basin are 

reduced ET and increased inflow at basin boundaries. 

Lowering the water table in areas where phreatophytes utilize water directly from the water table can 

cause a reduction in ET. This process is simulated in 

the model but is generally inconsequential. Increased 

boundary fluxes due to head changes are not simulated in this model. However, locations of probable 

pumping in the upper Deschutes Basin are sufficiently 

distant from boundaries (compared to the proximity to 

streams) that the probable contribution from boundaries is likely minuscule. The model should be used 

with caution when evaluating the effects of pumping 

of wells very close to basin boundaries, however, such 

as very near or on the Cascade Range crest, because 

simulated drawdowns may be larger than they would 

be in reality. 

Simulated Pumping in the Bend Area

The first pumping simulation is in model layer 2 

in Bend (row 87, column 50) (fig. 42). In this simulation, water comes from storage, increased stream 

losses, and decreased gains. Gaging-station data discussed in preceding sections show gaining and losing 

reaches near Bend, indicating shallow ground water 

in hydraulic connection with the Deschutes River. After 

2 years of pumping, about 2.3 ft3

/s or 23 percent of the 

pumping rate comes from diminished streamflow. 
69

Most of the effect is to the Deschutes River in the 

vicinity of Bend (fig. 42). The diminished reaches 

include both those affected by increased leakage 

(stream losses) and decreased spring discharge (stream 

gains). After 10 years of pumping, about 5.8 ft3/s, 

or 58 percent of the pumping rate, comes from diminished streamflow. The affected reaches are more widespread and include the confluence area, but are still 

concentrated in the Bend area. A graph depicting 

changes in the sources of water with time (fig. 43) 

shows that the amounts of water coming from storage 

and from diminished streamflow are equal after about 

7 years and that 90 percent of the pumped water comes 

from diminished streamflow after about 42 years. 

Maximum drawdown in model layer 2 after 2 years 

is 76 feet (fig. 42). It is much less in model layer 1 

(about 1 foot). The maximum drawdown in model 

layer 2 increases to only 80 feet after 10 years, but 

the effects are more widespread. Drawdown in model 

layer 1 is 4 feet after 10 years. 

Figure 42. Simulated drawdown in model layer 2 and net reduction in ground-water discharge to streams after 2 and 10 years of pumping 

10 ft3/s from model layer 2 in the Bend area. 

121o30' 121o00'

44o30'

44o00'

121o30' 121o00'

44o30'

44o00'

10 11 12 13 14 15

20

19

18

17

16

14

13

12

11

15

10 11 12 13 14 15

20

19

18

17

16

14

13

12

11

15

Millican

Sunriver

Bend

Redmond

Cline

Falls

Madras

Sisters

Metolius

Culver

Terrebonne

Lower

Bridge

Millican

Sunriver

Bend

Redmond

Cline

Falls

Madras

Sisters

Metolius

Culver

Terrebonne

Lower

Bridge

Arnold

Canal

Tumalo Feed Canal

Central Oregon 

Canal

North Unit

Main Canal

Squaw Creek

Canal

Arnold

Canal

Tumalo Feed Canal

Central Oregon 

Canal

North Unit

Main Canal

Squaw Creek

Canal

Lava

Butte

Cline

Buttes

Juniper

Butte

Osborne

Canyon

Gray

Butte

Lower Desert

Lava

Butte

Cline

Buttes

Juniper

Butte

Osborne

Canyon

Gray

Butte

Lower Desert

Study area

 boundary

Study area

 boundary

Spring

Haystack

Reservoir

Creek

Mud Springs

 Creek

Lake Billy

Willow

Chinook

Lake

Simtustus

Creek

Fly

Crooked

River

Metolius

Round

Butte

Dam

Houston

Lake

Dry

Creek

Cr

River

Squaw

River

River

Cr

Tumalo

DESCHUTES

RIVER

Spring

Haystack

Reservoir

Creek

Mud Springs

 

Willow

 Creek

Lake Billy

Chinook

Lake

Simtustus

Creek

Fly

Crooked

River

Metolius

Round

Butte

Dam

Houston

Lake

Dry

Creek

Cr

River

Squaw

River

River

Cr

Tumalo

DESCHUTES

RIVER

Indian

Ford

Indian

Ford

97

97

20

26

20

126

126

26

97

97

20

26

20

126

126

26

1

3

5

10

30

20

40

1

3

5

10

20

30

40

EXPLANATION

Streamflow reduction—In cubic feet per second per mile

 .01-.05 

 .05-.1

 More than .1

Line of equal drawdown—In feet. Contour intervals are 2, 5, and 10 feet.

Well location

10

0 5 10 MILES

0 5 10 KILOMETERS

2 YEARS 10 YEARS

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

 See table of contents for detailed information on sources.
70

Figure 43. Variations with time in the amounts of ground water 

coming from storage and from diminished streamflow in response 

to pumping 10 ft3

/s from model layer 2 in the Bend area. 

Simulated pumping of a deep well in Bend 

(model layer 7) shows how the response of the 

ground-water system can vary depending on well 

depth. After 2 years of pumping, most of the 

streamflow effects occur in the regional discharge 

area just above Lake Billy Chinook and in the 

area just above and below Benham Falls (fig. 44). 

The rate of streamflow depletion after 2 years is 

about 1.7 ft3/s, or 17 percent of the pumping rate.

0

2

4

6

8

10

0 10 20 30 40 50

TIME IN YEARS

CUBIC FEET PER SECOND

LAYER 2

From storage

Total decrease in streamflow

From decreased stream gain

From increased stream loss

Figure 44. Simulated drawdown in model layer 7 and net reduction in ground-water discharge to streams after 2 and 10 years of 

pumping 10 ft3/s from model layer 7 in the Bend area. 

121o30' 121o00'

44o30'

44o00'

121o30' 121o00'

44o30'

44o00'

10 11 12 13 14 15

20

19

18

17

16

14

13

12

11

15

10 11 12 13 14 15

20

19

18

17

16

14

13

12

11

15

Millican

Sunriver

Bend

Redmond

Cline

Falls

Madras

Sisters

Metolius

Culver

Terrebonne

Lower

Bridge

Millican

Sunriver

Bend

Redmond

Cline

Falls

Madras

Sisters

Metolius

Culver

Terrebonne

Lower

Bridge

Arnold

Canal

Tumalo Feed Canal

Central Oregon 

Canal

North Unit

Main Canal

Squaw Creek

Canal

Arnold

Canal

Tumalo Feed Canal

Central Oregon 

Canal

North Unit

Main Canal

Squaw Creek

Canal

Lava

Butte

Cline

Buttes

Juniper

Butte

Osborne

Canyon

Gray

Butte

Lower Desert

Lava

Butte

Cline

Buttes

Juniper

Butte

Osborne

Canyon

Gray

Butte

Lower Desert

Study area

 boundary

Study area

 boundary

1

3

5

10

20

Spring

Haystack

Reservoir

Creek

Mud Springs

 Creek

Lake Billy

Willow

Chinook

Lake

Simtustus

Creek

Fly

Crooked

River

Metolius

Round

Butte

Dam

Houston

Lake

Dry

Creek

Cr

River

Squaw

River

River

Cr

Tumalo

DESCHUTES

RIVER

Indian

Ford

Spring

Haystack

Reservoir

Creek

Mud Springs

 Creek

Lake Billy

Willow

Chinook

Lake

Simtustus

Creek

Fly

Crooked

River

Metolius

Round

Butte

Dam

Houston

Lake

Dry

Creek

Cr

River

Squaw

River

River

Cr

Tumalo

DESCHUTES

RIVER

Indian

Ford

Benham

Falls

Benham Falls

97

97

20

26

20

126

126

26

97

97

20

26

20

126

126

26

1

1

3

5

10

20

0 5 10 MILES

0 5 10 KILOMETERS

EXPLANATION

Streamflow reduction—In cubic feet per second per mile

 .01-.05 

 .05-.1

 More than .1

Line of equal drawdown—In feet. Contour intervals are 2, 5, and 10 feet.

Well location

10

2 YEARS 10 YEARS

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

 See table of contents for detailed information on sources.
71

After 10 years of pumping, effects have spread to a 

few more stream reaches, but most effects are still 

seen along the Deschutes and Crooked Rivers in 

the regional discharge area upstream of Lake Billy 

Chinook (fig. 44). Stream gains and losses in the Bend 

area are relatively unaffected by deep pumping. The 

rate of streamflow depletion after 10 years of pumping 

is about 5.7 ft3

/s, or 57 percent of the pumping rate. As 

in the case of the shallower well, the amounts of water 

coming from storage and from diminished streamflow 

are equal at about 7 years (fig. 45). About 90 percent 

of the pumped water comes from diminished streamflow after about 35 years, sooner than the shallower 

simulation because the smaller storage coefficients in 

deeper model layers cause the cone of depression to 

expand and capture streamflow more quickly. Maximum drawdown in model layer 7 is about 30 feet 

after 2 years and does not increase after 10 years, 

but the cone of depression is more widespread.

Figure 45. Variations with time in the amounts of ground water 

coming from storage and from diminished streamflow in response 

to pumping 10 ft3/s from model layer 7 in the Bend area. 

Simulated Pumping in the Redmond Area 

Simulations suggest that streamflow depletion 

caused by pumping in the Redmond area (row 63, 

column 67) is concentrated in the ground-water 

discharge areas along the Deschutes and Crooked 

Rivers upstream of Lake Billy Chinook. Because 

vertical anisotropy is less than in the Bend area, the 

ground-water system responds similarly to simulated 

pumping of either shallow or deep wells (figs. 46 and 

47). After 2 years of pumping, the rates of streamflow 

depletion are 1.6 and 1.9 ft3/s (16 and 19 percent) for 

the shallow and deep wells, respectively. Most of the 

depletion is due to diminished spring discharge just 

above Lake Billy Chinook (figs. 46 and 47). After 

10 years of pumping, streamflow depletion is about 

5.8 ft3/s (58 percent of the pumping rate) for both the 

shallow and deep simulations. Streamflow depletion is 

more widespread after 10 years, but most of the effect 

is still in the inflow area above Lake Billy Chinook. 

Both simulations in the Redmond area show that the 

amounts of ground water coming from storage and 

diminished streamflow are equal after about 7 years, 

and that 90 percent of the pumping rate comes from 

diminished streamflow after about 32 years (fig. 48). 

Simulated drawdowns are small in the Redmond area 

due to the high transmissivity; this is consistent with 

aquifer-test results. Simulated drawdown around the 

shallow well is only 4 feet after 10 years of pumping, 

and only 2 feet in the deep well after the same period 

(figs. 46 and 47).

0

2

4

6

8

10

0 10 20 30 40 50

TIME IN YEARS

CUBIC FEET PER SECOND

LAYER 7

From storage

Total decrease in streamflow

From decreased stream gain

From increased stream loss
72

Figure 46. Simulated drawdown and net reduction in ground-water discharge to streams after 2 and 10 years of pumping 10 ft3/s 

from model layer 2 in the Redmond area. 

44o30'

44o00'

121o30' 121o00'

44o30'

44o00'

121o30' 121o00'

10 11 12 13 14 15

20

19

18

17

16

14

13

12

11

15

10 11 12 13 14 15

20

19

18

17

16

14

13

12

11

15

Millican

Sunriver

Bend

Redmond

Cline

Falls

Madras

Sisters

Metolius

Culver

Terrebonne

Lower

Bridge

Millican

Sunriver

Bend

Redmond

Cline

Falls

Madras

Sisters

Metolius

Culver

Terrebonne

Lower

Bridge

Arnold

Canal

Tumalo Feed Canal

Central Oregon 

Canal

North Uni

Main Canal

Squaw Creek

Canal

t

Arnold

Canal

Tumalo Feed Canal

Central Oregon 

Canal

North Uni

Main Canal

Squaw Creek

Canal

t

Lava

Butte

Cline

Buttes

Juniper

Butte

Osborne

Canyon

Gray

Butte

Lower Desert

Lava

Butte

Cline

Buttes

Juniper

Butte

Osborne

Canyon

Gray

Butte

Lower Desert

Spring

Haystack

Reservoir

Creek

Mud Springs

 Creek

Lake Billy

Willow

Chinook

Lake

Simtustus

Creek

Fly

Crooked

River

Metolius

Round

Butte

Dam

Houston

Lake

Dry

Creek

Cr

River

Squaw

River

River

Cr

Tumalo

DESCHUTES

RIVER

Spring

Haystack

Reservoir

Creek

Mud Springs

 Creek

Lake Billy

Willow

Chinook

Lake

Simtustus

Creek

Fly

Crooked

River

Metolius

Round

Butte

Dam

Houston

Lake

Dry

Creek

Cr

River

Squaw

River

River

Cr

Tumalo

DESCHUTES

RIVER

Indian

Ford

Indian

Ford

97

97

20

26

20

126

126

26

97

97

20

26

20

126

126

26

3

1

3

1

1

0 5 10 MILES

0 5 10 KILOMETERS

EXPLANATION

Streamflow reduction—In cubic feet per second per mile.

 .01-.05 

 .05-.1

 More than .1

Line of equal drawdown—In feet. Contour intervals are 2, 5, and 10 feet.

Well location

10

2 YEARS 10 YEARS

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

 See table of contents for detailed information on sources.
73

Figure 47. Simulated drawdown and net reduction in ground-water discharge to streams after 2 and 10 years of pumping 10 ft3/s 

from model layer 7 in the Redmond area.

Figure 48. Variations with time in the amounts of ground water coming from storage and from diminished streamflow in response to 

pumping 10 ft3/s from model layers 2 and 7 in the Redmond area.

44o30'

44o00'

121o30' 121o00' 121o30' 121o00'

44o30'

44o00'

10 11 12 13 14 15

20

19

18

17

16

14

13

12

11

15

10 11 12 13 14 15

20

19

18

17

16

14

13

12

11

15

Millican

Sunriver

Bend

Redmond

Cline

Falls

Madras

Sisters

Metolius

Culver

Terrebonne

Lower

Bridge

Millican

Sunriver

Bend

Redmond

Cline

Falls

Madras

Sisters

Metolius

Culver

Terrebonne

Lower

Bridge

Arnold

Canal

Tumalo Feed Canal

Central Oregon 

Canal

North Unit

Main Canal

Squaw Creek

Canal

Arnold

Canal

Tumalo Feed Canal

Central Oregon 

Canal

North Unit

Main Canal

Squaw Creek

Canal

Lava

Butte

Cline

Buttes

Juniper

Butte

Osborne

Canyon

Gray

Butte

Lower Desert

Lava

Butte

Cline

Buttes

Juniper

Butte

Osborne

Canyon

Gray

Butte

Lower Desert

Spring

Haystack

Reservoir

Creek

Mud Springs

 Creek

Lake Billy

Willow

Chinook

Lake

Simtustus

Creek

Fly

Crooked

River

Metolius

Round

Butte

Dam

Houston

Lake

Dry

Creek

Cr

River

Squaw

River

River

Cr

Tumalo

DESCHUTES

RIVER

Spring

Haystack

Reservoir

Creek

Mud Springs

 Creek

Lake Billy

Willow

Chinook

Lake

Simtustus

Creek

Fly

Crooked

River

Metolius

Round

Butte

Dam

Houston

Lake

Dry

Creek

Cr

River

Squaw

River

River

Cr

Tumalo

DESCHUTES

RIVER

Indian

Ford

Indian

Ford

97

97

20

26

20

126

126

26

97

97

20

26

20

126

126

26

1

1

1

0 5 10 MILES

0 5 10 KILOMETERS

EXPLANATION

Streamflow reduction—In cubic feet per second per mile.

 .01-.05 

 .05-.1

 More than .1

Line of equal drawdown—In feet. Contour intervals are 2, 5, and 10 feet.

Well location

10

2 YEARS 10 YEARS

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

 See table of contents for detailed information on sources.

0

2

4

6

8

10

CUBIC FEET PER SECOND

From storage

Total decrease in streamflow

From decreased stream gain

From increased stream loss

LAYER 2

0 10 20 30 40 50

TIME IN YEARS

0

2

4

6

8

10

0 10 20 30 40 50

TIME IN YEARS

CUBIC FEET PER SECOND

LAYER 7
74

Simulated Pumping between the Deschutes 

and Crooked Rivers near Osborne Canyon

The last well simulations are at an arbitrary location between the Deschutes and Crooked Rivers near 

their confluence (row 37, column 53) and are intended to 

demonstrate the effects of pumping in close proximity to 

a major ground-water discharge area. After 2 years of 

pumping from model layer 2, about 7.3 ft3

/s (73 percent 

of the pumping rate) of the pumping comes from streamflow (fig. 49). Effects after 2 years are limited to the 

lower Crooked and the Deschutes Rivers just above Lake 

Billy Chinook and to lower Squaw Creek. After 10 years 

of pumping from model layer 2, about 8.5 ft3

/s (85 percent of the pumping rate) comes from streamflow 

depletion, predominantly from the same stream reaches. 

Simulation of pumping from model layer 2 at this location indicates that the amounts of water coming from 

storage and from diminished streamflow are equal after 

only about 2 weeks and, that 90 percent of the pumping 

comes from streamflow depletion after about 12 years 

(fig. 50). Maximum drawdown in model layer 2 is 

about 10 feet after 2 years of pumping and does not 

significantly increase by 10 years (fig. 49). Maximum 

drawdown in model layer 1 is about 2 feet.

Figure 49. Simulated drawdown in model layer 2 and net reduction in ground-water discharge to streams after 2 and 10 years of 

pumping 10 ft3/s from model layer 2 between the Deschutes and Crooked Rivers near Osborne Canyon.

121o30' 121o00'

44o30'

44o00'

121o30' 121o00'

44o30'

44o00'

10 11 12 13 14 15

20

19

18

17

16

14

13

12

11

15

10 11 12 13 14 15

20

19

18

17

16

14

13

12

11

15

Millican

Sunriver

Bend

Redmond

Cline

Falls

Madras

Sisters

Metolius

Culver

Terrebonne

Lower

Bridge

Millican

Sunriver

Bend

Redmond

Cline

Falls

Madras

Sisters

Metolius

Culver

Terrebonne

Lower

Bridge

Arnold

Canal

Tumalo Feed Canal

Central Oregon 

Canal

North Unit

Main Canal

Squaw Creek

Canal

Arnold

Canal

Tumalo Feed Canal

Central Oregon 

Canal

North Unit

Main Canal

Squaw Creek

Canal

Lava

Butte

Cline

Buttes

Juniper

Butte

Osborne

Canyon

Gray

Butte

Lower Desert

Lava

Butte

Cline

Buttes

Juniper

Butte

Osborne

Canyon

Gray

Butte

Lower Desert

Study area

 boundary

Study area

 boundary

Spring

Haystack

Reservoir

Creek

Mud Springs

 Creek

Lake Billy

Willow

Chinook

Lake

Simtustus

Creek

Fly

Crooked

River

Metolius

Round

Butte

Dam

Houston

Lake

Dry

Creek

Cr

River

Squaw

River

River

Cr

Tumalo

DESCHUTES

RIVER

Spring

Haystack

Reservoir

Creek

Mud Springs

 Creek

Lake Billy

Willow

Chinook

Lake

Simtustus

Creek

Fly

Crooked

River

Metolius

Round

Butte

Dam

Houston

Lake

Dry

Creek

Cr

River

Squaw

River

River

Cr

Tumalo

DESCHUTES

RIVER

Indian

Ford

Indian

Ford

97

97

20

26

20

126

126

26

97

97

20

26

20

126

126

26

1

3

5 5

3

1

0 5 10 MILES

0 5 10 KILOMETERS

EXPLANATION

Streamflow reduction—In cubic feet per second per mile.

 .01-.05 

 .05-.1

 More than .1

Line of equal drawdown—In feet. Contour intervals are 2, 5, and 10 feet.

Well location

10

2 YEARS 10 YEARS

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

 See table of contents for detailed information on sources.
75

Figure 50. Variations with time in the amounts of ground water 

coming from storage and from diminished streamflow in response 

to pumping 10 ft3/s from model layer 2 between the Deschutes and 

Crooked Rivers near Osborne Canyon. 

The response of streamflow to simulated pumping from model layer 7 near the confluence is only 

slightly slower then the response to pumping from 

model layer 2. After 2 years of pumping from model 

layer 7, about 4.7 ft3

/s (47 percent of the pumping 

rate) comes from diminished streamflow (fig. 51). The 

affected stream reaches are virtually identical to those 

affected during the model layer 2 simulation, lower 

Squaw Creek and the lower Crooked River and Deschutes River just above Lake Billy Chinook (fig. 51). 

After 10 years of pumping, about 7.2 ft3/s or 72 percent of the pumping rate is supplied by diminished 

streamflow. 

0

2

4

6

8

10

0 10 20 30 40 50

TIME IN YEARS

CUBIC FEET PER SECOND

LAYER 2

From storage

Total decrease in streamflow

From decreased stream gain

From increased stream loss

Figure 51. Simulated drawdown in model layer 7 and net reduction in ground-water discharge to streams after 2 and 10 years of pumping 

10 ft3/s from model layer 7 between the Deschutes and Crooked Rivers near Osborne Canyon. 

121o30' 121o00'

44o30'

44o00'

121o30' 121o00'

44o30'

44o00'

10 11 12 13 14 15

20

19

18

17

16

14

13

12

11

15

10 11 12 13 14 15

20

19

18

17

16

14

13

12

11

15

Millican

Sunriver

Bend

Redmond

Cline

Falls

Madras

Sisters

Metolius

Culver

Terrebonne

Lower

Bridge

Millican

Sunriver

Bend

Redmond

Cline

Falls

Madras

Sisters

Metolius

Culver

Terrebonne

Lower

Bridge

Arnold

Canal

Tumalo Feed Canal

Central Oregon 

Canal

North Unit

Main Canal

Squaw Creek

Canal

Arnold

Canal

Tumalo Feed Canal

Central Oregon 

Canal

North Unit

Main Canal

Squaw Creek

Canal

Lava

Butte

Cline

Buttes

Juniper

Butte

Osborne

Canyon

Gray

Butte

Lower Desert

Lava

Butte

Cline

Buttes

Juniper

Butte

Osborne

Canyon

Gray

Butte

Lower Desert

Study area

 boundary

Study area

 boundary

Spring

Haystack

Reservoir

Creek

Mud Springs

 Creek

Lake Billy

Willow

Chinook

Lake

Simtustus

Creek

Fly

Crooked

River

Metolius

Round

Butte

Dam

Houston

Lake

Dry

Creek

Cr

River

Squaw

River

River

Cr

Tumalo

DESCHUTES

RIVER

Spring

Haystack

Reservoir

Creek

Mud Springs

 Creek

Lake Billy

Willow

Chinook

Lake

Simtustus

Creek

Fly

Crooked

River

Metolius

Round

Butte

Dam

Houston

Lake

Dry

Creek

Cr

River

Squaw

River

River

Cr

Tumalo

DESCHUTES

RIVER

Indian

Ford

Indian

Ford

97

97

20

26

20

126

126

26

97

97

20

26

20

126

126

26

5

1

3

1

5

3

0 5 10 MILES

0 5 10 KILOMETERS

EXPLANATION

Streamflow reduction—In cubic feet per second per mile.

 .01-.05 

 .05-.1

 More than .1

Line of equal drawdown—In feet. Contour intervals are 2, 5, and 10 feet.

Well location

10

2 YEARS 10 YEARS

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

 See table of contents for detailed information on sources.
76

By 10 years, simulated effects expand up the Deschutes River beyond Bend and to the lower Metolius 

River (fig. 51), but are still concentrated along the 

Deschutes and Crooked Rivers and lower Squaw 

Creek just above Lake Billy Chinook. The amounts 

of water from storage and streamflow depletion are 

equal after about 2 years, and 90 percent of the pumping rate is supplied from streamflow depletion after 

about 21 years (fig. 52). After 2 years, maximum 

drawdown is about 7 feet in model layer 7 and less 

than 1 foot in model layer 1. 

Figure 52. Variations with time in the amounts of ground water 

coming from storage and from diminished streamflow in response 

to pumping 10 ft3/s from model layer 7 between the Deschutes and 

Crooked Rivers near Osborne Canyon. 

Simulated Effects of Canal Lining

Lining canals to reduce leakage affects the 

ground-water system by locally reducing artificial 

recharge. As a result, the water-table elevation drops 

in the vicinity of the lined canal and, eventually, 

streamflow diminishes by the amount artificial 

recharge is reduced. As with pumping, streamflow 

does not respond immediately to canal lining. Streamflow responds gradually as the water table drops and 

approaches a new equilibrium. As the water table 

drops, a certain amount of the reduced recharge is 

accommodated for by changes in aquifer storage. 

Two canal-lining scenarios are presented in this 

section of the report.

In the first scenario, leakage in all main canals 

in the vicinity of Bend was reduced by 50 percent 

(fig. 53). This results in an average annual reduction in 

leakage of 69 ft3/s (based on 1994 leakage estimates). 

After 2 years, net streamflow reduction is about 

30 ft3/s, mostly occurring in the vicinity of Bend. 

The changes in streamflow include both increased 

stream leakage and diminished stream gains. Gagingstation data discussed in preceding sections indicate 

that both gaining and losing stream reaches occur in 

the vicinity of Bend; this result is reasonable and 

expected. After 10 years, streamflow depletion is 

about 47 ft3/s, and the effects expand to the major 

ground-water discharge area just above Lake Billy 

Chinook. After roughly 3 years, about one-half of 

the reduced canal leakage is reflected as diminished 

streamflow (fig. 54); the other one-half is accommodated by changes in aquifer storage. After 31 years, 

streamflow depletion reflects about 90 percent of the 

reduced canal-leakage rate. Maximum head changes 

due to the simulated canal lining are as much as 60 feet 

near the canal diversions north of Bend, but range 

between 5 and 30 feet over a larger area (fig. 53).

0

2

4

6

8

10

0 10 20 30 40 50

TIME IN YEARS

CUBIC FEET PER SECOND

LAYER 7

From storage

Total decrease in streamflow

From decreased stream gain

From increased stream loss
77

121o30' 121o00'

30'

15'

44o45'

44o00'

121o30' 121o00'

30'

15'

44o45'

44o00'

10 11 12 13 14 15

20

19

18

17

16

14

13

12

11

15

10 11 12 13 14 15

20

19

18

17

16

14

13

12

11

15

Millican

Sunriver

Bend

Redmond

Cline

Falls

Madras

Sisters

Metolius

Culver

Terrebonne

Lower

Bridge

Millican

Sunriver

Bend

Redmond

Cline

Falls

Madras

Sisters

Metolius

Culver

Terrebonne

Lower

Bridge

Arnold

Canal

Tumalo Feed Canal

Central Oregon 

Canal

North Unit

Main Canal

Squaw Creek

Canal

Arnold

Canal

Tumalo Feed Canal

Central Oregon 

Canal

North Unit

Main Canal

Squaw Creek

Canal

Lava

Butte

Cline

Buttes

Juniper

Butte

Osborne

Canyon

Gray

Butte

Lower Desert

Lava

Butte

Cline

Buttes

Juniper

Butte

Osborne

Canyon

Gray

Butte

Lower Desert

Study area

 boundary

Study area

 boundary

Spring

Haystack

Reservoir

Creek

Mud Springs

 Creek

Lake Billy

Willow

Chinook

Lake

Simtustus

Creek

Fly

Crooked

River

Metolius

Round

Butte

Dam

Houston

Lake

Dry

Creek

Cr

River

Squaw

River

River

Cr

Tumalo

DESCHUTES

RIVER

Pelton

Dam

Spring

Haystack

Reservoir

Creek

Mud Springs

 Creek

Lake Billy

Willow

Chinook

Lake

Simtustus

Creek

Fly

Crooked

River

Metolius

Round

Butte

Dam

Houston

Lake

Dry

Creek

Cr

River

Squaw

River

River

Cr

Tumalo

DESCHUTES

RIVER

Pelton

Dam

Indian

Ford

Indian

Ford

97

97

20

26

20

126

126

26

97

97

20

26

20

126

126

26

1

5

10

20

30

40 50

60

1

1

5

10

30

40

1

5

20

10

20

30

40

50

0 5 10 MILES

0 5 10 KILOMETERS

10

EXPLANATION

Streamflow reduction—In cubic feet per second per mile.

 .05-.1

 .1-.5

 .5-1.0

 More than 1.0

Line of equal head drop—In feet. Contour intervals are 4, 5, and 10 feet.

Canal with leakage reduced—69 cubic feet per second reduction of 

 annual mean leakage from canal

2 YEARS 10 YEARS

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

 See table of contents for detailed information on sources.

Figure 53. Simulated head decline in model layer 1 and diminished streamflow due to a 50 percent reduction in leakage of main 

irrigation canals in the vicinity of Bend after 2 and 10 years. 

Figure 54. Variations with time in the amounts of ground 

water coming from storage and from diminished streamflow in response to a 50 percent reduction in leakage of 

0 main irrigation canals in the vicinity of Bend. 

10

20

30

40

50

60

70

0 10 20 30 40 50

TIME IN YEARS

CUBIC FEET PER SECOND

From storage

Total decrease in streamflow

From decreased stream gain

From increased stream loss
78

The second canal-lining scenario involved 

reducing the leakage in the northern sections of main 

canals (between Redmond and Madras) by 50 percent 

(fig. 55). This results in an average annual reduction in 

canal leakage of about 19 ft3

/s (based on 1994 leakage 

estimates). The response of streamflow to the reduction in recharge is slower than the previous canallining scenario. After 2 years, the net reduction in 

streamflow due to the reduced recharge is only about 

2.1 ft3/s, or about 11 percent of the reduced recharge. 

Ground-water discharge is affected along the lower 

Crooked River above Osborne Canyon, the Deschutes 

River between Lower Bridge and Squaw Creek, and at 

Mud Springs Creek (fig. 55). After 10 years, about 

8.2 ft3/s, or about 43 percent of the reduced recharge, 

is reflected as diminished streamflow. Similar stream 

reaches are affected, although the effects spread 

to include Willow Creek west of Madras (fig. 55). 

Figure 55. Simulated head decline in model layer 1 and diminished streamflow due to a 50 percent reduction in leakage of main 

irrigation canals north of Redmond after 2 and 10 years. 

121o30' 121o00'

30'

15'

44o45'

44o00'

121o30' 121o00'

30'

15'

44o45'

44o00'

10 11 12 13 14 15

20

19

18

17

16

14

13

12

11

15

10 11 12 13 14 15

20

19

18

17

16

14

13

12

11

15

Millican

Sunriver

Bend

Redmond

Cline

Falls

Madras

Sisters

Metolius

Culver

Terrebonne

Lower

Bridge

Millican

Sunriver

Bend

Redmond

Cline

Falls

Madras

Sisters

Metolius

Culver

Terrebonne

Lower

Bridge

Arnold

Canal

Tumalo Feed Canal

Central Oregon 

Canal

North Unit

Main Canal

Squaw Creek

Canal

Arnold

Canal

Tumalo Feed Canal

Central Oregon 

Canal

North Unit

Main Canal

Squaw Creek

Canal

Lava

Butte

Cline

Buttes

Juniper

Butte

Osborne

Canyon

Gray

Butte

Lower Desert

Lava

Butte

Cline

Buttes

Juniper

Butte

Osborne

Canyon

Gray

Butte

Lower Desert

Study area

 boundary

Study area

 boundary

Spring

Haystack

Reservoir

Creek

Mud Springs

 

Willow

 Creek

Lake Billy

Chinook

Lake

Simtustus

Creek

Fly

Crooked

River

Metolius

Round

Butte

Dam

Houston

Lake

Dry

Creek

Cr

River

Squaw

River

River

Cr

Tumalo

DESCHUTES

RIVER

Pelton

Dam

Spring

Haystack

Reservoir

Creek

Mud Springs

 

Willow

 Creek

Lake Billy

Chinook

Lake

Simtustus

Creek

Fly

Crooked

River

Metolius

Round

Butte

Dam

Houston

Lake

Dry

Creek

Cr

River

Squaw

River

River

Cr

Tumalo

DESCHUTES

RIVER

Pelton

Dam

Indian

Ford

Indian

Ford

97

97

20

26

20

126

126

26

97

97

20

26

20

126

126

26

1

5

10

10

10

10

1

5

50

50

50

20

30 40

60

90

1

5

10

5

1

20

1

5

1

5

5

5

10

1

5 10

30

30

1

5

1

40

10

10

20

0 5 10 MILES

0 5 10 KILOMETERS

EXPLANATION

Streamflow reduction—In cubic feet per second per mile.

 .05-.1

 .1-.5

 .5-1.0

 More than 1.0

Line of equal head decline—In feet. Contour intervals are 4, 5, and 10 feet.

Canal with leakage reduced—19 cubic feet per second reduction of 

 annual mean leakage from canal

10

2 YEARS 10 YEARS

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

 See table of contents for detailed information on sources.
79

It takes 13 years until one-half of the reduction in 

canal loss is reflected as diminished streamflow 

(fig. 56). Only 80 percent of the loss rate reduction is 

reflected as diminished streamflow after 50 years, the 

end of the simulation. The increased time to equilibrium in this simulation compared to the canal-lining 

simulation near Bend is due to the larger storage coefficients in the northern part of the model area and the 

larger distance to the affected stream reaches. Maximum head changes are only about 5 to 10 feet and 

occur mostly along the canals (fig. 55). Head changes 

of up to 10 feet are localized along the North Unit 

canal in the Madras area. Head changes of 5 to 10 feet 

are widespread east of Redmond, between the North 

Unit and Central Oregon Canals. 

Figure 56. Variations with time in the amounts of ground water 

coming from storage and from diminished streamflow in response 

to a 50 percent reduction in leakage of main irrigation canals north 

of Redmond. 

The simulations in this section show how the 

upper Deschutes Basin ground-water model can be 

used. The ground-water model could also be used to 

evaluate other types of pumping scenarios, such as 

with multiple wells or with intermittent or varying 

pumping schedules. Other possible canal scenarios 

could include different locations for lining or the use 

of canals for artificial ground-water recharge. Another 

potential use of the model would be to evaluate the 

effects of long-term climate change on the regional 

ground-water system. 

SUMMARY AND CONCLUSIONS

Rapid population growth is increasing demand 

for water in the upper Deschutes Basin. Surface-water 

resources are fully allocated, so the increased demand 

is being supplied chiefly by ground water. Numerous 

hydrologic investigators over the past 100 years have 

recognized that ground water and surface water are 

intimately connected in the upper Deschutes Basin 

and that most of the flow of the Deschutes River is 

supplied by ground water (Russell, 1905; Stearns, 

1931; Sceva, 1960, 1968; Gannett and others, 2001). 

Because of this connection, pumping and consumptively using ground water in the upper Deschutes 

Basin reduces the amount of ground water discharging 

to streams and, consequently, streamflow. Reducing 

artificial recharge by lining leaking irrigation canals 

will affect streams in a similar manner. Until now, 

regulatory agencies and the public have had no means 

to evaluate the effects of pumping, canal lining, or 

other activities on the ground-water and surface-water 

systems. The model described in this report provides 

a means for evaluating the possible effects of such 

activities, particularly with regard to water-level 

changes and the timing and distribution of the effects 

on streamflow.

The model covers most of the Deschutes Basin 

upstream of the mouth of Trout Creek north of 

Madras. Recharge is principally from precipitation in 

the Cascade Range, with a lesser amount of artificial 

recharge from leaking irrigation canals. Ground water 

moves from the Cascade Range east and north and 

ultimately discharges to streams in the basin. There are 

two principal settings where ground water discharges: 

to streams along the topographic margin of the 

Cascade Range, and to the Deschutes, Crooked, and 

Metolius Rivers near their confluence in the vicinity 

of Lake Billy Chinook. Approximately one-half of 

the ground-water discharge occurs at each of these 

settings. Geology is the principal factor controlling 

the discharge of ground water in the confluence area. 

It is not far north of the confluence area where the permeable Deschutes Formation strata pinch out against 

uplands of the low-permeability John Day Formation. 

The northward-flowing ground water hits this low 

permeability boundary and discharges to the deeply 

incised streams. 

Ground-water flow is simulated by the finitedifference method using the programs MODFLOW 

(McDonald and Harbaugh, 1988) and MODFLOWP, 

a version of MODFLOW modified by Hill (1992) 

to estimate parameters by nonlinear regression. The 

finite-difference grid consists of 8 layers, 127 rows, 

and 87 columns. Geographic boundaries are simulated 

as no-flow except in two areas where there is inferred 

subsurface inflow across topographic divides. All 

major streams and most principal tributaries in the 

upper Deschutes Basin are included in the model 

as head-dependent flux boundaries. Ground-water 

0

5

10

15

20

0 10 20 30 40 50

TIME IN YEARS

CUBIC FEET PER SECOND

From storage

Total decrease in streamflow

From decreased stream gain

From increased stream loss
80

recharge from precipitation was estimated using the 

Deep Percolation Model of Bauer and Vaccaro (1987). 

Artificial recharge from leaking irrigation canals 

and on-farm losses was estimated from diversion 

and delivery records, seepage studies, and crop data. 

Details of recharge estimation can be found in Boyd 

(1996) and Gannett and others (2001). Evapotranspiration from the water table in areas of shallow depths 

to water is simulated as a head-dependent process. 

Ground-water pumpage for irrigation and public water 

supplies is included in the model. Irrigation pumpage 

was estimated using information on water rights and 

crop data. Public-supply pumpage was estimated from 

pumping records provided by water suppliers and 

from population data. Hydraulic conductivities were 

estimated from well-yield tests, aquifer tests, and 

literature values. Data are too sparse and the geology 

too heterogeneous to allow continuous cell-by-cell 

mapping of the hydraulic-conductivity distribution. 

The distribution of hydraulic conductivity is represented by a set of discrete zones within which conductivity values were uniform. The zonation scheme is 

based on hydrogeologic units defined by Lite and 

Gannett (2002), which were further subdivided during 

model calibration.

The model was calibrated to time-averaged 

steady-state conditions from 1993 to 1995 using 

parameter-estimation programs and techniques of 

Hill (1992, 1998). These methods involve determining 

the set of optimal parameter values that best fit the 

observed hydraulic heads and flows by using nonlinear regression to minimize a weighted sum-of-squares 

objective function. Fourteen hydraulic-conductivity 

parameters and two vertical conductance parameters 

were determined using nonlinear regression. Other 

parameter values mostly had sensitivities too low to 

allow optimization. Parameter values not determined 

by regression were set to reasonable and expected 

values and adjusted judiciously by trial and error during model calibration. Final parameter values are all 

within reasonable and expected ranges. 

The steady-state model fit to measured heads 

and flows is good. Residuals are reasonably random 

and normally distributed. The general shape and slope 

of the simulated water-table surface and overall 

hydraulic-head distribution match the geometry determined from field measurements. The fitted standard 

deviation for hydraulic head is about 76 feet, which is 

considered good given that heads in the upper basin 

span a range of over 4,500 feet and field measurements range from 1,797 to 5,586 feet. The only area 

of significant systematic error in simulated heads is 

in the La Pine subbasin, where simulated heads are 

too high. The general magnitude and distribution of 

ground-water discharge to streams is also reasonably 

well simulated throughout the model. Ground-water 

discharge to streams in the area of the confluence 

of the Deschutes, Crooked, and Metolius Rivers is 

closely matched. Significant differences between 

simulated and measured ground-water discharge to 

streams, where they occur, are generally attributable to 

scale problems where highly heterogeneous hydrogeologic conditions occur in areas where subsurface data 

are sparse and the model grid is coarse. Systematic 

errors in ground-water discharge to streams occur in 

high-elevation tributary streams in the Cascade Range, 

where simulated discharge is less than observed. 

Overall, the steady-state model reasonably reproduces 

the behavior of the regional ground-water system.

The model was calibrated to transient conditions 

from 1978 to 1997 (using 40 semiannual stress periods) by traditional trial-and-error methods. Climatic 

cycles during this period provided an excellent regional hydrologic signal for calibration. Model parameterization is the same as used in the steady-state 

calibration, with the addition of storage coefficient 

parameters that were set up to follow the hydraulicconductivity zonation scheme. Unlike hydraulic 

conductivity, storage coefficients were varied from 

model layer to model layer, usually decreasing with 

depth. Final storage coefficient values are all within 

reasonable and expected ranges, and are consistent 

with aquifer-test results. Other model parameters are 

the same as for the steady-state calibration.

Simulated variations in ground-water discharge 

to streams were evaluated by comparing simulated 

values and measured values in six areas. Overall, the 

transient model simulated the volumetric distribution 

and temporal variations in ground-water discharge 

reasonably well. The match between simulated and 

measured volume of and variations in ground-water 

discharge was, however, somewhat dependent on 

geographic scale. Simulated discharge fluxes and 

variations matched observed conditions best at 

regional scales. At smaller scales (as would be represented by first-order tributaries) geologic heterogeneities, topographic complexity, and model 

discretization combine to reduce the accuracy of 

matches between observed and simulated discharge 

fluxes and variations. Nevertheless, simulated fluxes 
81

and variations are very close for certain small streams, 

such as Odell Creek. The overall magnitude and temporal variations of stream discharge in the confluence 

area are reproduced by the transient model with reasonable accuracy.

Simulated variations in hydraulic head were 

evaluated by comparing simulated and measured 

water-level fluctuations in 64 wells in 12 subareas. 

The transient model simulates fluctuations caused 

by climate cycles, annual recharge cycles, and canal 

leakage and irrigation. Water-level fluctuations caused 

by stream-stage variations are not simulated. Climatedriven water-level fluctuations are simulated with 

reasonable accuracy over most of the model area. 

The timing of the simulated water-level response to 

the change to wet conditions starting in 1996 matches 

the observed response within months in most cases 

and generally within a year. The timing and magnitude 

of simulated water-level fluctuations caused by annual 

pulses of recharge from precipitation match those 

observed reasonably well, given the limitations of the 

time discretization in the model. Water-level fluctuations caused by annual canal leakage are simulated 

very well over most of the area where this occurs. 

Example simulations were run to demonstrate 

the model usefulness for predicting the effects of 

changes to the hydrologic stresses on the system. The 

examples involve increased ground-water pumping 

or canal lining. All example pumping simulations 

assumed continuous pumping, but this is not required. 

All pumping simulations show that pumped water 

comes largely from aquifer storage when pumping 

begins, but that as the water table stabilizes the 

pumping increasingly diminishes the discharge to 

streams and, hence, streamflow. The time it takes for 

pumping to affect streamflow varies geographically. 

Simulations involving continuous pumping of wells 

near Bend and Redmond show that, in general, about 

one-half of the pumping comes from diminished 

streamflow after 5 to 10 years and 90 percent of the 

pumping comes from diminished streamflow after 

30 to 40 years. Pumping in the Bend and Redmond 

areas generally affects streamflow in the regional 

discharge area near the confluence of the Deschutes, 

Crooked, and Metolius Rivers. However, shallow 

wells in the Bend area tend to affect streamflow in and 

around Bend. Wells placed very close to streams have 

a more rapid effect because of the lack of intervening 

aquifer storage. Simulated pumping of a well between 

the Crooked and Deschutes Rivers at about the latitude 

of Osborne Canyon showed that 50 percent of the 

pumping rate was supplied by diminished streamflow 

after a few weeks.

Two example simulations were run to demonstrate the model’s ability to estimate the effects of 

canal lining and consequent reduction in artificial 

recharge. The simulations show that canal lining has 

effects similar to pumping wells. The reduction in artificial recharge is at first offset by releases from aquifer 

storage as the water table drops. With time, however, 

ground-water discharge to streams diminishes in 

response to the reduced recharge and head changes. 

Simulated canal lining in the Bend area affected 

streamflow in much the same way as pumping shallow 

wells in the area, with the effects most prominent 

along the Deschutes River in the vicinity of Bend. 

Simulated ground-water levels dropped tens of feet 

in the Bend area in response to this canal-lining 

scenario. A second canal-lining simulation in the area 

from Redmond to Madras showed that one-half of the 

amount of reduced recharge is supplied by diminished 

streamflow after about 12 years, and that streamflow is 

diminished by about 80 percent of the amount after 50 

years (the length of the simulation). In this simulation, 

the stream reaches most affected include the lower 

Crooked River, the Deschutes River between Lower 

Bridge and Squaw Creek, and Mud Springs Creek.

The model described in this report simulates 

the behavior of the regional ground-water flow system 

in the upper Deschutes Basin and its response to new 

conditions. Examples in this report show how simulations can be used to understand the effects of new 

pumping or reduced recharge due to canal lining. 

The model could also be used to evaluate the effects 

of extended drought or climate change.

The model is intended to simulate the regional 

ground-water flow system and is most useful at 

regional scales. The model accurately simulated 

hydrologic conditions in specific local areas (usually 

corresponding to areas where data were plentiful); 

however, there are also areas where the model does 

not accurately simulate local or small-scale conditions. 

Simulation results must be evaluated within the 

context of overall hydrologic knowledge of the area. 

All model predictions have inherent error and uncertainty that is difficult to quantify. Nonetheless, this 

model should be a useful tool for resource managers 

and others for evaluating the response of the regional 

ground-water flow system to a range of management 

strategies and conditions. 
82
83

REFERENCES CITED

Anderson, M.P., and Woessner, W.W., 1992, Applied 

ground-water modeling: Academic Press Inc., 

San Diego, California, 381 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.

Blackwell, D.D., and Priest, G.R., 1996, Comment on “Rates 

and patterns of groundwater flow in the Cascade Range 

volcanic arc and the effect on subsurface temperatures” 

by S.E. Ingebritsen, D.R. Sherrod, and R.H. Mariner: 

Journal of Geophysical Research, v. 101, no. B8, 

p. 17561–17568.

Boyd, T.G., 1996, Groundwater recharge of the middle 

Deschutes Basin, Oregon: M.S. thesis, Portland State 

University, 86 p.

Caldwell, R.R., and Truini, Margot, 1997, Ground-water and 

water-chemistry data for the upper Deschutes Basin, 

Oregon: U.S. Geological Survey Open-File Report 

97–197, 77 p.

Century West Engineering Corporation, 1982, La Pine 

aquifer management plan: Bend, Oregon, prepared 

for Deschutes County, Oregon, variously paginated, 

several appendices.

Cooley, R.L., and Naff, R.L., 1990, Regression modeling 

of ground-water flow: U.S. Geological Survey 

Techniques of Water-Resources Investigations, 

book 3, chap. B4, 232 p.

Draper, N.R., and Smith, Harry, 1998, Applied regression 

analysis (3d ed.): John Wiley and Sons, New York, 

New York, 706 p.

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

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

Gannett, M.W., Lite, K.E., Jr., Morgan, D.S., and Collins, 

C.A., 2001, Ground-water hydrology of the upper 

Deschutes Basin, Oregon: U.S. Geological Survey 

Water-Resources Investigations Report 00–4162, 77 p.

[http://or.water.usgs.gov/pubs/WRIR00-4162/]

Gonthier, J.B., 1985, A description of aquifer units in 

eastern Oregon: U.S. Geological Survey WaterResources Investigations Report 84–4095, 39 p.

Hansen, A.J., Vaccaro, J.J., and Bauer, H.H., 1994, 

Ground-water 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.

Harbaugh, A.W., 1990, A simple contouring program for 

gridded data: U.S. Geological Survey Open-File Report 

90–144, 37 p.

Hill, M.C., 1990, Preconditioned conjugate-gradient 2 

(PCG2), a computer program for solving groundwater flow equations: U.S. Geological Survey WaterResources Investigations Report 90–4048, 43 p.

Hill, M.C., 1992, A computer program (MODFLOWP) for 

estimating parameters of a transient, three-dimensional, 

ground-water flow model using nonlinear regression: 

U.S. Geological Survey Open-File Report 91–484, 

358 p.

Hill, M.C., 1994, Five computer programs for testing 

weighted residuals and calculating linear confidence 

and prediction intervals on results from the groundwater parameter-estimation program MODFLOWP: 

U.S. Geological Survey Open-File Report 93–481, 

81 p.

Hill, M.C., 1998, Methods and guidelines for effective 

model calibration: U.S. Geological Survey WaterResources Investigations Report 98–4005, 90 p.

Hill, M.C., Cooley, R.L., and Pollock, D.W., 1998, A 

controlled experiment in ground-water flow model 

calibration using nonlinear regression: Ground Water, 

v. 36, p. 520–535.

Ingebritsen, S.E., Sherrod, D.R., and Mariner, R.H., 

1992, Rates and patterns of groundwater flow in 

the Cascade Range volcanic arc, and the effect on 

subsurface temperatures: Journal of Geophysical 

Research, v. 97, no. B4, p. 4599–4627.

Lite, K.E., Jr., and Gannett, M.W., 2002, Geologic framework of the regional ground-water flow system in the 

upper Deschutes Basin Oregon: U.S. Geological Survey Water-Resources Investigations Report 02–4015, 

44 p. [http://or.water.usgs.gov/pubs/WRIR02-4015/]

Mac Nish, R.D., and Barker, R.A., 1976, Digital simulation 

of a basalt aquifer system, Walla Walla River Basin, 

Washington and Oregon: State of Washington 

Department of Ecology Water-Supply Bulletin 44, 

51 p.

Manga, Michael, 1996, Hydrology of spring-dominated 

streams in the Oregon Cascades: Water Resources 

Research, v. 32, no. 8, p. 2435–2439.

Manga, Michael, 1997, A model for discharge in springdominated streams and implications for the transmissivity and recharge of Quaternary volcanics in 

the Oregon Cascades: Water Resources Research, 

v. 33, no. 8, p. 1813–1822.

Mantua N.J., Hare, S.J., Zhang, Yuan, Wallace, J.M., and 

Francis, R.C., 1997, A Pacific interdecadal climate 

oscillation with impacts on salmon production: Bulletin 

of the American Meteorology Society, v. 78, no. 6, 

p. 1069–1079.
84

McDonald, M.G., and Harbaugh, A.W., 1988, A modular 

three-dimensional finite-difference ground-water 

flow model: U.S. Geological Survey, Techniques of 

Water-Resources Investigations, book 6, chap. A1, 

586 p.

McFarland, W.D., and Ryals, G.N., 1991, Adequacy of 

available hydrologic data for evaluation of declining 

ground-water levels in the Fort Rock Basin, southcentral Oregon: U.S. Geological Survey WaterResources Investigations Report 89–4057, 47 p.

Miller, D.W., 1986, Ground-water conditions in the Fort 

Rock Basin, northern Lake County, Oregon: State of 

Oregon, Water Resources Department, Ground-Water 

Report No. 31, 196 p.

Newcomb, R.C., 1953, Ground water available for irrigation 

in the Fort Rock Basin, northern Lake County, Oregon: 

U.S. Geological Unpublished Open-File Report, 5 p.

Oregon Water Resources Department, 1965, Compilation 

of surface water records of Oregon: 219 p.

Packard, F.A., Hansen, A.J., 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.

Prudic, D.E., 1989, Documentation of a computer program 

to simulate stream-aquifer relations using a modular, 

finite-difference, ground-water flow model: 

U.S. Geological Survey Open-File Report 88–729, 

113 p.

Redmond, K.T., and Koch, R.W., 1991, Surface climate 

and streamflow variability in the western United States 

and their relationship to large-scale circulation indices: 

Water Resources Research, v. 27, no. 9, p. 2381–2399. 

Russell, I.C., 1905, Geology and water resources of central 

Oregon: U.S. Geological Survey Bulletin 252, 138 p.

Sceva, J.E., 1960, A brief description of the ground water 

resources of the Deschutes River Basin, Oregon: 

Oregon State Engineer (now Oregon Water Resources 

Department), Salem, Oregon, 55 p.

Sceva, J.E., 1968, Liquid waste disposal in the lava terrane 

of central Oregon: U.S. Department of the Interior, 

Federal Water Pollution Control Administration, 

Technical Projects Branch Report No. FR-4, 66 p., 

plus a 96 page appendix.

Smith, G.A., 1986, Stratigraphy, sedimentology, and 

petrology of the Neogene rocks in the Deschutes basin, 

central Oregon: A record of continental-margin 

volcanism and its influence on fluvial sedimentation in 

an arc-adjacent basin: Corvallis, Oregon, Oregon State 

University doctoral dissertation, 464 p. (also Richland, 

Washington, Rockwell Hanford Operations Report 

RHO-SA-555P, 260 p.)

Stearns, H.T., 1931, Geology and water resources of the 

middle Deschutes River Basin, Oregon: U.S. Geological Survey Water-Supply Paper 637–D, 220 p.

Tanaka, H.H., Hansen, A J., and Skrivan, J.A., 1974, 

Digital-model study of ground-water hydrology, 

Columbia Basin Irrigation Project area, Washington: 

State of Washington Department of Ecology WaterSupply Bulletin 40, 60 p.

Taylor, G.H., 1993, Normal annual precipitation, State of 

Oregon: Corvallis, Oregon State University, Oregon 

Climate Service, map.