Evaluation of Approaches for Managing Nitrate 

Loading from On-Site Wastewater Systems near 

La Pine, Oregon

Scientific Investigations Report 2007–5237

Prepared in cooperation with the Oregon Department of Environmental Quality and 

Deschutes County

U.S. Department of the Interior

U.S. Geological Survey
Cover: Photograph of Little Deschutes River near La Pine, Oregon. (Photograph taken 

by David Morgan, U.S. Geological Survey, 2000.)
Evaluation of Approaches for Managing 

Nitrate Loading from On-Site Wastewater 

Systems near La Pine, Oregon

By David S. Morgan, Stephen R. Hinkle, U.S. Geological Survey, and 

Rodney J. Weick, Oregon Department of Environmental Quality

Prepared in cooperation with the Oregon Department of Environmental Quality and 

Deschutes County

Scientific Investigations Report 2007–5237

U.S. Department of the Interior

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

DIRK KEMPTHORNE, Secretary

U.S. Geological Survey

Mark D. Myers, Director

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

For product and ordering information:World Wide Web: http://www.usgs.gov/pubprod

Telephone: 1-888-ASK-USGS

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

natural hazards, and the environment:

World Wide Web: http://www.usgs.govTelephone: 1-888-ASK-USGS

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

by the U.S. Government.

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

reproduce any copyrighted materials contained within this report.

Suggested citation:

Morgan, D.S., Hinkle, S.R., and Weick, R.J., 2007, Evaluation of approaches for managing nitrate loading from on-site 

wastewater systems near La Pine, Oregon: U.S. Geological Survey Scientific Investigations Report 2007-5237, 66 p.
iii

Contents

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

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

Purpose and Scope ..............................................................................................................................4

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

Conceptual Model of Ground-Water System............................................................................................5

Geologic Setting....................................................................................................................................5

Three-Dimensional Hydrogeologic Model .......................................................................................7

Subsurface Geologic Data .........................................................................................................7

Lithologic Unit Sections and Fluvial Sediment Mapping.......................................................7

Transition Probability Geostatistical Modeling.......................................................................8

Hydraulic Conductivity................................................................................................................9

Ground-Water Recharge ...................................................................................................................12

Ground-Water Flow............................................................................................................................14

Ground-Water Discharge ..................................................................................................................15

Streams and Springs.................................................................................................................15

Evapotranspiration ....................................................................................................................19

Wells. .........................................................................................................................................19

Nitrogen Fate and Transport.............................................................................................................20

Nitrogen Loading from On-Site Wastewater Systems .................................................................22

Nitrate Fate and Transport Simulation Models.......................................................................................24

Modeling Approach............................................................................................................................24

Model Parameters..............................................................................................................................26

Hydraulic Conductivity..............................................................................................................26

Porosity........................................................................................................................................26

Dispersion ...................................................................................................................................26

Model Stresses ...................................................................................................................................26

Recharge.....................................................................................................................................27

Rivers. .........................................................................................................................................27

Springs.........................................................................................................................................27

Evapotranspiration ....................................................................................................................27

Transect Model ...................................................................................................................................28

Boundaries and Discretization ................................................................................................28

Calibration and Sensitivity........................................................................................................28

Simulated Travel Time and Comparison with Ground-Water Age.....................................31

Study-Area Model...............................................................................................................................31

Boundaries and Discretization ................................................................................................32

Calibration and Sensitivity........................................................................................................32

Calibration of Ground-Water Flow.................................................................................32

Comparison of Simulated and Measured Nitrate Concentrations ...........................35

Management Scenario Simulations ................................................................................................36
iv

Nitrate Loading Management Model .......................................................................................................43

Formulation of Nitrate Loading Management Model ...................................................................43

Constraints on Ground-Water Nitrate Concentration..........................................................43

Constraints on Discharge of Nitrate from Ground Water to Streams...............................45

Constraints on Reduction of Nitrate Loading........................................................................45

Response-Matrix Technique for Solution of Nitrate Loading Management Model ................46

Application of Model..........................................................................................................................48

Sensitivity of Optimal Solution to Water-Quality Constraints.............................................49

Ground-Water Nitrate Concentration............................................................................49

Ground-Water Nitrate Discharge Loading to Streams...............................................50

Sensitivity of Optimal Solution to Nitrate Loading Constraints ..........................................50

Sensitivity of Optimal Solution to Cost Factors.....................................................................51

Spatial Distribution of Loading for Optimal Solution............................................................52

Comparison of Scenario Simulations and Optimal Solution...............................................53

Limitations and Appropriate Use of Models............................................................................................56

Simulation Models..............................................................................................................................56

Management Model...........................................................................................................................57

Summary and Conclusions.........................................................................................................................57

Acknowledgments.......................................................................................................................................59

References Cited..........................................................................................................................................59

Appendix A. Vertical Hydraulic Head Gradient Data from Measurements Made 

on the Deschutes and Little Deschutes Rivers, in the La Pine, Oregon, 

Study Area, October 23– November 4, 2000 ..............................................................................63

Contents—Continued

Plate

Plate 1. Selected cross sections showing lithologic units of the La Pine region, 

Oregon ..............................................................................[In pocket]
v

Figures

Figure 1. Map showing location of La Pine study area, Oregon, and extent of the nitrate

fate and transport model .................................................................. 3

Figure 2. Map showing generalized geology and hydrogeologic units of the La Pine

region, Oregon .............................................................................. 6

Figure 3. Cross sections showing three dimensional hydrogeologic model of the La

Pine, Oregon, study area .................................................................. 10

Figure 4. Graphs showing distributions of estimated hydraulic conductivity values from

slug tests and well-yield data ............................................................ 12

Figure 5. Map showing distribution of mean annual recharge, water-table elevation

contours (June 2000), and locations of monitoring wells in the La Pine,

Oregon, study area ........................................................................ 13

Figure 6. Hydrographs showing water levels in selected wells in the La Pine, Oregon,

study area .................................................................................... 14

Figure 7. Map showing gaining and losing reaches of the Deschutes and Little

Deschutes Rivers and locations of measurement sites for gain-loss surveys

between October 1995 and October 2000 in the La Pine, Oregon, study area ...... 17

Figure 8. Map showing estimated thickness of oxic ground-water layer in the shallow

aquifer in the La Pine, Oregon, study area................................................ 21

Figure 9. Graph showing annual and cumulative estimated nitrate loading from on-site

wastewater systems in the La Pine, Oregon study area, 1960–2005 .................. 23

Figure 10. Map showing spatial relations between the upper Deschutes Basin regional

ground-water model, the La Pine study area model, and the transect model ...... 25

Figure 11. Plan and section views of the transect model showing simulated

water-levels, ground-water travel time, and particle paths in the La Pine,

Oregon, study area ........................................................................ 29

Figure 12. Map showing contours of simulated and observed heads (June 2000) for the

La Pine, Oregon, study area ............................................................... 34

Figure 13. Graph showing simulated head residuals and observed heads (June 2000)

from the La Pine, Oregon, study area ................................................... 35

Figure 14. Boxplots showing measured and simulated nitrate concentrations in the La

Pine, Oregon, study area, 1999 ............................................................ 36

Figure 15. Map showing simulated nitrate concentrations near the water table in the La

Pine, Oregon, study area, 1999 ............................................................ 37

Figure 16. Graph showing historical nitrate loading from on-site wastewater systems

and eight nitrate loading scenarios tested with the study-area model ............... 39

Figure 17. Map showing simulated equilibrium ground-water nitrate concentrations

near the water table for the base scenario in the La Pine, Oregon, study

area .......................................................................................... 41

Figure 18. Map showing simulated equilibrium ground-water nitrate concentrations

near the water table for 20 milligrams N per liter advanced treatment on-site

wastewater systems in the La Pine, Oregon, study area .............................. 42

Figure 19. Map showing locations of management areas and ground-water nitrate

concentration constraint locations in the Nitrate Loading Management

Model for the La Pine, Oregon, study area ............................................. 44
vi

Figure 20. Map showing locations of management areas near Burgess Road and

management area 31 in the La Pine, Oregon, study area .............................. 47

Figure 21. Graph showing sensitivity of optimal loading solutions to ground-water

nitrate concentration constraints in the La Pine, Oregon, study area ............... 49

Figure 22. Graph showing sensitivity of optimal loading solutions to constraints on the

minimum reduction in ground-water discharge loading to streams in the La

Pine, Oregon, study area .................................................................. 50

Figure 23. Graph showing sensitivity of optimal solution to minimum decentralized

wastewater treatment performance standards for future homes in the La

Pine, Oregon, study area .................................................................. 51

Figure 24. Diagram showing sensitivity of optimal solutions to relative cost difference of

nitrate loading reduction for existing and future homes in the La Pine,

Oregon, study area ........................................................................ 52

Figure 25. Map showing optimal reduction in nitrate loading from existing homes in the

La Pine, Oregon, study area ............................................................... 54

Figure 26. Map showing optimal reduction in nitrate loading from future homes in the La

Pine, Oregon, study area .................................................................. 55

Figure 27. Graphs showing comparison of loading and water quality between optimal

and nonoptimal management scenarios for the La Pine, Oregon, study area ...... 56

Figures—Continued

Tables

Table 1. Markov chain model parameters and transition probabilities in the final

hydrofacies model........................................................................... 9

Table 2. Hydraulic conductivity estimates from slug test data from wells near LaPine,

Oregon ....................................................................................... 11

Table 3. Basin information for water-level monitoring wells in the La Pine, Oregon,

study area .................................................................................... 15

Table 4. Summary of data from stream gain-loss surveys on the Deschutes and Little

Deschutes Rivers in the La Pine, Oregon, study area, 1995–2000 ..................... 18

Table 5. Values of horizontal and vertical hydraulic conductivity for hydrofacies based

on field data and model calibration in the La Pine, Oregon, study area............... 30

Table 6. Summary of eight on-site wastewater management scenarios tested with the

study-area model in the La Pine, Oregon, study area ................................. 39

Table 7. Summary of model simulation results for eight on-site wastewater

management scenarios tested with the study-area model in the La Pine,

Oregon, study area ........................................................................ 40

Table 8. Response coefficients relating the effects of loading in nearby areas to the

nitrate concentration at a constraint location in management area 31 in the

La Pine, Oregon, study area ............................................................... 48
vii

Conversion Factors, Datums, and Abbreviations and 

Acronyms

Conversion Factors

Multiply By To obtain

acre 0.4047 hectare (ha)

cubic foot per second per mile (ft3

/s/mi) 0.04560 cubic meter per second per kilometer

(m3

/s/km)

cubic foot per second (ft3

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

cubic foot per second per square mile

[(ft3

/s)/mi2]

0.01093 cubic meter per second per square kilometer

[(m3/s)/km2]

foot (ft) 0.3048 meter (m)

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

foot squared per day (ft2

/d)* 0.09290 meter squared per day (m2

/d) 

gallon (gal) 3.785 liter (L) 

gallon per day (gal/d) 0.003785 cubic meter per day (m3

/d)

inch (in.) 2.54 centimeter (cm)

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

kilograms per day (kg/d) 2.205 pound per day (lb/d)

kilogram per year (kg/yr) 2.205 pound per year (lb/yr)

mile (mi) 1.609 kilometer (km)

million gallons per day (Mgal/d) 0.04381 cubic meter per second (m3

/s)

pound, avoirdupois (lb) 0.4536 kilogram (kg) 

pound per year (lb/yr) 0.4536 kilogram year (kg/yr)

square mile (mi2

) 2.590 square kilometer (km2) 

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

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

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

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

or micrograms per liter (μg/L).

Datums

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

(NAVD 29).

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

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

Abbreviations and Acronyms

Abbreviations and Acronyms Meaning

kg N/yr kilograms of nitrogen per year

lb N/yr pounds of nitrogen per year

mg N/L milligrams of nitrogen per liter

my million years

NDP La Pine National On-Site Wastewater Demonstration Project

NLMM Nitrate Loading Management Model

ODEQ Oregon Department of Environmental Quality

PLSS Public Land Survey System

RMSE root mean square error

TDC Transferable Development Credit

USEPA U.S. Environmental Protection Agency

USGS U.S. Geological Survey

Conversion Factors, Datums, and Abbreviations and 

Acronyms—Continued
Abstract

This report presents the results of a study by the U.S. 

Geological Survey, done in cooperation with the Oregon 

Department of Environmental Quality and Deschutes County, 

to develop a better understanding of the effects of nitrogen 

from on-site wastewater disposal systems on the quality of 

ground water near La Pine in southern Deschutes County 

and northern Klamath County, Oregon. Simulation models 

were used to test the conceptual understanding of the system 

and were coupled with optimization methods to develop the 

Nitrate Loading Management Model, a decision-support tool 

that can be used to efficiently evaluate alternative approaches 

for managing nitrate loading from on-site wastewater systems. 

The conceptual model of the system is based on geologic, 

hydrologic, and geochemical data collected for this study, as 

well as previous hydrogeologic and water quality studies and 

field testing of on-site wastewater systems in the area by other 

agencies. 

On-site wastewater systems are the only significant 

source of anthropogenic nitrogen to shallow ground water 

in the study area. Between 1960 and 2005 estimated nitrate 

loading from on-site wastewater systems increased from 3,900 

to 91,000 pounds of nitrogen per year. When all remaining 

lots are developed (in 2019 at current building rates), nitrate 

loading is projected to reach nearly 150,000 pounds of 

nitrogen per year. Low recharge rates (2–3 inches per year) 

and ground-water flow velocities generally have limited the 

extent of nitrate occurrence to discrete plumes within 20–30 

feet of the water table; however, hydraulic-gradient and age 

data indicate that, given sufficient time and additional loading, 

nitrate will migrate to depths where many domestic wells 

currently obtain water. In 2000, nitrate concentrations greater 

than 4 milligrams nitrogen per liter (mg N/L) were detected in 

10 percent of domestic wells sampled by Oregon Department 

of Environmental Quality. 

Numerical simulation models were constructed at transect 

(2.4 square miles) and study-area (247 square miles) scales to 

test the conceptual model and evaluate processes controlling 

nitrate concentrations in ground water and potential 

ground-water discharge of nitrate to streams. Simulation 

of water-quality conditions for a projected future build-out 

(base) scenario in which all existing lots are developed using 

conventional on-site wastewater systems indicates that, at 

equilibrium, average nitrate concentrations near the water 

table will exceed 10 mg N/L over areas totaling 9,400 acres. 

Other scenarios were simulated where future nitrate loading 

was reduced using advanced treatment on-site systems and 

a development transfer program. Seven other scenarios were 

simulated with total nitrate loading reductions ranging from 15 

to 94 percent; simulated reductions in the area where average 

nitrate concentrations near the water table exceed 10 mg N/L 

range from 22 to 99 percent at equilibrium. Simulations also 

show that the ground-water system responds slowly to changes 

in nitrate loading due to low recharge rates and ground-water 

flow velocity. Consequently, reductions in nitrate loading 

will not immediately reduce average nitrate concentrations 

and the average concentration in the aquifer will continue to 

increase for 25–50 years depending on the level and timing of 

loading reduction. The capacity of the ground-water system 

to receive on-site wastewater system effluent, which is related 

to the density of homes, presence of upgradient residential 

development, ground-water recharge rate, ground-water flow 

velocity, and thickness of the oxic part of the aquifer, varies 

within the study area. 

Optimization capability was added to the study-area 

simulation model and the combined simulation-optimization 

model was used to evaluate alternative approaches to 

management of nitrate loading from on-site wastewater 

systems to the shallow alluvial aquifer. The Nitrate Loading 

Management Model (NLMM) was formulated to find the 

minimum reductions from projected future loading required 

to maintain or restore ground-water nitrate concentrations 

or ground-water discharge of nitrate to streams below userspecified levels. Sensitivity analysis of the NLMM showed 

that loading primarily is constrained by nitrate concentration 

in the shallow part of the oxic ground-water system, within 

5–10 feet of the water table. 

Evaluation of Approaches for Managing Nitrate 

Loading from On-Site Wastewater Systems near 

La Pine, Oregon

By David S. Morgan1, Stephen R. Hinkle1, and Rodney J. Weick2

1 U.S. Geological Survey

2 Oregon Department of Environmental Quality
2 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

Introduction

Rural residential areas near La Pine in southern 

Deschutes County and northern Klamath County, Oregon 

(fig. 1) have experienced rapid growth in recent years. 

More than 9,000 residential lots, ranging in size from 0.5 

to 10 acres, are in the area where existing and future homes 

will rely on individual on-site wastewater systems and wells 

for wastewater disposal and water supply. Most existing 

wells have screened intervals within 50 ft of land surface and 

extract water from alluvial sands and gravels that constitute 

the primary aquifer in the area. The water table is shallow, 

typically is less than 20 ft below land surface, and in some 

low-lying areas rises seasonally to within 2 ft of land surface. 

Sandy soils derived from pumice contain little organic matter 

and allow rapid infiltration of on-site wastewater effluent. 

The vulnerability of the shallow aquifer has led to 

concern by County and State land-use and environmentalhealth regulators that ground-water quality may be impaired if 

development continues at planned densities using conventional 

on-site wastewater systems. (In this report, conventional 

on-site systems include standard, pressure-distribution, and 

packed-bed [sand] filter systems). Another potential concern 

is the quality of local streams. The Deschutes and Little 

Deschutes Rivers, which flow through developed areas near 

La Pine, have been listed as “water-quality impaired” for 

temperature and turbidity; nutrient loading from ground water 

has been identified as a potential contributor to excessive algal 

growth in some reaches (Anderson, 2000; Jones, 2003) that 

may exacerbate water-quality concerns.

On-site wastewater systems are the principal source of 

nitrogen to the shallow ground-water system in the La Pine 

area (Century West Engineering, 1982; Oregon Department 

of Environmental Quality, 1994; Hinkle and others, 2007a). 

On-site wastewater systems do not remove nitrogen from 

wastewater; however, most nitrogen is converted from organic 

nitrogen and ammonium to nitrate before it reaches the 

saturated part of the ground-water system (water table). Once 

in the saturated zone, nitrate generally is stable in the presence 

of oxic ground water (ground water that contains dissolved 

oxygen). Adsorption or chemical reactions within the flow 

system do not readily remove nitrate. Nitrate is a human health 

concern because it can cause methemoglobinemia (BlueBaby Syndrome) in infants (http://www.atsdr.cdc.gov/HEC/

CSEM/nitrate/). Nitrogen also is an environmental concern 

as a potential source of nutrient enrichment to streams. 

Nutrient enrichment contributes to algal blooms detected 

in the Deschutes and Little Deschutes Rivers. The U.S. 

Environmental Protection Agency (USEPA) has established 

10 mg N/L as the maximum allowable nitrate concentration 

in drinking water for public water supply systems. Oregon, by 

statute, has established a nitrate concentration of 7 mg N/L as 

the value at which action may be taken to control water-quality 

degradation by regulatory means. 

The present location of the city of La Pine (fig. 1) was the 

first concentrated development within the study area. The first 

building permits recorded in what was then called the “core 

area” date from 1910. In 2006, the core area was incorporated 

as the city of La Pine. Degradation of ground-water quality 

by on-site wastewater systems was first documented in the 

core area in 1979 by Oregon Department of Environmental 

Quality (ODEQ) after samples from 46 wells revealed that 

water from 8 wells contained nitrate concentrations greater 

than 10 mg N/L (Cole, 2006) and a maximum of 26 mg N/L. 

Nitrate concentrations in ground water as much as 41 mg N/L 

were detected in a follow-up study in 1982 (Century 

West Engineering, 1982). These results led to an ODEQ 

administrative rule requiring community sewage collection, 

treatment, and disposal for the core area. In 1993, ODEQ 

sampled 36 wells in residential areas near La Pine as part 

of its statewide Ambient Groundwater Monitoring Program 

and detected nitrate concentrations greater than 2 mg N/L in 

19 wells. They concluded that elevated concentrations were 

caused by anthropogenic influences (on-site wastewater 

systems) (Cole, 2006). Concentrations were greater than 

10 mg N/L in 4 (11 percent) of the 36 wells sampled. In 1994, 

Deschutes County requested that ODEQ further evaluate 

ground-water quality in the area. Water from more than 120 

domestic and public water-supply wells was sampled in 1994 

and 1995 and ODEQ delineated several areas of elevated 

nitrate concentrations underlying the most densely developed 

parts of the region (R.J. Weick, ODEQ, written commun., 

1998; Cole, 2006). As part of the 1995 assessment, ODEQ 

constructed simplified nitrate transport models which included 

basic assumptions on nitrogen chemistry. The ODEQ transport 

models predicted that nitrate concentrations would exceed 

drinking water standards within 20 years of full build-out (R.J. 

Weick, ODEQ, written commun., 1998; Cole, 2006). 

In 1999, a nitrate concentration threshold of 5 mg N/L 

was selected for the La Pine area by the Deschutes County 

Working Group on Groundwater Issues for the South 

Deschutes Basin to serve as a “proactive target to protect 

the Basin’s groundwater quality.” The working group also 

recommended that ODEQ “...address these problems and 

concerns by considering the adoption of a geographic rule...

to protect the Basin’s groundwater resource.” (http://www.

co.deschutes.or.us/download.cfm?DownloadFile=507DDA68-

BDBD-57C1-902813CCF281DDF0). 
Introduction 3

Figure 1. Location of La Pine study area, Oregon, and extent of the nitrate fate and transport model.

OR19-0049_fig 01

Base 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

Publication projection is 

Lambert Conformal Conic

Standard parallels 42o20' and 44o40', 

central meridian -120o30'

Datum is NAD83

T 19 S

T 20 S

T 24 S

T 23 S

T 22 S

T 21 S

R 9 E R 10 E R 11 E

DESCHUTES CO.

L A K E C O .

KLAMATH CO.

EXPLANATION

Sewered area

Alluvium

Basalt, andesite

Existing home within model 

area and not served by 

sewer

Transect well locations 

Geologic unit

Paulina

River

RIVER

DESCHUTES

Deschutes

Little

Creek

Creek Crescent

River

Fall

Spring

River

Long Creek

Wickiup

Reservoir

97

97

31

31

Masten Road

River

Spring

Century

South

Road

Road

Butte

Finley

Drive

Burgess Road

5 KILOMETERS

0 5 MILES

0

43°35'

40'

121°40' 35' 30' 25' 121°20'

45'

50'

43°55'

Extent of study area model

Deschutes

Basin

Cascade

Range

Newberry Upper Volcano Deschutes

Basin

La Pine

region

OREGON

Bend

SunriverSunriver

La 

Pine

Long Prairie

La Pine 

core 

area

Wickiup 

Junction Burgess

transect

Century 

transect
4 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

In 1999, Deschutes County and ODEQ identified the 

need for an improved understanding of the processes that 

affect the transport and fate of nitrogen in the La Pine area 

before making decisions among alternatives for managing 

ground-water quality. To help achieve that understanding, 

Deschutes County and ODEQ applied for and received 

funding from the USEPA to evaluate methods to protect 

ground-water quality in the area (Oregon Department 

of Environmental Quality, 2004a) as part of the La Pine 

National On-Site Wastewater Demonstration Project (NDP). 

The objectives of the NDP also included (1) assessing the 

effectiveness of advanced treatment (denitrifying) systems for 

on-site wastewater and (2) developing a more complete and 

useful understanding of processes that affect nitrogen in the 

ground-water system. The second objective was the subject 

of a cooperative study by the ODEQ, Deschutes County 

Community Development Department (CDD) and the Oregon 

Water Science Center of the U.S. Geological Survey (USGS). 

The primary objective of this study was to develop 

a thorough understanding of the hydrologic and chemical 

processes that affect the movement and fate of nitrogen within 

the shallow aquifers of the La Pine region. A secondary 

objective was to provide a method for analyzing the effects 

of existing and future development on water quality. 

This understanding will provide local and State resource 

management agencies with information and tools needed 

to determine the probable effects of present and future land 

use on nitrogen concentrations in the shallow aquifer and on 

nitrogen loading from ground water to the Deschutes and 

Little Deschutes Rivers. 

Purpose and Scope

This report describes conceptual, simulation, and 

management models of the ground-water system near the 

community of La Pine in central Oregon. The description 

of the conceptual model provides a background for the 

development of the simulation and management models and 

includes the geologic framework, hydrologic processes, and 

processes affecting nitrate transport and fate. The description 

of the simulation model includes a discussion of how the 

processes and boundaries of the system were represented in 

a computer simulation model, how the model was calibrated, 

and the results of predictive simulations. The simulation model 

was enhanced as a decision-support tool (management model) 

by incorporating optimization techniques that allow users to 

identify management solutions that will meet water quality 

goals. The discussion of the management model includes 

a description of the method of incorporating optimization 

techniques with the simulation model, the formulation of 

the management problem, and analysis of the sensitivity 

of optimal management solutions to the values of various 

constraints.

Description of Study Area

The La Pine study area encompasses about 250 mi2

within the Deschutes River drainage basin in central Oregon 

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

tributaries including the Little Deschutes, Spring, and Fall 

Rivers. Land-surface elevation ranges from about 4,000 ft near 

the Deschutes River at the northern boundary of the study 

area, to nearly 5,700 ft at the peaks of volcanic buttes in the 

northwestern corner of the study area. Most of the study area 

lies in the relatively low-relief alluvial plain of the Deschutes 

and Little Deschutes Rivers at elevations between 4,150 and 

4,300 ft. The populated areas in the study area include low and 

medium density rural-residential subdivisions adjacent to and 

between the Deschutes and Little Deschutes Rivers extending 

northward from La Pine to the community of Sunriver (fig. 1). 

The subdivisions are surrounded by Federal lands managed by 

the U.S. Forest Service and Bureau of Land Management. 

The study area boundaries were selected to include the 

most densely populated parts of southern Deschutes County 

and northern Klamath County, where on-site wastewater 

systems are the predominant method of wastewater disposal 

for existing and future residential development. The 

boundaries also were selected to include the principal area of 

the shallow alluvial aquifer that provides drinking water for 

most of the population in the La Pine area. The data collection 

and analysis for this study was focused in the 247 mi2

 area of 

the study-area flow and transport simulation model (fig. 1). 

Residential development in the area began to accelerate in 

the 1960s (Century West Engineering, 1982) and almost 3,000 

new residential lots were created in the 1970s in response 

to demand for vacation homes and full-time residences. 

Although lots range in size from 0.5 to more than 10 acres, 

58 percent of lots in the study area in 2000 were less than 

1 acre and 82 percent were less than 2 acres. The population 

of the area was less than 1,000 in 1960, and increased to about 

5,600 in 1981 (Century West Engineering, 1982), and was 

approximately 14,000 in 2000 (U.S. Bureau of the Census, 

2000). 

Land use is primarily low- to medium-density residential. 

Commercial and medium density residential areas are in the 

incorporated City of La Pine (including the La Pine “core 

area” and Wickiup Junction area) and industrial development 

is immediately east of the core area. Agricultural lands cover 

less than 4 percent of the study area and most of that area is 

nonirrigated pasture (Tim Berg, Deschutes County Community 

Development Department, written commun., 2004).

The climate of the area primarily is controlled by 

eastward moving air masses from the Pacific Ocean. 

Orographic precipitation in the Cascade Range results in more 

than 200 in/yr in some locations, although precipitation rates 

in the lower part of the Deschutes River basin average as little 

as 10 in/yr. Mean annual (1971–2000) precipitation in the La 

Pine area ranges from 16 to 24 in. (Daly and Gibson, 2002). 
Conceptual Model of Ground-Water System 5

Conceptual Model of 

Ground-Water System

Geologic Setting

The La Pine study area of the upper Deschutes River 

basin is a sediment-filled graben feature that lies within the 

geologically complex transition area between three geologic 

provinces: the Cascade Range, the Basin and Range Province, 

and the High Lava Plains. The La Pine study area shares 

attributes characteristic of each of these geologic provinces. 

The Cascade Range forms the western margin of the La Pine 

study area. Newberry volcano forms the eastern study area 

margin and transition between the extensional Basin and 

Range Province and the High Lava Plains. North to northnorthwest trending and northeast trending en echelon faults 

define the eastern boundary of the graben forming the La Pine 

study area. 

The Cascade Range is a constructional feature of northsouth trending eruptive centers that extends from northern 

California to southern British Columbia and has been 

volcanically active for the past 35 million years (Sherrod and 

Pickthorn, 1989). In central Oregon, these volcanic eruptive 

centers are stratovolcanoes, such as the North, Middle, and 

South Sister, and Mount Jefferson, all with elevations greater 

than 10,000 ft (Lite and Gannett, 2002). 

The Basin and Range Province is a region of crustal 

extension that covers most of the western United States and 

is characterized by north to northwest trending sub parallel 

fault-bounded down-dropped grabens forming fault-block 

ranges with basins typically 10–20 mi wide. In central 

Oregon, the La Pine graben is defined by north-northwest 

trending and northeast trending down-dropped faults (Allen, 

1966; MacLeod and Sherrod, 1992) forming the sedimentfilled 10–15-mi wide subbasin with characteristics similar 

to those of Basin and Range extensional features along the 

northwestern transitional boundary in northeastern California 

and northwestern Nevada. 

The north–south trending graben that underlies the La 

Pine study area is estimated to have down-dropped 1,800 to 

2,400 ft based on gravity and aeromagnetic anomalies (Couch 

and Foote, 1985; MacLeod and Sherrod, 1992) Formation 

of the basin began during the mid-Pleistocene between 

0.6 and 1 my (million years) ago (Couch and Foote, 1985; 

Gettings and Griscom, 1988; Sherrod and Pickthorn, 1989; 

MacLeod and Sherrod, 1992; MacLeod and others, 1995) 

and subsequently has been filled with several hundred feet of 

sediment. Nearly 1,400 ft of sediment have been penetrated by 

deep water-supply wells drilled in the study area. Descriptions 

of sediment lithology in drillers’ reports for water wells 

indicate that graben formation was concurrent with volcanism 

and late Pleistocene glacial outwash deposition. Periods of 

active volcanism and quiescence and fluvial and lacustrine 

deposition during graben development created a complex 

sequence of intercalated lava flows, ignimbrites, and alluvial 

deposits. The depositional history is further complicated by 

fluvial-lacustrine deposition during the Pleistocene.

Many basin-fill sediments in the La Pine study area are 

fine-grained lacustrine silt and clay. Within these fine-grained 

deposits are fine to coarse fluvial sand and gravel channelfill deposits and discontinuous cinder, pumice- and ash-fall 

beds. Methane, ammonium, and reduced iron detected in 

water well samples (Hinkle and others, 2007a) from these 

fine-grained deposits indicate a predominant quiescent marsh 

and lake depositional environment with episodic volcanic 

deposition. This low-energy depositional environment may 

be related to the onset and development of Newberry volcano 

about 0.7 my ago. MacLeod and others (1995) indicated that 

Newberry lavas backed up against Cascade Range lava flows, 

blocking the channel of the Deschutes River which created a 

lake and marsh environment over much of the study area. A 

pumice bed exposed near the top of the lacustrine deposits 

and about 35 ft (10.8 m) below the fluvial/lacustrine contact 

at Pringle Falls on the Deschutes River is 0.22 to 0.17 my old 

(Herrero-Bervera and others, 1994).

About 0.2 my ago the depositional environment probably 

underwent an abrupt change from lacustrine to predominantly 

fluvial with deposition of heterogeneous silt, fine to coarse 

sand, gravel, and pumaceous sand and gravel. These 

deposits are likely associated with Pleistocene glaciation 

of the Cascade Range and Newberry Volcano. These highenergy deposits are capped by a study area wide 3–5 ft-thick 

pumice- and ash-fall deposit from the Mt. Mazama eruption 

of 7,627 ± 150 years ago (Zdanowicz and others, 1999). The 

Deschutes, Little Deschutes, and Fall Rivers have reworked 

and down-cut through the Mt. Mazama pumice and ash-fall 

deposit. The high degree of heterogeneity noted in lithologic 

descriptions in drillers’ reports indicate an active fluvial 

depositional environment in the central and southern parts of 

the La Pine study area and deposition more characteristic of a 

lacustrine environment in the northern part of the study area. 

Based on interpretation of more than 460 water-well logs, 

fluvial silt, sand, and gravel deposits were determined to range 

in thickness from less than 10 ft in the northern study area, 

to as much as 100 ft in the central and southern parts of the 

study area (fig. 2). A thin veneer of gravel overlying a paleosol 

was observed in several monitoring wells in the study area. 

The presence of the paleosol may represent post-Pleistocene 

glacial soil development. The thin sand and gravel veneer 

overlying the paleosol may represent the brief early Holocene 

glacial event that predates the mid-Holocene Mt. Mazama 

eruption. 
6 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

Figure 2. Generalized geology and hydrogeologic units of the La Pine region, Oregon.

OR19-0049_fig 02

Base 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

Publication projection is 

Lambert Conformal Conic

Standard parallels 42o20' and 44o40', 

central meridian -120o30'

Datum is NAD83

R 9 E R 10 E R 11 E

T 21 S

T 22 S

T 23 S

T 24 S

T 20 S

T 19 S

DESCHUTES CO.

L A K E C O .

KLAMATH CO.

Spring

River

Fall

River

Crescent Creek

Creek

Little

Deschutes

DESCHUTES

RIVER

River

Paulina

Long Creek

5 KILOMETERS

0 5 MILES

0

Wickiup

Reservoir

97

97

31

31

Burgess Road

Drive

Finley

Butte

Road

Road

South

Century

Spring

River

Masten Road

SunriverSunriver

La 

Pine

43°55'

50'

45'

121°40' 35' 30' 25' 121°20'

40'

43°35'

A’

D

C

B

D’

C’

B’

A

Extent of study area model

EXPLANATION

Geologic unit

Alluvium

 Basalt, andesite

Thickness of fluvial

 deposits, in feet within 

 model boundary

0-20 

20-40 

40-60 

60-80

80-100

 Trace of section shown in 

 plate 1

 Other sections

Well used to construct 

 sections

A A’

Deschutes

Basin

Cascade

Range

Newberry Upper Volcano Deschutes

Basin

La Pine

region

OREGON

Bend
Conceptual Model of Ground-Water System 7

Three-Dimensional Hydrogeologic Model

A detailed characterization of subsurface heterogeneity 

can substantially improve the reliability of models of 

ground-water contaminant transport (Fogg, 1986; Anderson 

1987; Johnson and Dreiss, 1989). In this study, standard 

hydrogeologic interpretation and analysis techniques, 

including construction of two-dimensional geologic sections 

and surface maps, were used in conjunction with transition 

probability geostatistics to develop a three-dimensional 

hydrogeologic model that represents the heterogeneity of the 

complex glaciofluvial system. 

Subsurface Geologic Data

The primary source of subsurface geologic data was 

the nearly 5,400 drillers’ reports available (as of 1999) for 

domestic and public water supply wells in the study area. 

Although these reports are plentiful, the quality of lithologic 

descriptions in the reports is inconsistent. The spatial 

distribution of domestic wells also is not ideal for developing 

a detailed, three-dimensional geologic model of the groundwater system; as would be expected, domestic wells are 

concentrated in areas of residential development and tend to be 

completed at the shallowest depth that a satisfactory yield can 

be obtained. A subset of 346 wells was visited as part of this 

and previous studies (Gannett and others, 2001) to determine 

accurate locations and collect water-level and other data. The 

subset of wells visited was selected to provide the best spatial 

distribution within the study area and good descriptions of 

geologic materials penetrated by the well. Drillers’ reports for 

an additional 118 wells were selected to provide information 

in areas where field-located wells were not available. Although 

these wells were not visited, the drillers’ reports included 

superior descriptions of geologic materials and well locations 

could be accurately estimated using tax lot information. 

Descriptions of geologic materials were transcribed to 

a database from the drillers’ reports using a standardized set 

of lithologic descriptors developed for this study. A twoletter descriptor was assigned to the primary lithology, and, 

if needed, secondary lithology reported by the driller for 

each depth interval. For example, if the driller reported a 

layer of “gravel with sand” between 40 and 80 ft, the interval 

was assigned a primary descriptor of GR for gravel and a 

secondary descriptor of SA for sand. This system retained 

most of the detail of the original description by the driller and 

allowed for comparison and analysis of lithology between 

wells. Thirty-six descriptors were used to describe primary 

and secondary lithology. Combinations of the primary and 

secondary descriptors resulted in 157 unique lithologic 

descriptors for the 464 wells used in the analysis. 

Lithologic Unit Sections and Fluvial Sediment 

Mapping

The 157 unique lithologic descriptors were grouped 

into six lithologic units. Five lithologic units comprised 

unconsolidated sediments: clay-silt, pumice-sand, sand, 

sand-gravel, and gravel. The sixth lithologic unit represented 

consolidated and semi-consolidated volcanic rocks such 

as basalt, basaltic andesite, and tuff. The lithologic unit 

data for each well were stored in a Geographic Information 

System (GIS) database that facilitated data interpretation and 

construction of 34 two-dimensional cross sections through the 

study area showing the thickness and extent of the primary 

lithologic units. Locations of the sections are shown in figure 2

and selected sections are shown on plate 1. 

Although deep wells are somewhat sparse in the study 

area, the sections show that as much as 100 ft of fluvial silt, 

sand, and gravel overlie predominately fine-grained lacustrine 

sediments throughout the study area (plate 1). A relatively 

sharp transition exists between the overlying fluvial sediments 

and the basal lacustrine silt and clay, basalt, and volcaniclastic 

deposits. The elevation of the top of the lacustrine sediments 

was mapped using data from the 464 wells (fig. 2). Infrequent 

sand and gravel lenses as well as pumice and ash beds are 

present within the basal lacustrine sediments, but because of 

their discontinuous nature, they are not regionally significant 

sources of water to wells. The thickness of the fluvial 

sediments that constitute the primary aquifers in the study area 

was computed by subtracting the elevation of the top of the 

basal lacustrine sediments from land surface elevation (fig. 2). 

The greatest thicknesses of fluvial sediments are in the central, 

east-central, and south-central parts of the study area where 

higher elevations occur within the study area. Deposition 

of Pleistocene alluvial sediments from fans emanating from 

Newberry volcano contributes to the thickness of coarse 

sediments on the eastern margin of the study area. The total 

thickness (unsaturated and saturated) of fluvial sediments 

exceeds 60 ft over much of this area. The fluvial sediments 

are thinner in three areas: (1) within the floodplain of the 

Little Deschutes River where they have been eroded, (2) at the 

margins of the study area where they overlie shallow tuffs and 

basaltic and andesitic volcanic rocks, and (3) in the northern 

part of the study area where lacustrine sediments are closer to 

the surface (fig. 2). 

Interpretation of the geologic data during construction 

of the two-dimensional sections revealed that there is a 

greater degree of heterogeneity within the fluvial part of the 

system than could be delineated using traditional methods 

of interpretation using drillers’ reports as the primary data 

source. The drillers’ reports were useful for defining the 

vertical heterogeneity at each well location; however, because 
8 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

of the distance between wells and the uncertainty in the 

lithologic descriptions, it was not possible to accurately 

represent the lateral heterogeneity of the fluvial hydrofacies. 

Another consideration was the difficulty in constructing a 

fully three-dimensional hydrogeologic model from the twodimensional representations in the lithofacies sections. 

Important information on the nature and proportion 

of lithologic units and dimensions and interconnectedness 

of depositional features was gained through the process of 

constructing the lithologic unit sections. This information 

and the location of the lower boundary of the fluvial aquifer 

system were used as supplemental information in the 

development of the three-dimensional hydrogeologic model 

described in the next section.

Transition Probability Geostatistical Modeling

A transition probability geostatistical approach was 

used to model the heterogeneity of the fluvial sediments that 

constitute the shallow aquifer near La Pine. The method was 

applied using the suite of programs, Transition Probability 

Geostatistical Software (T-PROGS), documented by Carle 

(1999). T-PROGS implements a transition probability/Markov 

approach to geostatistical analysis and simulation of spatial 

distributions of categorical variables (such as hydrofacies). 

The approach uses the probabilities of lateral and vertical 

transitions between hydrofacies to define spatial variability. 

The transition probabilities, volumetric proportion of facies, 

and mean lengths of depositional features are the parameters 

used to fit a three-dimensional Markov chain model to 

geologic data from drillers’ logs and subjective knowledge 

based on an understanding of the depositional environment. 

Once a satisfactory Markov chain model is obtained, it is used 

in sequential indicator simulation with simulated annealing 

to produce realizations of the subsurface facies distributions 

that are conditioned using the observed geologic data. One 

advantage of the transition probability method is the ability to 

easily incorporate subjective or “soft” geologic information 

into the Markov chain model. Details on the theoretical 

development of the transition probability method, as wells as 

examples and comparisons with other methods, are given in 

Carle and Fogg (1996, 1997), Carle and others (1998), and 

Weissmann and others (1999).

The transition probability method was used to model the 

distribution of facies within the fluvial sediments only. The 

lacustrine sediments were not included because the transition 

probability method is predicated on the assumption that all 

the sediments within the modeled area were emplaced under 

similar depositional conditions. 

 A hydrofacies is defined as one or more lithologic units 

that have similar hydraulic characteristics. Three hydrofacies 

were used in the hydrogeologic model: (1) clay-silt, (2) sand, 

and (3) gravel. The five lithologic units used to subdivide the 

fluvial sediments in the sections on plate 1 were reduced to 

three hydrofacies by combining units with similar hydraulic 

characteristics. The clay-silt hydrofacies represent flood 

plain deposits with low hydraulic conductivity. The sand 

hydrofacies include varying amounts of silt and represent 

levee and proximal overbank deposits with moderate hydraulic 

conductivity. The gravel hydrofacies includes varying amounts 

of sand and represents channel deposits with high hydraulic 

conductivity. 

In the rare case where there are abundant, high 

quality core data, geophysical logs, and detailed lithologic 

descriptions, the parameters for the Markov model can be 

estimated directly from bivariate statistics (variograms) for 

the hydrofacies. More commonly, as in the La Pine study 

area, most available geologic data are obtained from drillers’ 

reports. Another limitation of the well data stems from the lack 

of correspondence between the scale of depositional features 

(for example, channel deposits) that have mean lengths of 

10s or 100s of feet, and the spacing of data points (wells) 

that may be 100s or 1,000s of feet. The disparity in the scales 

of the system and the data available to characterize it makes 

quantifying lateral transition probabilities based solely on well 

data difficult. However, the transition probability approach 

provided a means of including the geologic knowledge gained 

through the process of developing the two-dimensional 

lithologic unit sections. 

The mean lengths of the high-hydraulic conductivity 

gravel hydrofacies and the low hydraulic conductivity claysilt hydrofacies were estimated by manually measuring the 

length of lithologic units interpreted in the 34 lithologic unit 

sections (examples shown on plate 1) and computing average 

lengths for each hydrofacies in the lateral (strike and dip) and 

vertical dimensions. The depositional strike was aligned with 

the northeast-southwest strike of the study area (fig. 2). The 

estimates of lateral mean lengths for the gravel hydrofacies 

based on the lithologic unit sections were larger than expected 

for channel deposits typically found in fluvial systems 

analogous to the Little Deschutes and Deschutes River 

systems in the La Pine study area. The lateral mean lengths for 

the gravel hydrofacies were reduced based on measurements 

of present channel geometry, and meander and oxbow features 

identified from aerial photography. The estimated mean 

lengths for the gravel hydrofacies were modified through 

a trial and error process until the Markov chain model 

produced a hydrofacies realization that was consistent with 

geomorphologic observations. 
Conceptual Model of Ground-Water System 9

The volumetric proportions of the hydrofacies were 10, 

50, and 40-percent for gravel, sand, and clay-silt, respectively. 

The sand hydrofacies was selected as the “background” 

category because it represented the highest proportion of 

materials in the system. The transition probabilities (table 1) 

need only be specified for the nonbackground hydrofacies. 

The probabilities are expressed as decimal percentages; for 

example, there is a 30-percent probability of transitioning 

from the gravel facies to the clay-silt facies in the strike or dip 

direction, but only a 15-percent probability of transitioning 

from clay-silt to gravel. 

The grid used to create the fluvial hydrofacies model had 

cell dimensions of 500 ft horizontally along the depositional 

strike and dip and 5 ft vertically. These cell dimensions 

were sufficiently small to allow good representation of the 

gravel and clay-silt facies (within the sand matrix) with mean 

horizontal and vertical lengths of 1,500–8,100 ft and 8–10 ft, 

respectively (table 1). The overall dimensions of the modeled 

region are 138,000 by 50,000 ft along the depositional 

strike and dip, and 120 ft along the vertical dimension. 

The T-PROGS program was used with the Markov chain 

parameters to generate the three-dimensional distribution of 

the fluvial hydrofacies. The elevation of the base of the fluvial 

sediments was contoured using drillers’ log data; and cells 

within the hydrogeologic model below the base of the fluvial 

sediments were assigned to the lacustrine hydrofacies. The 

surficial geologic map (fig. 2) was used to assign cells to the 

basalt hydrofacies and land surface elevation from a DEM 

was used to “remove” cells from the hydrogeologic model that 

were above land surface. Sections through the complete threedimensional hydrogeologic model are shown in figure 3.

Once the appropriate parameters for the Markov Chain 

model are determined, many equally probable realizations 

of the hydrostratigraphy of the fluvial hydrofacies can be 

simulated using the T-PROGS program. Each realization 

would be consistent with the data from drillers’ logs 

and interpretations made during construction of the twodimensional sections used to estimate the mean lengths. 

The capability to generate many realizations of the 

hydrostratigraphy represents a potentially powerful tool 

for evaluating the effects of hydrogeologic uncertainty on 

simulation models. Using multiple realizations, uncertainty 

in simulated nitrate concentrations resulting from 

uncertainty in the hydrogeologic model could be evaluated. 

While time consuming, this capability could provide 

important information on the sensitivity of model results to 

hydrogeologic uncertainty. 

This capability would be especially useful for 

contaminant transport models in highly heterogeneous 

systems with large contrasts in hydraulic conductivity between 

hydrofacies where the objective was to simulate individual 

contaminant plumes. This approach was not used in the 

current study because (1) although the fluvial sediments in the 

La Pine study area are somewhat heterogeneous, the hydraulic 

conductivities of the hydrofacies that compose them fall within 

a relatively small range; and (2) the purpose of the model in 

this study was to estimate mean concentrations over relatively 

large areas as a tool for evaluating watershed-scale wastewater 

management approaches. Variability at the scale of individual 

model cells caused by uncertainty in the hydrogeologic model 

would not impose significant limitations on the use of the 

model. 

Hydraulic Conductivity

Hydraulic conductivity is the characteristic of a porous 

medium that describes its ability to transmit water, and is 

expressed in units of length per unit time (for example, feet per 

day). The hydraulic conductivity of most geologic materials 

is vertically anisotropic; that is, hydraulic conductivity is 

greater in the horizontal direction than in the vertical direction. 

Simulation of three dimensional ground-water flow, as was 

done for this study, requires estimates of both the horizontal 

and vertical components of hydraulic conductivity. Horizontal 

hydraulic conductivity estimates were made for alluvial 

sediments and volcanic rocks using data from well-yield tests 

reported in 221 drillers’ reports and slug tests done in 24 

monitoring wells constructed for this study. 

Transmissivity is another measure of the ability of 

porous materials to transmit water, and is equal to the product 

of hydraulic conductivity and the saturated thickness of the 

material. Transmissivity has units of length squared per time 

(for example, feet squared per day) and can be estimated using 

the Theis nonequilibrium equation and solving as reported in 

Vorhis (1979). Minimum data required from the well-yield 

test include pumping rate, drawdown, and time pumped. Only 

field-located wells with complete construction information 

that included well depth, casing diameter, and open or 

screened intervals were used in the analysis. The method of 

Vorhis (1979) yields a solution to the Theis equation in which 

transmissivity is expressed as a function of storage coefficient. 

Table 1. Markov chain model parameters and transition 

probabilities in the final hydrofacies model.

Hydrofacies Mean length (feet)

Strike Dip Vertical

Gravel 3,000 1,500 10

Clay-silt 8,100 5,600 8

Hydrofacies 

(From To)

Transition probability 

(decimal percentage)

Gravel Clay-silt 0.3 0.3 0.4

Clay-silt Gravel .15 .15 .25
10 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

OR19-0049_fig 03

Above ground

Lacustrine clay-silt

Volcanics: basalt, andesite

Fluvial gravel

Fluvial sand

Fluvial clay-silt

EXPLANATION

Hydrofacies

4,350

4,300

4,250

4,200

4,150

4,100

ELEVATION, IN FEET

4,350

4,300

4,250

4,200

4,150

4,100

ELEVATION, IN FEET

4,350

4,300

4,250

4,200

4,150

4,100

ELEVATION, IN FEET

4,350

4,300

4,250

4,200

4,150

4,100

ELEVATION, IN FEET

4,350

4,300

4,250

4,200

4,150

4,100

ELEVATION, IN FEET

0

20

40

60

80

100

120

140

160

180

200

220

240

260

275

MODEL ROW NUMBER

4,350

4,300

4,250

4,200

4,150

4,100

ELEVATION, IN FEET

MODEL COLUMN NUMBER

0 10 20 30 40 50 60 70 80 90 100

Notes: Each row and column is 

500 feet wide. Model is truncated 

at elevation of 4,350 feet. 

Figure 3. Three dimensional hydrogeologic model of the La Pine, Oregon, study area.
Conceptual Model of Ground-Water System 11

Appropriate transmissivity values were selected by assuming a 

storage coefficient in the range normally assumed for confined 

aquifers (0.001) for deep wells and a storage coefficient 

(specific yield) in the range normally assumed for unconfined 

aquifers (0.10) for shallow wells. Wells were assigned 

unconfined storage coefficients if their open intervals did not 

extend more than 30 ft below land surface in the northern part 

of the study area or 50 ft in the southern part (south of the 

confluence of Paulina Creek and the Little Deschutes River) 

(fig. 1). Horizontal hydraulic conductivity was estimated at 

each well by dividing the estimated transmissivity by the 

length of the perforated and (or) uncased intervals of the well. 

Horizontal hydraulic conductivities computed using 

the well-yield test data range from 0.15 to 2,500 ft/d, with a 

geometric mean of 60 ft/d and median of 21 ft/d. Gannet and 

others (2001) reported similar results from an analysis of 175 

well-yield tests in the La Pine basin, in which the median 

transmissivity was 901 ft2

/d. If a typical effective aquifer 

thickness of 50 ft is assumed, the transmissivity reported 

by Gannett and others (2001) is equivalent to a hydraulic 

conductivity of 18 ft/d. If well-yield data are reported 

accurately, the principal sources of error using this approach 

are head losses due to poor well construction and uncertainty 

in storage coefficient estimates. Sensitivity analysis showed 

that transmissivity, and therefore hydraulic conductivity, 

computed for a range of storage coefficients spanning three 

orders of magnitude (10-1 to 10-4) differed by an average of 

50 percent. Uncertainty related to storage coefficient is not 

considered an important source of error. Overestimation of 

hydraulic conductivity can occur using this method if a well 

only partially penetrates a thick, homogeneous aquifer in 

which substantial vertical flow toward the well can occur. The 

effects of partial penetration are not likely to be important 

in data from wells in the La Pine study area because the 

heterogeneity of the alluvial sediments would limit vertical 

flow near the well.

Slug tests were done at 24 monitoring wells installed 

along two transects for this study. One transect was located in 

the northern part of the study area near South Century Drive 

and the other was located in the south-central part of the 

area near Burgess Road (fig. 1). The Burgess Road transect 

included tests of 17 wells at 6 sites and the South Century 

Drive transect included tests of 8 wells at 6 sites (table 2). 

Table 2. Hydraulic conductivity estimates from slug test data from wells near LaPine, Oregon.

[Slug tests were conducted in April and September 2000. Latitude and Longitude: in degrees, minutes, seconds; North American Datum of 1927. Altitude: 

referenced to National Geodetic Vertical Datum of 1929. Abbreviations: USGS: U.S. Geological Survey. Symbols: –, no data] 

USGS site

identification

No.

Project well 

name Local No.

State

identification

No.

Latitude Longitude Altitude

(feet)

Depth

(feet)

Hydraulic 

conductivity

(feet per day)

434202121301501 Burgess 5.1 22.00S/10.00E-02BBC01 DESC 52678 43°42'01.89" 121°30'15.21" 4,210.9 23.9 15

434202121301502 Burgess 5.2 22.00S/10.00E-02BBC02 DESC 52677 43°42'01.89" 121°30'15.21" 4,211.1 36.9 24

434202121301503 Burgess 5.3 22.00S/10.00E-02BBC03 DESC 52676 43°42'01.89" 121°30'15.21" 4,210.9 49.0 .6

434203121305401 Burgess 4.1 22.00S/10.00E-03ABC01 DESC 52681 43°42'02.68" 121°30'53.68" 4,223.6 26.5 18

434203121305402 Burgess 4.2 22.00S/10.00E-03ABC02 DESC 52679 43°42'02.68" 121°30'53.68" 4,223.3 33.2 58

434203121305403 Burgess 4.3 22.00S/10.00E-03ABC03 DESC 52680 43°42'02.68" 121°30'53.68" 4,223.3 49.3 6.4

434206121311701 Burgess 3.1 22.00S/10.00E-03BBD01 DESC 52684 43°42'05.78" 121°31'17.04" 4,233.4 22.5 11

434206121311702 Burgess 3.2 22.00S/10.00E-03BBD02 DESC 52682 43°42'05.78" 121°31'17.04" 4,233.3 30.5 45

434206121311703 Burgess 3.3 22.00S/10.00E-03BBD03 DESC 52683 43°42'05.78" 121°31'17.04" 4,233.4 38.6 55

434206121312901 Burgess 2.1 22.00S/10.00E-03BBB01 DESC 52688 43°42'05.67" 121°31'28.92" 4,239.2 23.1 17

434206121312902 Burgess 2.2 22.00S/10.00E-03BBB02 DESC 52686 43°42'05.67" 121°31'28.92" 4,239.0 30.0 23

434206121312903 Burgess 2.3 22.00S/10.00E-03BBB03 DESC 52687 43°42'05.67" 121°31'28.92" 4,239.0 36.3 –

434206121312904 Burgess 2.4 22.00S/10.00E-03BBB04 DESC 52685 43°42'05.67" 121°31'28.92" 4,239.0 52.3 7.7

434213121295901 Burgess 6.1 21.00S/10.00E-35CDC1 DESC 53163 43°42'12.64" 121°29'59.34" 4,198.6 14.7 9.0

434213121295902 Burgess 6.2 21.00S/10.00E-35CDC2 DESC 53159 43°42'12.64" 121°29'59.34" 4,198.7 45.0 1.3

434213121324101 Burgess 1.1 22.00S/10.00E-05AAA01 DESC 52690 43°42'13.02" 121°32'41.42" 4,257.7 34.7 64

434213121324102 Burgess 1.2 22.00S/10.00E-05AAA02 DESC 52689 43°42'13.02" 121°32'41.42" 4,257.6 50.5 33

434910121275501 Century 1.1 20.00S/11.00E-19CCC1 DESC 53166 43°49'10.08" 121°27'55.27" 4,178.2 12.1 12

434910121275502 Century 1.2 20.00S/11.00E-19CCC2 DESC 53178 43°49'10.08" 121°27'55.27" 4,178.0 68.3 .3

434945121273501 Century 2 20.00S/11.00E-19BDB DESC 53184 43°49'44.61" 121°27'35.32" 4,175.4 11.1 8.8

435007121272701 Century 3 20.00S/11.00E-18CDC1 DESC 53175 43°50'07.04" 121°27'27.31" 4,173.9 8.9 97

435026121272101 Century 4 20.00S/11.00E-18CAA DESC 53185 43°50'25.67" 121°27'20.64" 4,174.1 10.7 57

435030121271401 Century 5 20.00S/11.00E-18ACC1 DESC 53160 43°50'29.75" 121°27'14.15" 4,172.9 11.7 5.2

435035121270801 Century 6.1 20.00S/11.00E-18ACA1 DESC 53164 43°50'34.68" 121°27'08.15" 4,171.9 10.7 31

435035121270802 Century 6.2 20.00S/11.00E-18ACA2 DESC 53158 43°50'34.68" 121°27'08.15" 4,171.8 64.0 .8
12 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

The wells ranged from 8.9 to 68 ft deep and were completed 

with 2-ft long screens. Falling- and rising-head tests were 

run with initial-head displacements of 1–2 ft. Test data were 

analyzed using the method of Bower and Rice (1976), with 

the assumptions that the effects of storage can be ignored and 

the change in saturated thickness of the aquifer is negligible. 

Horizontal hydraulic conductivities computed using the slug 

test data range from 0.3 to 98 ft/d, with a geometric mean of 

25 ft/d and median of 15 ft/d (fig. 4, table 2). The 10th and 

90th percentile values of hydraulic conductivity (0.6 and 

60 ft/d) span only two orders of magnitude which, even when 

considering that the tests selectively target coarser materials, 

suggest that hydraulic conductivity is relatively uniform in the 

area compared to some alluvial environments. The 10th–90th 

percentile range in hydraulic conductivity for the well-yield 

test data was 4 to 94 ft/d. The 10th percentile hydraulic 

conductivity is greater for this data set because the slug tests 

are conducted on water-supply wells selectively screened 

to the most permeable sediments and more thoroughly 

developed.

Median values from well-yield and slug test data were 

similar (21 and 15 ft/d, respectively). Geometric mean 

hydraulic conductivity values were greater for the well-yield 

test data (60 ft/d) because several basalt wells included large 

(> 100 ft/d) hydraulic conductivity values that skewed the 

mean (fig. 4). The agreement between the median values from 

both data sets indicates that both estimation methods provide 

good information on the hydraulic characteristics of the basinfill sediments. 

Ground-Water Recharge

The shallow alluvial ground-water system in the study 

area is recharged primarily by infiltration of precipitation 

(rainfall and snowmelt), with lesser amounts resulting 

from lateral ground-water inflow and effluent from on-site 

wastewater systems. Gannett and others (2001) estimated 

ground-water recharge to the upper Deschutes River basin 

from infiltration of precipitation using a water balance model 

(Deep Percolation Model) developed by Bauer and Vaccaro 

(1987). The water balance model is based on empirical 

relations that quantify processes such as interception and 

evaporation, snow accumulation and melt, plant transpiration, 

and runoff. The model computes a complete daily water 

balance for the soil zone using measured precipitation and 

temperature and data describing land cover, vegetation, and 

soil properties. The water balance was computed within 

rectangular areas (cells), each with an area of 1.3 mi2. A 

detailed description of application of the water balance model 

to the upper Deschutes River basin, including model input, is 

available in Boyd (1996) and a summary of the results for the 

entire upper basin is given in Gannett and others (2001). 

Computed recharge from the water balance model was 

used to specify the initial spatial distribution of recharge 

within the boundaries of the simulation model. Based on 

calibration of simulation models in the La Pine area (see 

“Nitrate Fate and Transport Simulation Models”), the 

recharge distribution was reduced slightly; the final calibrated 

distribution is shown in figure 5. The mean annual recharge 

distribution ranges from about 2 to 20 in/yr with the greatest 

recharge occurring in the upland areas; however, most of the 

model area lies at lower elevations adjacent to the Deschutes 

and Little Deschutes Rivers where mean annual recharge is 

about 2 in/yr (fig. 5). The mean annual recharge rate for the 

simulation model is 3.2 in/yr (58 ft3

/s). 

100

PERCENTAGE OF OBSERVATIONS WITH 

HYDRAULIC CONDUCTIVITY LESS THAN INDICATED VALUE

0 10 20 30 40 50 60 70 80 90

HYDRAULIC CONDUCTIVITY. IN FEET PER DAY

0.1

 1

10

100

 1

10

100

1,000

10,000

Slug test data

 (n=24)

Well yield test data

 (n=221)

Figure 4. Distributions of estimated hydraulic 

conductivity values from slug tests and well-yield data.
Conceptual Model of Ground-Water System 13

Figure 5. Distribution of mean annual recharge, water-table elevation contours (June 2000), and locations of monitoring wells in 

the La Pine, Oregon, study area. 

OR19-0049_fig 05

Base 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

Publication projection is 

Lambert Conformal Conic

Standard parallels 42o20' and 44o40', 

central meridian -120o30'

Datum is NAD83

R 9 E R 10 E R 11 E

T 21 S

T 22 S

T 23 S

T 24 S

T 20 S

T 19 S

DESCHUTES CO.

L A K E

KLAMATH CO.

Spring

River

Fall

River

Crescent Creek

Creek

Little

Deschutes

DESCHUTES

RIVER

River

Paulina

Long Creek

5 KILOMETERS

0 5 MILES

0

Wickiup

Reservoir

97

97

31

31

Burgess Road

Drive

Finley

Butte

Road

Road

South

Century

Spring

River

Masten Road

SunriverSunriver

La 

Pine

43°55'

50'

45'

121°40' 35' 30' 25' 121°20'

40'

43°35'

Extent of study-area model

EXPLANATION

Mean annual ground-water 

recharge, in inches per 

year 

1-2 

2-4 

4-6 

6-8

8-12

12-20

Line of equal water table 

elevation, in feet above 

North American Vertical 

Datum of 1929

Well location and 

measurement frequency

2-hour

Bimonthly 

4220

4210 4220

4170

4230

4190

4180

4240

4160

4260

4210

4220

4170

4230

4190

4180

4240

4160

4260

4200

4250

4200

4250

Deschutes

Basin

Cascade

Range

Newberry Upper Volcano Deschutes

Basin

La Pine

region

OREGON

Bend
14 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

Effluent from on-site wastewater systems also contributes 

to recharge of the shallow ground-water system. The annual 

rate was estimated for 1999 when there were about 5,200 

on-site wastewater systems in the study area. The average 

daily volume of effluent produced per system depends on the 

number of occupants and their water-use habits. Deschutes 

County monitored several on-site systems as part of the La 

Pine National On-Site Demonstration Project and determined 

that the average daily effluent volume was 45 gallons per day 

per person (B. Rich, written commun., Deschutes County, 

2003). La Pine households averaged 2.55 persons in 2000 

(U.S. Bureau of the Census, 2000), so each household would 

produce an average of 115 gal/d. In 1999, about 5,200 on-site 

systems would have contributed 0.6 Mgal/d, (0.9 ft3/s), of 

recharge to the ground-water system. 

Ground-Water Flow

Hydraulic head, or simply head, is a measure of the force 

that drives ground-water movement. Ground water flows from 

areas of high head to areas of low head. In an unconfined

aquifer, such as the shallow alluvial deposits in the La Pine 

study area, the elevation of the water table represents the 

head at the upper surface of the aquifer. The rate of change in 

head with distance is called the hydraulic gradient, which is 

the slope of the water table. Ground water in deeper aquifers 

may be confined by lower permeability layers (for example, 

silt or clay). Ground water in confined aquifers may be under 

pressure, such that when a well penetrates the aquifer, ground 

water will rise in the well casing to levels above the top of 

the aquifer. The hydraulic gradient in the vertical direction is 

positive if head increases with depth, which indicates ground 

water is flowing upward; conversely, it is negative if head 

decreases with depth, and ground water is flowing downward.

Water-level measurements were made in 192 wells in 

June 2000. Water levels from 170 wells with screened intervals 

in the upper 100 ft were used to construct a contour map of 

the water-table surface (fig. 5) used to calibrate the simulation 

model. Ground water in the shallow part of the system in the 

La Pine study area generally flows toward the rivers, which 

generally are gaining in the study area. Water-table contours 

also indicate that shallow ground water moves toward lowlying Long Prairie south of the La Pine core area (fig. 1). 

The water table is typically within 20 ft of land surface 

in the study area. The shallowest depths to water in the study 

area are along the lower reach of the Deschutes River and the 

entire Little Deschutes River. In many areas, the water table is 

within 5 ft of land surface and seasonally inundates low lying 

areas adjacent to the streams during high flows in the spring. 

Ground-water levels vary with time in response to 

changes in rates of recharge to and discharge from the groundwater system. To characterize seasonal and shorter variations, 

bi-monthly water-level measurements were made in 48 wells 

and digital recorders measured water levels every 2 hours 

in 9 additional wells. Most wells were monitored between 

March 2000 and December 2001; however, several wells in 

the bi-monthly network also were monitored between 1994 

and 1998 as part of a previous study of the upper Deschutes 

basin (Gannett and others, 2001). Two long-term observation 

wells also are in the study area and are measured by the 

Oregon Water Resources Department (OWRD). These wells, 

one of which has been measured quarterly since 1945, provide 

insights on the effects of withdrawals (pumping), decadalscale climate variation, and other influences. 

Three well hydrographs are shown in figure 6 to illustrate 

general observations regarding the response of ground-water 

levels to seasonal and long-term variation in recharge. Basic 

information for all wells included in the monitoring network 

is listed in table 3. Additional information for these wells, 

including the complete history of water-level data, is available 

from the USGS National Water Information System (NWIS) at 

http://waterdata.usgs.gov/or/nwis. 

USGS Site ID 434400121275801

21S/11E-19ccc 

OWRD Well ID DESC 7620

CALENDAR YEAR

CALENDAR YEAR

1945 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 

10

15

20

25

30

35

40

USGS Site ID 433231121360301

23S/09E-36bbc 

OWRD Well ID KLAM 136

30

35

40

USGS Site ID 434328121304601

21S/10E-27dbb1

OWRD Well ID DESC 53655

1994 1995 1996 1997 1998 1999 2000 2001 

WATER LEVEL, IN FEET BELOW LAND SURFACE

 5

10

15

Depth: 100 feet 

(deepened from 37 feet in 1964)

Open interval: 70-100 feet

Depth: 21 feet

Open interval: 18-21 feet

Depth: 470 feet

Open interval: 145-470 feet

Figure 6. Water levels in selected wells in the La Pine, 

Oregon, study area.
Conceptual Model of Ground-Water System 15

Recharge to the ground-water system from infiltration 

of precipitation and snowmelt occurs primarily in winter and 

early spring when evaporation and transpiration losses are at 

a minimum and available moisture is at a maximum. When 

recharge to the ground-water system exceeds discharge, 

ground water is added to storage and water levels rise. The 

response of hydraulic head to the annual recharge pulse in 

the winter and early spring is attenuated in deeper wells, such 

as KLAM 136 (fig. 6); however, the annual recharge pulse 

causes seasonal water-level increases of as much as 5 ft in 

shallow wells like DESC 53655 (fig. 6). The effects of longerterm variation in climate are evident in shallow and deep 

wells. Shallow (DESC 53655) and deep (KLAM 136) wells 

responded to the wet years from 1996 to 1998, although the 

response in the deep well was delayed and more gradual. The 

long-term monitoring well DESC 7620 is relatively shallow 

(100 ft) and shows that several wet and dry periods have 

occurred in the basin since the mid-1940s (fig. 6). The dry 

periods have an average length of 11 years. Based on the longterm hydrograph at well DESC 7620, water levels measured at 

wells in June 2000 and used to construct the water-table map 

in figure 5 are representative of long-term (1945–2005) mean 

water levels in the shallow part of the ground-water system. 

Ground-Water Discharge

Ground water discharges to streams, springs, and wells, 

and by evapotranspiration. Discharge to streams and springs 

is the primary pathway for ground water to leave the shallow 

aquifers in the study area. Discharge to streams occurs as 

diffuse seepage through the streambed (below the water line), 

as discrete springs or spring complexes that form tributaries 

(for example, Fall River), and as seepage along the stream 

bank (above the water line). Evapotranspiration, which 

includes bare soil evaporation as well as transpiration by 

plants, is another mechanism for discharge where the water 

table is shallow. Pumping by wells also accounts for part of 

the overall discharge from the shallow aquifers in the study 

area. 

Streams and Springs

The direction and rate of flow between ground water and 

streams is directly related to the direction and magnitude of 

the hydraulic gradient between the stream and the groundwater system. For ground water to discharge to a stream, 

the ground-water level must be greater than the stage of the 

stream. The rate of discharge to the stream is proportional to 

the difference between the ground-water level and the stream 

stage.

Table 3. Basin information for water-level monitoring wells in the 

La Pine, Oregon, study area.

[Abbreviations: USGS, U.S. Geological Survey. Frequency: B, bi-monthly; R, 

recorder (every 2 hours)]

State 

identification 

No.

USGS site 

identification No. Local No. Frequency

DESC 52676 434202121301503 22.00S/10.00E-02BBC03 R

DESC 52677 434202121301502 22.00S/10.00E-02BBC02 R

DESC 52678 434202121301501 22.00S/10.00E-02BBC01 R

DESC 52684 434206121311701 22.00S/10.00E-03BBD01 R

DESC 52690 434213121324101 22.00S/10.00E-05AAA01 R

DESC 53158 435035121270802 20.00S/11.00E-18ACA2 R

DESC 53164 435035121270801 20.00S/11.00E-18ACA1 R

DESC 53656 434328121304603 21.00S/10.00E-27DBB3 R

DESC 53660 434657121274101 21.00S/11.00E-06CBA R

DESC 937 433756121335501 22S/10E-30DDD B 

DESC 968 434810121305201 20S/10E-34ABC B 

DESC 1450 433850121302301 22S/10E-22DDD B 

DESC 1529 434700121285701 21S/10E-01BCD B 

DESC 1727 434731121341001 20S/10E-31DDC B 

DESC 6593 434516121295101 21.00S/10.00E-14CDD B 

DESC 6959 434400121291401 21S/10E-23DDD B 

DESC 7618 434501121272401 21S/11E-18CDA3 B 

DESC 7620 434400121275801 21S/11E-19CCC B 

DESC 7740 434217121263801 21.00S/11.00E-32CCC B 

DESC 9173 433948121300101 22S/10E-14CCA B 

DESC 9303 434752121343101 20.00S/10.00E-31DBB B 

DESC 9314 434534121323801 21.00S/10.00E-16BBB B 

DESC 9655 433959121300801 22S/10E-14CBB02 B 

DESC 9889 434041121330701 22.00S/10.00E-08DCC B 

DESC 50418 435104121282301 20.00S/10.00E-12DCB B 

DESC 52533 433847121302301 22.00S/10.00E-22DAA2 B 

DESC 52679 434203121305402 22.00S/10.00E-03ABC02 B 

DESC 52680 434203121305403 22.00S/10.00E-03ABC03 B 

DESC 52681 434203121305401 22.00S/10.00E-03ABC01 B 

DESC 52682 434206121311702 22.00S/10.00E-03BBD02 B 

DESC 52683 434206121311703 22.00S/10.00E-03BBD03 B 

DESC 52685 434206121312904 22.00S/10.00E-03BBB04 B 

DESC 52686 434206121312902 22.00S/10.00E-03BBB02 B 

DESC 52687 434206121312903 22.00S/10.00E-03BBB03 B 

DESC 52688 434206121312901 22.00S/10.00E-03BBB01 B 

DESC 52689 434213121324102 22.00S/10.00E-05AAA02 B 

DESC 53159 434213121295902 21.00S/10.00E-35CDC2 B 

DESC 53160 435030121271401 20.00S/11.00E-18ACC1 B 

DESC 53163 434213121295901 21.00S/10.00E-35CDC1 B 

DESC 53166 434910121275501 20.00S/11.00E-19CCC1 B 

DESC 53174 435007121272702 20.00S/11.00E-18CDC2 B 

DESC 53175 435007121272701 20.00S/11.00E-18CDC1 B 

DESC 53177 435030121271402 20.00S/11.00E-18ACC2 B 

DESC 53178 434910121275502 20.00S/11.00E-19CCC2 B 

DESC 53184 434945121273501 20.00S/11.00E-19BDB B 

DESC 53185 435026121272101 20.00S/11.00E-18CAA B 

DESC 53651 434029121292401 22.00S/10.00E-11DCC B 

DESC 53652 434026121314301 22.00S/10.00E-16AAB1 B 

DESC 53653 434026121314302 22.00S/10.00E-16AAB2 B 

DESC 53655 434328121304601 21.00S/10.00E-27DBB1 B 

DESC 53657 434328121304602 21.00S/10.00E-27DBB2 B 

DESC 53658 434212121335301 22.00S/10.00E-06AAB B 

DESC 53659 433703121305101 22.00S/10.00E-34DCC B 

DESC 53661 434226121282302 21.00S/10.00E-36ACD2 B 

KLAM 136 433231121360301 23S/09E-36BBC B 

KLAM 138 433152121281301 23S/10E-36DDC B 

KLAM 51668 433625121354401 23.00S/09.00E-01DAC B 
16 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

Ground water discharge to streams and springs was 

estimated from gain-loss surveys of the Deschutes and Little 

Deschutes Rivers between October 1995 and October 2000. 

Gain-loss surveys, sometimes termed “seepage studies”, are 

used to measure increases or decreases in streamflow between 

upstream and downstream measurement sites that can be 

attributed to exchange of water between the stream channel 

and the underlying ground-water system. The gain-loss 

surveys were conducted by measuring streamflow at intervals 

of 1 to 10 mi on the Deschutes River between Wickiup Dam 

(river mile [RM] 226.7) and Benham Falls (RM 181.4) and 

on the Little Deschutes River between Wagon Trail Ranch 

(RM 43.0) and immediately upstream of the mouth (RM 

0.5). Gaging stations are located below Wickiup Dam and 

at Benham Falls on the Deschutes River and at RM 26.7 on 

the Little Deschutes River (fig. 7). Tributary inflow to the 

Deschutes River was measured from Fall, Spring, and Little 

Deschutes Rivers, and tributary inflow to the Little Deschutes 

River was measured from Paulina Creek and Long Creek. The 

gain or loss to or from each reach was computed by summing 

the upstream and tributary inflow to the reach and deducting 

the downstream outflow and diversions. Any positive residual 

was assumed to be discharge from ground water; negative 

residuals were assumed to represent stream leakage (recharge) 

to ground water. Streamflow measurements and estimated 

gains and losses for surveys in October 1995, February 1996, 

March 2000, and October 2000, are listed in table 4.

The gain-loss surveys of the Deschutes River show 

net gains of 143 and 166 ft3

/s in 1996 and 2000. Most gains 

occur downstream of Harper’s Bridge, RM 191.7, where large 

springs discharge near river level from basalt aquifers. Spring 

discharge to this reach was about 170–180 ft3/s (15–20 percent 

of flow) in the two gain-loss surveys. Between Wickiup Dam 

(RM 226.7) and Harper’s Bridge the 1996 and 2000 surveys 

show a possible net loss of 10–40 ft3/s. Between Harper’s 

Bridge and General Patch Bridge, surveys showed net losses 

of 54 and 24 ft3/s in 1996 and 2000, respectively. Although 

the seepage results show apparent losses in this reach they are 

similar in magnitude to the estimated error of measurement 

(32–38 ft3

/s). Other data (contours showing horizontal 

hydraulic head gradients toward the river and upward vertical 

gradients measured below the river bed) indicate that ground 

water discharges to the river in this reach and seepage data 

lack the accuracy to measure the gain in flow. Upstream of 

General Patch Bridge, the Deschutes River receives groundwater discharge of as much as 25 ft3

/s, with most entering 

the river downstream of the La Pine State Recreation Area. 

The 10-mi reach upstream of Pringle Falls showed little or no 

ground-water interaction in 1996 and 2000.

Overall, little net gain or loss was measured in flow of 

the Deschutes River upstream of Harper’s Bridge, relative 

to the total discharge of the river. This is probably due to 

the low permeability of the ash and tephra deposits that 

constitute much of the streambed and lack of major springs 

discharging directly to the channel in this reach. However, 

exchange between the river and the ground-water system may 

be occurring at a smaller scale that cannot be resolved by 

gain-loss surveys in which measurements are several miles 

apart. The multitude of meander loops likely are hydraulically 

connected by old channel deposits that may provide groundwater flow pathways where stream water can leave the channel 

and re-enter downstream. 

The Little Deschutes River downstream of RM 43 

showed net gains ranging from 6.7 to 19.8 ft3

/s in the gain-loss 

surveys conducted in October 1995, March 2000, and October 

2000. All three surveys showed that the reach between Bridge 

Drive (RM 15.9) and South Century Drive (RM 5.5) received 

ground-water discharge; the gains in the reach ranged from 

5.9 to 25.2 ft3/s (0.6 to 2.4 ft3/s/mi), although the measurement 

error was nearly equal to the apparent 25.2 ft3

/s gain in March 

2000. Gains ranging from 2.9 to 28.5 ft3/s (0.2 to 1.8 ft3/s/mi) 

also were measured in the 16.3 mi reach upstream of the 

gaging station (at RM 26.7). The reach from South Century 

Drive (RM 5.5) to just upstream of the mouth (RM 0.5) had 

little or no gain or loss in the October surveys (1995 and 

2000); however, the March 2000 survey indicated a 35 ft3/s 

loss in this reach. 

Long Creek is an intermittent tributary to the Little 

Deschutes that drains the area between Highways 97 and 31 to 

the south of La Pine. This topographically low area probably 

is an abandoned channel of the Little Deschutes River. Aerial 

photography showing lush vegetation throughout summer 

indicates that this area has a shallow water table and likely is 

a ground-water discharge area. In a typical year, Long Creek 

flows until the water-table declines below the streambed in 

June–July and then remains dry until the water table rises in 

response to winter-spring recharge.

The October surveys probably were representative of 

average conditions (neutral or slight gains) on the lowermost 

reach of the Little Deschutes River, whereas the March 

2000 survey probably was affected by high flows (about 

300 ft3

/s) that occurred during the seepage measurements. 

The area between the Deschutes River and Little Deschutes 

immediately upstream of their confluence is underlain by 

10–15 ft of pumice, sand, and fine gravel that overlie a dense, 

low-permeability clay. The top of the clay is exposed in the 

banks of the Deschutes and Little Deschutes Rivers a few feet 

above the water surface at low flow. In late spring, groundwater seepage is visible on the bank at the contact between the 

sand and gravel and the clay. At high stream stage, as during 

the March 2000 survey, surface water would move from the 

channel into the exposed sand and gravel while the hydraulic 

gradient between the stream and ground water is temporarily 

reversed (stream stage is greater than the ground-water level). 
Conceptual Model of Ground-Water System 17

Figure 7. Gaining and losing reaches of the Deschutes and Little Deschutes Rivers and locations of measurement sites for 

gain-loss surveys between October 1995 and October 2000 in the La Pine, Oregon, study area. OR19-0049_fig 07

Base 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

Publication projection is 

Lambert Conformal Conic

Standard parallels 42o20' and 44o40', 

central meridian -120o30'

Datum is NAD83

R 9 E R 10 E R 11 E

T 21 S

T 22 S

T 23 S

T 24 S

T 20 S

T 19 S

DESCHUTES

L A K E

KLAMATH

Spring

River

Fall

River

Crescent Creek

Creek

Little

Deschutes

DESCHUTES

RIVER

River

Paulina

Long Creek

5 KILOMETERS

0 5 MILES

0

Wickiup

Reservoir

97

97

31

31

Burgess Road

Drive

Finley

Butte

Road

Road

South

Century

Spring

River

Masten Road

SunriverSunriver

La 

Pine

43°55'

50'

45'

121°40' 35' 30' 25' 121°20'

40'

43°35'

Extent of study area model

0-35

24-54

1-2

11-24

1-6

170-180

0-12

3-29

6-25

RM 13.7

14063000

14056500

14057500

14064500

RM43.0

RM26.7

RM226.7

RM

217.6

RM

208.6

RM

199.7

RM

192.7

RM

191.7

RM 181.4

RM

5.5

RM

0.5

EXPLANATION

Gain or loss for stream 

reach—Number is gain 

or loss in cubic feet per 

second

Gain—more certain

Gain—less certain

Loss—less certain

Loss—more certain

Gaging station with USGS 

site number—See table 4 

for details

Discharge measurement 

site with river mile—See 

table 4 for details

6-25

RM43.0

14063000

Long Prairie
18 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon Table 4. Summary of data from stream gain-loss surveys on the Deschutes and Little Deschutes Rivers in the La Pine, Oregon, study area, 1995–2000. [All values are in cubic feet per second unless otherwise noted. Q is measured discharge in the mainstem river; Tributary Q is measured tributary inflow above the indicated river mile; Seepage is the unmeasured gain (+) or loss (-) above the indicated river mile; Error is the estimated error measurement (+/-); Standard deviation of measurements is ± 5 percent; Error is computed as square root of the sum of the variances for measurements used to compute seepage. Data sources: source 1 is from K. Gorman, Oregon Water Resources Department, Bend, OR, written commun.; source 2 is from K. Lee, U.S. Geological Survey, written commun. Data for gaging stations at Benham Falls and near La Pine after 1991 are from Bureau of Reclamation. Abbreviation: USGS, U.S. Geological Survey] Measurement site (USGS station no.) River mile

Q Tributary 

Q

Estimated 

seepage

Error

Q Tributary 

Q

Apparent 

seepage

Error

Q Tributary 

Q

Estimated 

seepage

Error

Deschutes River

02-21-1996 (source: 1) 10-24-2000 (source: 1)

Below Wickiup Dam (14056500) 226.7 128 274

At Pringle Falls 217.6 126 -2 9 273 -1 19

At Lapine State Recreation Road 208.6 132 6 9 274 1 19

At General Patch Bridge 199.7 268 112 24 16 427 142 11 26

Above Little Deschutes 192.7 394 -33 29

At Harper’s Bridge 191.7 758 544 -54 38 490 87.0 9 32

At Benham Falls (14064500) 181.4 1,050 120 172 65 881 212 179 52

Total 45.3 146 53 166 46

Little Deschutes River

10-10-1995 (source: 1) 03-20-2000 (source: 2) 10-23-2000 (source: 1)

Wagon Trail Ranch 43.0 50.4 70.7

Masten Road 37.9 275

Deer Road 33.5 53.3 2.9 3.7

Dorance Meadow Road 32.6 71.4 0.7 5

Near La Pine (14063000) 26.7 53.3 0 3.8 308 4.50 29 21 81.2 9.8 5

Bridge Drive 15.9 296 -12 21

State Recreation Area Road 13.7 55.9 2.6 3.9 79.6 -1.6 6

South Century Drive 5.5 59.2 3.3 4.1 324 2.80 25 22 89.4 9.8 6

At Crosswater Golf Course bridge 1.4 57.4 -1.8 4.1

At Crosswater Golf Course bridge 

(mouth)

0.5 289 -35 22 90.5 1.1 6

Total 42.5 7.0 3.8 7 20 19.8 6
Conceptual Model of Ground-Water System 19

The surveys on the Little Deschutes River in October 

1995 and October 2000 are surprisingly dissimilar considering 

that the discharge and measurement error were relatively 

low for each survey. The apparent net gain in October 2000 

was nearly three times that measured in October 1995. The 

difference between the results of the two surveys likely is 

explained by the water-table elevation at the time of each 

survey. In 1995, the basin had undergone several years of 

less than-normal precipitation and ground-water levels had 

fallen by nearly 20 ft over the previous decade in long-term 

observation well DESC 7620 (fig. 6). The lowering of the 

water table reduced the hydraulic gradient between the shallow 

ground water and the stream and, therefore, the ground-water 

discharge to the stream also was reduced. Between 1995 and 

2000, ground-water levels recovered almost 15 ft, and based 

on the gain-loss surveys, ground-water discharge to the Little 

Deschutes River increased from 7 to nearly 20 ft3/s.

In October–November 2000, a survey was made of the 

vertical hydraulic head gradient between the shallow groundwater system and the Deschutes and Little Deschutes Rivers. 

The purpose of the survey was to identify the locations 

of gaining and losing reaches of the Little Deschutes and 

Deschutes Rivers at a more detailed scale than possible using 

the gain-loss surveys in which discharge measurements were 

made at locations separated by as much as 10 mi. Hydraulic 

gradient alone cannot be used to determine the volumetric 

rate of flow between a stream and the aquifer. Stream-channel 

profiles and bed lithology data (appendix A) also were 

collected during the survey and later used to estimate the 

stream bed conductance parameters required by the simulation 

model. The simulation model was used to estimate the rate 

of ground-water flow into or out of the stream; the simulated 

rates were compared with estimates from the seepage surveys 

during model calibration. Head gradient measurements were 

made at 12 sites on the Deschutes River between RM 192.5 

and RM 218.4 and 20 sites on the Little Deschutes River 

between RM 2.0 and RM 42.6 (appendix A). The average 

distance between head gradient measurements was 2 river 

miles. Measurements were made using a portable hydraulic 

potentiometer that was driven 1–3 ft into the streambed. The 

relative head difference between the stream and ground water 

was measured to the nearest 0.01 ft using a manometer. The 

design of the potentiometer and details on the procedure 

for measuring head gradients is given in Winter and others 

(1988). The dimensionless vertical hydraulic head gradient 

was computed at each location by dividing the head difference 

(ground-water level minus stream stage) by the depth of the 

center of the potentiometer screen below the stream bed. 

Ground-water discharge was indicated at 18 of 20 

measurement sites on the Little Deschutes River and at all 12 

sites on the Deschutes River (appendix A). Vertical hydraulic 

gradients ranged from -0.003 to 0.075 on the Little Deschutes 

River (negative gradients indicate downward flow from the 

river) and 0.010 to 0.060 on the Deschutes River. All gradients 

were upward (indicating ground-water discharge to the river) 

except for two sites on the Little Deschutes River (RM 28.1 

and 29.5; fig. A1, table A1) where neutral or slight downward 

gradients were measured. 

Evapotranspiration

Ground-water discharge to the atmosphere can occur 

by evaporation from bare soil and transpiration from the 

leaves of phreatophytes (plants whose roots draw water from 

the saturated zone). The rate of evaporation from bare soil 

diminishes rapidly with increasing depth to the water table 

and is negligible if the water table is more than 10 ft below 

land surface (Brutsaert, 1982, p. 236). Transpiration rates are 

dependent on the type and density of phreatophytes, climatic 

conditions, quality of water, and depth to water (Robinson, 

1958, p. 16). These processes usually are considered together 

and referred to as evapotranspiration, or ET. 

In June 2000, the water table was within 10 ft of land 

surface over an area of about 22 mi2

. The potential ET (PET) 

rate is the maximum possible rate for the local climate 

conditions if water is nonlimiting. As part of a regional 

ground-water budget, Gannett and others (2001) estimated that 

the PET for the saturated zone in the La Pine study area was 

22 in/yr, which is equivalent to about 1.6 ft3/s/mi2. This rate 

would occur only if the water table was at or near land surface 

with full coverage of phreatophytes. Under these conditions, 

the maximum rate of ground-water discharge from the study 

area by ET would be 35 ft3/s. Because the water table is not at 

land surface and the coverage of phreatophytes is considerably 

less than 100 percent, a more reasonable estimate of total ET 

from the study area is in the range of 10 to 20 ft3/s.

Wells

Eleven public-water supply systems (Oregon Department 

of Human Services, 2005) served about 1,650 people in the 

study area in 2000. The largest water purveyor, the La Pine 

Water District, serves 800 residents in the central core area of 

La Pine. A private water system serves about 350 residents 

of neighborhoods north of South Century Drive, between the 

Deschutes and Little Deschutes Rivers. Several other small 

water systems serve about 500 residents of mobile home and 

recreational vehicle parks scattered throughout the study area. 

About 12,350 residents, about 90 percent of the population, 

obtained their water from individual, privately-owned wells in 

2000. Local water use data for self-supplied rural residents in 

the study area was not available. Mean per capita withdrawals 

for self-supplied rural domestic use in Oregon is estimated 

to be 110 gal/d (Hutton and others, 2004); however, given 

the limited amount of irrigated landscaping in the La Pine 
20 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

area, this value probably is not applicable. Deschutes County 

monitored wastewater effluent discharge at homes in the area 

as part of the La Pine Demonstration Project and determined 

the average per capita rate was 45 gal/d (Oregon Department 

of Environmental Quality, 2004b). To estimate domestic 

well withdrawals, it was assumed that outside water use was 

equal to 25 percent of total use, making the total per capita 

withdrawal rate equal to 60 gal/d. Assuming this rate for 

each of the 14,000 residents, about 0.84 Mgal/d (1.3 ft3/s) 

were withdrawn in 2000. Infiltration of effluent from on-site 

wastewater systems returns is about 45 gal/d per capita 

(0.63 Mgal/d total) to the ground-water system as recharge at 

the water table, resulting in net withdrawals of 0.21 Mgal/d 

(0.32 ft3/s). 

Nitrogen Fate and Transport

As part of this study, Hinkle and others (2007a) applied 

geochemical and isotopic tools at various scales to provide 

a framework for understanding aquifer-scale nitrate source, 

transport, and fate. (Other wastewater contaminants, such as 

viruses and pharmaceuticals, were evaluated by Hinkle and 

others [2005].) The conceptual model resulting from this work 

was used to develop the transport simulation model.

Nitrogen concentration data, tracer-based apparent 

ground-water ages, ratios of N/Cl-, N isotope data, and 

ground-water flow directions indicate that on-site wastewater 

effluent is the only significant anthropogenic nitrogen source 

to shallow ground water in the area. High concentrations of 

ammonium were measured in samples from deep ground 

water. Nitrogen isotopes, N/Cl- and N/C ratios, tritiumage data, and hydraulic-head gradients support a natural, 

sedimentary organic matter source for the high concentrations 

of ammonium in the deep aquifer, as opposed to an origin 

from on-site wastewater derived nitrogen. 

Most residential development in the La Pine area 

has occurred since 1960, with accelerated growth during 

1990–2005. As a result, loading of nitrate from on-site 

wastewater systems is a relatively recent phenomenon that, 

in combination with low ground-water recharge rates, has 

resulted in high concentrations of nitrate near the water 

table. Low recharge rates and flow velocities have, for now, 

generally restricted nitrate occurrence to discrete plumes 

within 20–30 ft of the water table. Concentrations of nitrate 

typically are low in deeper, older ground water due to the 

nature and timing of nitrate loading and transport, and to loss 

by denitrification. Denitrification is the process of reducing

nitrate into gaseous nitrogen. The process is performed by 

bacteria when dissolved oxygen (which is a more favorable 

electron acceptor) is depleted, and bacteria turn to nitrate in 

order to oxidize organic carbon or other electron donors. 

Denitrification was identified as an important 

geochemical process in the study area by Hinkle and others 

(2007a), who reported various supporting evidence, including 

nitrate/chloride relations, presence of excess N2 enriched in 15N, presence of N2

O, and nitrate losses within a framework 

of progressively reduced ground water, reflecting classical 

geochemical evolution. Ground water in the study area evolves 

from oxic to increasingly reduced conditions with increasing 

depth below the water table. Suboxic conditions are achieved 

in 15–30 years, and the boundary between oxic and suboxic 

ground water is sharp. Nitrate is denitrified near the oxicsuboxic boundary. 

The denitrification process was simulated by specifying 

the location of the oxic-suboxic boundary within the groundwater system and applying the assumption that nitrate is 

rapidly converted to nitrogen gas as it is transported across 

the boundary. The location of the boundary was mapped 

within the study-area using data from 256 wells where 

dissolved oxygen concentrations were available (Hinkle and 

others, 2007a). Oxic ground water was defined as water with 

dissolved oxygen concentrations greater than to 0.5 mg/L. The 

boundary between oxic and suboxic water was constrained at 

each sampling site using information on the depth to the water 

table and well construction. The elevation of the oxic/suboxic 

boundary cannot be exactly determined at wells; however, 

the elevation can be constrained if depth to the water table 

and the top of the screened interval of the well are known. 

The presence of suboxic water indicates that the boundary 

lies above the top of the screened interval and the presence of 

oxic water indicates that the boundary lies below the top of 

the screened interval. In the latter case, the well may yield a 

mixture of oxic and suboxic water if the boundary lies within 

the screened interval, or the water may be entirely oxic if the 

boundary lies below the screened interval. In areas where 

data were sparse, an understanding of general patterns of 

occurrence of oxic and suboxic water was used to infer the 

location of the boundary. The general absence of oxic water 

below a depth of 50 ft below the water table helped constrain 

the depth over large parts of the study area. Similarly, the 

general absence of suboxic water in ground water within 10 ft 

of the water table constrained the minimum thickness of the 

oxic part of the aquifer. However, one limitation of the data 

was the absence of wells in near-stream areas to constrain 

the boundary location. The thickness of the oxic part of the 

aquifer was contoured (fig. 8) and used to specify the location 

of the oxic-suboxic boundary in the transport model. The 

thickness of the oxic part of the aquifer decreases as ground 

water moves toward discharge areas along the Deschutes 

River, Little Deschutes River, and Long and Paulina Creeks. 

This is consistent with the aging of ground water as it moves 

along flow paths terminating at surface-water features in 

the area. Details of the implementation of the boundary are 

presented in the description of the simulation model.
Conceptual Model of Ground-Water System 21

Figure 8. Estimated thickness of oxic ground-water layer in the shallow aquifer in the La Pine, Oregon, study area.

OR19-0049_fig08

Base 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

Publication projection is 

Lambert Conformal Conic

Standard parallels 42o20' and 44o40', 

central meridian -120o30'

Datum is NAD83

R 9 E R 10 E R 11 E

T 21 S

T 22 S

T 23 S

T 24 S

T 20 S

T 19 S

DESCHUTES CO.

L A K E

KLAMATH CO.

EXPLANATION

Estimated thickness of the 

oxic layer—In feet, 

where oxic ground water 

is defined as water with a 

dissolved oxygen 

concentration greater 

than 0.5 milligrams per 

liter.

10-12 

20 

40

50

Spring

River

Fall

River

Crescent Creek

Creek

Little

Deschutes

DESCHUTES

RIVER

River

Paulina

Long Creek

5 KILOMETERS

0 5 MILES

0

Wickiup

Reservoir

97

97

31

31

Burgess Road

Drive

Finley

Butte

Road

Road

South

Century

Spring

River

Masten Road

SunriverSunriver

La 

Pine

43°55'

50'

45'

121°40' 35' 30' 25' 121°20'

40'

43°35'

Extent of study-area model
22 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

Nitrogen Loading from On-Site Wastewater 

Systems

The annual rate of nitrogen loading to the shallow 

aquifer system from on-site wastewater systems was 

estimated for 1960–2005. The simulation model was run 

for 1960–99 using the estimated annual loading rates as 

input and simulated nitrate concentrations were compared 

with measured concentrations. Previous investigations 

(Century West Engineering, 1982; Oregon Department of 

Environmental Quality, 1994) have concluded that the only 

significant anthropogenic source of nitrogen to the groundwater system in the La Pine study area is effluent from on-site 

wastewater systems. As previously described, Hinkle and 

others (2007a) also cite several lines of geochemical evidence 

supporting the conclusion that on-site wastewater systems are 

the dominant anthropogenic nitrogen source. This conclusion 

also is consistent with the low level of agricultural activity, 

limited number of livestock, and infrequent use of turf (and 

therefore fertilizer) in residential landscaping in the study 

area (Century West Engineering, 1982; Oregon Department of 

Environmental Quality, 1994; Hinkle and others, 2007a).

Annual nitrogen loading to ground water from on-site 

wastewater-system effluent was estimated for each tax lot in 

the study area annually from 1960 to 1999. The population 

of the study area prior to 1960 was sufficiently small that 

nitrogen loading was considered negligible. Nitrogen loading 

estimates were compiled by tax lot, but later aggregated 

across larger areas for input to the simulation model. 

Sources considered for nitrogen loading estimates included 

residential, schools, motels, restaurants, and recreational 

vehicle and trailer parks. Tax lot data (location, property 

class, improvement class, year built) were provided by 

Deschutes County (Tim Berg, Deschutes County Department 

of Community Development, written commun., 2001) 

and Klamath County (Jim McClellan, Klamath County 

Department of Management Information Services, written 

commun., 2001). Nitrogen loading was assumed to begin 

when lots were developed and end when lots were connected 

to a centralized sewer system (see areas served by centralized 

sewers in figure 1). Dates of sewer installation were provided 

by ODEQ (Dick Nichols, ODEQ, written commun., 2001) 

and La Pine Special Sewer District (Andy Newton, La Pine 

Special Sewer District, oral commun., 2001). Nitrogen loads 

were calculated for each of the approximately 5,200 lots that 

were developed and used on-site wastewater systems in 1999. 

Nitrogen loading for homes was assumed to equal the 

product of the average number of persons per household 

(2.55 persons/household in the study area [U.S. Bureau 

of the Census, 1960, 1970, 1980, 1990, 2000]); average 

volume of effluent discharged to the on-site system by 

each person (45 gal/person/d or 170 L/person/d [U.S. 

Environmental Protection Agency, 1980; Oregon Department 

of Environmental Quality, 2004b]), and average nitrogen 

concentration in on-site wastewater system effluent (61 mg 

N/L in the study area [B. Rich, ODEQ, written commun., 

2001]). Nitrogen discharged from on-site wastewater systems 

also was assumed to be completely oxidized to nitrate in 

the unsaturated zone beneath the drain field. Using these 

assumptions, the estimated potential average nitrogen loading 

per home is 21 lb N/yr (9.7 kg N/yr). 

The potential for nitrogen loss in the unsaturated zone 

was evaluated by sampling packed-bed (sand) filter on-site 

systems in the upper Deschutes Basin. Inflow to and outflow 

from the sand filters was sampled at 15 sites in two networks. 

The first network consisted of 10 sand filters, installed 0.9 to 

9.1 years prior to sampling. The second network consisted of 5 

new sand filters. The first network was sampled once (October 

2001) and the second network was sampled approximately 

bi-monthly for 3 years after installation (October 2000 – 

November 2003). Samples were analyzed by ODEQ as part 

of the La Pine NDP for nitrate, ammonium, total Kjeldahl 

nitrogen, and chloride. Nitrogen concentrations were adjusted 

for dilution and evaporation using chloride concentrations. 

The adjusted nitrogen concentrations from the first sand-filter 

network indicated median and mean nitrogen losses of 31 

percent and 29 percent, respectively (ODEQ, unpub. data, 

2004). Data from the second network showed that the sand 

filters underwent a period of maturation after installation. 

Using data from samples collected more than one year after 

installation, adjusted nitrogen concentrations from the second 

sand-filter network indicated median and mean nitrogen 

losses of 11 percent (Oregon Department of Environmental 

Quality, 2004b). The artificial unsaturated zone created by a 

sand-filter on-site system would be expected to function much 

like the natural unsaturated zone beneath standard on-site 

systems in the study area and similar rates of nitrogen loss 

would be expected beneath standard and sand-filter on-site 

wastewater systems. Based on the range of nitrogen loss 

indicated in the sand-filter sampling study (approximately 

10 to 30 percent), nitrogen losses in the unsaturated zone 

were assumed to average 25 percent for conventional on-site 

systems (including standard, pressure-distribution, and sandfilter systems) for computing loading rates. Nitrogen loss of 

25 percent in the unsaturated zone would reduce the assumed 

nitrogen concentration of effluent as it enters the saturated part 

of the ground-water system (water table) from 61 mg N/L to 

46 mg N/L. For the same effluent volume, the resulting annual 

nitrogen loading to the ground-water system would be reduced 

from 21 lb N/yr (9.7 kg N/yr) to 16 lb N/yr (7.3 kg N/yr). 

Nitrogen loading estimates were further adjusted to 

account for seasonal residency. For many years a significant 

number of homes in the study area were second homes or 

vacation rentals occupied only part of the year. In 1980, about 

46 percent of the population in the study area was seasonal 

(Century West Engineering, 1982). Seasonal residents made 

up an estimated 20 percent of study area population in 2000 

(Keven Bryan, La Pine Postmaster, oral commun., 2000; 

Carla Crume, La Pine Chamber of Commerce, oral commun., 

2000). To account for seasonal residents in the estimates of 

nitrate loading it was assumed that (1) a seasonal resident did 

not occupy their home for 6 months of the year, (2) seasonal 
Conceptual Model of Ground-Water System 23

residents were 46 percent of the population from 1960 to 

1980, and (3) the percentage of seasonal residents decreased 

from 46 to 20 percent between 1980 and 1999. Adjusted for 

seasonal residency, residential nitrate loading estimates ranged 

from 12 to 14 lb/yr (5.6 to 6.5 kg/yr) per home between 1960 

and 1999.

About 9,300 residential lots in the study area depend 

(or will depend) on on-site systems for wastewater disposal. 

Nearly 8,200 lots are in Deschutes County and the remaining 

1,100 lots are in Klamath County. About 1,800 additional 

undeveloped lots in Deschutes County may not be suitable for 

on-site wastewater systems due to extremely shallow watertable conditions (Tim Berg, Deschutes County Community 

Development Department, written commun., November 2005). 

Another 485 lots in Deschutes County will be served by a 

central sewer system scheduled to be completed in 2008 and 

are not included in these totals (B. Rich, Deschutes County 

Community Development Department, written commun., 

November 2005). As of 1999, homes with on-site systems 

had been built on 4,800 and 390 lots in Deschutes County and 

Klamath County, respectively. Assuming an average nitrogen 

loading rate of 14 lb/yr (6.5 kg/yr) per home, total residential 

loading in 1999 was 74,000 lb/yr (33,750 kg/yr). According 

to Deschutes County records, about 1,000 new homes were 

built between 2000 and 2005; this represents about 26 percent 

of the 2000 inventory of 3,900 lots in the Deschutes County 

part of the study area. Building data were not available for 

the Klamath County part of the study area for 2000–2005, 

however if the same growth rate is assumed, an additional 

180 new homes would have been built during this period. The 

estimated number of homes in the study area in 2005 was 

6,370 with nitrogen loading of 91,000 lb/yr (41,400 kg/yr). 

The remaining inventory of lots suitable for on-site systems 

in 2005 was 2,930, which could potentially add 42,000 lb/yr 

(19,000 kg/yr) of nitrogen loading from on-site systems. These 

residential loading estimates are based on the assumption that 

20 percent of residents were seasonal in 1999 and 2005. If 

all residents were full-time, the average nitrogen loading per 

home would be 16 lb/yr (7.3 kg/yr) and the total loading from 

on-site systems after all lots (9,300) were developed would be 

about 150,000 lb/yr (67,900 kg/yr), or about 164 percent of the 

2005 estimate.

As with homes, nitrogen loading for schools, motels, 

restaurants, and recreational vehicle and trailer parks was 

assumed to equal the product of the number of people using 

the facilities, average volume of effluent generated per person, 

and the average nitrogen concentration in on-site system 

effluent (also assumed to be 46 mg N/L for these facilities). 

In most instances, estimates of the number of people served 

by on-site systems in study area schools, motels, restaurants 

(data in terms of restaurant seats), and recreational vehicle 

and trailer parks were obtained from employees at the 

facilities and from the Deschutes County Department of 

Environmental Health. Data could not be obtained directly 

from some motels and recreational vehicle and trailer park 

sites; in these instances, data obtained from similar sites were 

applied. Average volumes of effluent generated in schools 

were reported to be 3.2 gal/d (12 L/d) per elementary school 

student and 5.3 gal/d (20 L/d) per middle school student (John 

Rexford, Operations Manager, Bend-La Pine Public School 

District, written commun., 2001); volumes for high school 

students were assumed to be the same as for middle school 

students. A volume of 14 gal/d (53 L/d) per restaurant seat 

(obtained by adjusting the value reported by Frimpter and 

others [1990], to estimated study area restaurant use) was used 

in estimates for restaurants. A volume of 34 gal/d (130 L/d) 

per person was used for motels (Tchobanoglous and Burton, 

1991), and 45 gal/d (170 L/d) per person (same as for homes) 

for recreational vehicle and trailer parks. In 1999, nitrogen 

loading from nonresidential sources totaled about 6,600 lb/yr 

(3,000 kg/yr). 

Annual and cumulative nitrogen loading to the groundwater system from 1960 through 2005 are shown in figure 9. 

Annual rates for 1960–99 were estimated using the method 

described above. The 2005 loading rate was based on updated 

information on the number of homes built as of 2005; the 

2000–2004 loading rates were estimated by assuming that 

the number of homes (and loading) increased linearly during 

the period and nonresidential loading remained constant at 

the 1999 rate of 6,600 lb/yr (3,000 kg/yr). Loading increased 

at a moderate rate as the population grew from 1960 to the 

mid-1970s. In the late 1970s, a surge in building (and nitrogen 

loading) occurred prior to enactment of more restrictive rules 

for on-site systems. Construction of a sewer system for the 

central business district of La Pine in 1989 caused only a 

minor reduction in overall nitrogen loading that was followed 

by rapid population growth and concurrent loading into the 

2000s. In 2005, the cumulative mass of nitrogen added to 

the ground-water system since 1960 was more than 2 million 

pounds (900,000 kg). As the number of homes contributing 

nitrogen to ground water increases, the cumulative mass of 

nitrogen added to ground water rapidly increases. One-half 

of the total mass of nitrogen added during the 45-year period 

(1960–2005) was added during the last 12 years (1993–2005) 

(fig. 9).

1960 1965 1970 1975 1980 1985 1990 1995 2000 2005

CUMULATIVE MASS, IN TONS

0

200

400

600

800

1,200

1,000

LOADING RATE, IN POUNDS PER YEAR TIMES 1,000

 0

40

20

60

80

120

100 Cumulative mass loaded

Loading rate

Figure 9. Annual and cumulative estimated nitrate 

loading from on-site wastewater systems in the La Pine, 

Oregon study area, 1960–2005.
24 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

Nitrate Fate and Transport Simulation 

Models

Computer models were used to simulate the physical and 

chemical processes affecting the fate of nitrate in the shallow 

aquifer system near La Pine. The purpose of developing the 

models was to (1) test concepts of how hydrogeologic and 

geochemical processes affect the movement of nitrate through 

the ground-water system, and (2) provide tools that could be 

used to evaluate future water-quality conditions and alternative 

management options. This section describes the process of 

constructing and calibrating simulation models at two scales: 

(1) a transect model used to estimate model parameters and 

boundary conditions along a ground-water flow path where 

there is a detailed understanding of the flow system and 

geochemical evolution of ground-water, and (2) a studyarea model covering the area where there are water-quality 

management concerns. 

This section includes: (1) discussion of the overall 

modeling approach used to construct and calibrate the transect 

and study-area models, (2) discussion of model parameters 

that must be specified for both models and how the values 

were estimated and modified during calibration, (3) discussion 

of hydrologic stresses specified for both models, (4) a detailed 

description of boundaries, discretization, and calibration of the 

transect model, and (5) a detailed description of boundaries, 

discretization, and calibration of the study-area model. 

Alternative management options were evaluated first by using 

the study-area model to simulate prescribed scenarios and then 

by coupling the study-area model with optimization methods. 

The coupled simulation and optimization model is referred to 

as the Nitrate Loading Management Model (NLMM). 

Modeling Approach

The study-area model includes 247 mi2 of southern 

Deschutes and northern Klamath Counties (fig. 10) where 

either existing or planned future development uses on-site 

wastewater systems. The study-area model was used to 

simulate the long-term effects of historical and future nitrate 

loading from on-site wastewater systems on ground-water 

quality in the study area. 

The transect model is more detailed and includes 2.4 mi2

of the shallow aquifer system near Burgess Road (fig. 10). 

The transect model was used during the calibration process to 

(1) test basic concepts of ground-water flow and nitrate fate 

and transport, and (2) refine estimates of model parameters 

through calibration in part of the study area where abundant 

geologic, hydrologic, and chemical data were available to 

constrain the model.

A regional ground-water flow model of the 4,500 mi2

upper Deschutes River basin was developed by Gannet and 

Lite, (2004) to address issues of water supply and groundwater/surface-water interaction. The model does not have the 

capability to simulate transport of nitrogen and the resolution 

(grid cell size) of the model grid is not adequate near La Pine 

to make this model an effective tool for addressing issues 

related to on-site wastewater systems. The regional model 

was used, however, to specify initial values of parameters and 

boundary conditions for the study area and transect models.

The first step in developing the transect and study-area 

models was to construct a ground-water flow model and 

simulate the velocity and direction of ground-water movement. 

This velocity distribution then was used in a solute-transport 

model to simulate the advection and dispersion of nitrate 

in the ground-water system. Steady-state flow conditions 

were assumed for the ground-water velocity distribution; 

that is, it was assumed that there were no changes in storage 

and the direction of ground-water flow and velocity did not 

significantly change during the simulation period. The steadystate flow models (transect and study area) were calibrated 

using observations (measurements) of hydraulic head, groundwater discharge rates to rivers, and chlorofluorocarbon (CFC) 

based time-of-travel estimates (Hinkle and others, 2007a). The 

accuracy of the transport model was evaluated by simulating 

nitrate loading during the 40-year period 1960–99 and 

comparing statistical distributions of observed and simulated 

ground-water nitrate concentrations. 

Ground-water flow and nitrate transport were simulated 

using the programs MODFLOW-96 and MT3DMS. The 

U.S. Geological Survey’s MODFLOW model (Harbaugh and 

McDonald, 1996) is a numerical simulation program capable 

of simulating hydraulic heads and ground-water flow within a 

three-dimensional system. The MT3DMS simulation program 

(Zheng and Wang, 1999) uses the velocity field computed by 

MODFLOW to simulate the movement of dissolved chemical 

species in ground water due to advection and dispersion. 

MT3DMS also is capable of simulating the effects of chemical 

reactions on the concentration of dissolved species.

The study-area and transect models developed for this 

study consist of a set of data input files that specify the 

boundary and initial conditions, the hydraulic and transport 

parameters, and the hydraulic and chemical sources and 

sinks specific to the study area. These input files provide 

the data required by the model programs (MODFLOW-96 

and MT3DMS) to solve a complex set of partial differential 

equations that compute ground-water head and flux as well as 

solute concentration and flux in the subsurface.
Nitrate Fate and Transport Simulation Models 25

OR19-0049_fig10

Base 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

Publication projection is 

Lambert Conformal Conic

Standard parallels 42o20' and 44o40', 

central meridian -120o30'

Datum is NAD83

R. 9 E. R. 10 E. R. 11 E.

 T.

20

 S.

 T.

21

 S.

 T.

22

 S.

 T.

19

 S.

 T.

23

 S.

 T.

24

 S.

DESCHUTES CO.

L A K E C O .

KLAMATH CO.

EXPLANATION

Boundary conditions

River cell 

Spring cell

Evapotranspiration cell

Spring

River

Fall

River

Crescent Creek

Creek

Little

Deschutes

DESCHUTES

RIVER

River

Paulina

Long Creek

5 KILOMETERS

0 5 MILES

0

Wickiup

Reservoir

97

97

31

31

Burgess Road

Drive

Finley

Butte

Road

Road

South

Century

Spring

River

Masten Road

SunriverSunriver

La 

Pine

43°

55'

50'

45'

121°40' 35' 30' 25' 121°20'

40'

43°

35'

Extent of study area model

10 20 30 40 50 60 70 80 90 100

10

20

30

40

50

60

70

80

90

100

110

120

130

140

150

160

170

180

190

200

210

220

230

240

250

260

270

Columns

Rows

Extent of 

transect model

Upper

Deschutes

Basin

Regional

Model

La Pine

 study area 

model

Long Prairie

Figure 10. Spatial relations between the upper Deschutes Basin regional ground-water model, the La Pine study area model, and 

the transect model.
26 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

Model Parameters

Several parameters must be specified for each cell in a 

model. These parameters, hydraulic conductivity, porosity, and 

dispersion coefficient, define the characteristics of the porous 

medium (geologic materials) that control the movement of 

ground water and the fate and transport of nitrate. 

Hydraulic Conductivity

Initial values of hydraulic conductivity for the five major 

lithofacies in the hydrogeologic model were estimated using 

(1) results from 24 slug tests done for this study, (2) analysis 

of 221 well-yield tests reported by drillers, and (3) literature 

values. The three fluvial hydrofacies, clay-silt, sand, and 

gravel, were assigned values of 3, 25, and 57 ft/d, respectively, 

based on mean values from the slug tests. The lacustrine 

clay-silt hydrofacies was assigned a value of 1 ft/d. The initial 

value for the basalt hydrofacies was 800 ft/d based on values 

used by Gannett and Lite (2004) in the La Pine study area. The 

initial values of vertical hydraulic conductivity were derived 

by adopting the horizontal to vertical anisotropy ratio of 1,000 

used in the regional model for the La Pine study area (Gannett 

and Lite, 2004). The resulting initial values of vertical 

hydraulic conductivity were 0.001, 0.003, 0.025, 0.057, and 

0.8 ft/d for the lacustrine clay-silt, fluvial clay-silt, fluvial 

sand, fluvial gravel, and basalt hydrofacies, respectively. 

Initial values of horizontal and vertical hydraulic conductivity 

were assigned to each cell in the model on the basis of the 

hydrofacies represented in the cell.

Porosity

Effective porosity is defined as the amount of 

interconnected pore space and fracture openings available for 

the transmission of fluids, expressed as the ratio of the volume 

of interconnected pores and openings to the volume of rock 

(Lohman, 1972). Effective porosity generally is smaller than 

the total porosity of the porous medium, reflecting that some 

pore spaces may contain immobile water with zero groundwater seepage velocity (Zheng and Wang, 1999). However, as 

discussed in some detail by Zheng and Bennett (1995), this 

so-called effective porosity cannot be readily measured in 

the field due to the complexity of the pore structure. Rather, 

it generally must be interpreted as a lumped parameter and 

estimated during calibration by comparison of simulated and 

observed solute movement or ground-water travel time. An 

initial value of 0.30 was used for the effective porosity of all 

five hydrofacies.

Dispersion

The solute transport model includes the process of 

dispersion to account for factors such as diffusion, flow path 

tortuosity, and heterogeneities in the porous media that cause 

velocities to deviate from the average seepage velocity. The 

geostatistical method used to simulate the distribution of 

hydrofacies allowed the model to represent a great deal of the 

large scale heterogeneity in hydraulic properties. The scale of 

heterogeneities incorporated in the hydrogeologic model was, 

however, limited to the scale of the model grid (200 and 500 ft 

laterally for the transect and study-area models, respectively, 

and 5 ft vertically for both models). Because much of the 

heterogeneity occurs at scales that cannot be represented by 

a grid cell, the dispersion coefficient was used to represent 

the smaller (sub-grid cell) scale heterogeneity of the system. 

Additional factors, such as temporal variations in recharge 

rates and flow directions that affect solute transport often 

are incorporated into simulations as dispersion (Goode and 

Konikow, 1990). 

Dispersion is represented by three dispersivity 

coefficients—one for longitudinal dispersivity, which 

represents dispersion along the primary flow axis, and two 

for transverse dispersivity values, which represent dispersion 

in the horizontal and vertical directions normal to the axis of 

flow. MT3DMS allows the user to specify the longitudinal 

dispersivity for each model layer and to set the transverse 

horizontal and transverse vertical dispersivity values as a 

fraction of the longitudinal dispersivity. The initial value for 

the longitudinal dispersivity was set to 50 ft to represent subgrid cell heterogeneity. Typically, the horizontal transverse 

dispersivity is much smaller than the longitudinal dispersivity 

and the vertical transverse dispersivity is smaller still 

(Zheng and Bennett, 1995). Initial values for the transverse 

dispersivity were set to 5 and 0.5 ft for horizontal and vertical, 

respectively. The final, calibrated longitudinal, transverse, 

and vertical dispersivity values were 60, 6, and 0.06 ft, 

respectively.

Model Stresses

Model stresses are the flows in and out of the groundwater system that are specified or simulated by the model. 

The stresses include recharge, discharge to springs and by 

evapotranspiration (ET), and exchange between streams 

and the ground-water system. Discharge by pumping wells 

is another stress, however it was not included in the model 

because seventy-five percent of domestic well pumping 

is returned as recharge from on-site wastewater systems. 
Nitrate Fate and Transport Simulation Models 27

Recharge is the only specified stress in the models; other 

stresses are simulated by the model based on relations between 

head in the ground-water system and stream or springs 

elevation or the rooting depth of plants that discharge ground 

water by ET.

Specified flux and head-dependent flux boundary 

conditions were used to simulate movement of water 

between the ground-water system and streams, springs, and 

the atmosphere. Recharge by infiltration of precipitation 

was represented by specifying flux to the water table using 

the recharge (RCH) package (Harbaugh and McDonald, 

1996). Streams and rivers, springs, and evapotranspiration 

were represented using head-dependent flow boundary 

options in MODFLOW-96 (RIV, DRN, and EVT packages). 

The locations of cells designated as head-dependent flux 

boundaries are shown in figure 10.

Recharge

The distribution of recharge by infiltration of 

precipitation and snowmelt was estimated for the upper 

Deschutes basin by Gannett and others (2001) using the Deep 

Percolation Model (DPM) of Bauer and Vaccaro (1987). 

The basin was modeled using a grid with cell dimensions 

of 6,000 ft per side. Gannett and Lite (2004) used estimated 

1993–95 mean annual recharge to simulate steady-state 

conditions. The 1993–95 recharge distribution was used as 

the initial estimate of recharge to the study-area model; values 

for each cell were derived by overlaying the study-area model 

grid with the DPM model grid using a GIS (fig. 5). Neither 

recharge from on-site systems nor withdrawals by domestic 

wells were included in the model because, in the context of 

the overall water budget, these components are essentially 

equal and offsetting. Reductions in water level near pumping 

wells are not represented in the study-area model using this 

approach. Individual domestic wells may influence groundwater flow within localized areas (a few cells), but they 

would not have a significant effect on the direction and rate 

of ground-water flow at the scale the model is intended to be 

used.

Rivers

Ground-water discharge to and recharge from streams 

was simulated with the river (RIV) package of MODFLOW-96 

(Harbaugh and McDonald, 1996). Data required by the RIV 

package include: mean stream stage, stream-channel elevation, 

and streambed conductance for each cell that contains a 

stream. All major streams in the study area were simulated 

using this method (fig. 10).

All stream-channel elevations were estimated using a 

10-m digital elevation model (DEM). The mean stage was 

estimated to be 4 ft above the streambed in the Deschutes 

River, 2 ft above the streambed in Paulina Creek and Little 

Deschutes River, and 1 ft above the streambed in Long Creek. 

Streambed conductance (SC) is calculated using the equation:

SC K A

b

K

v

v

= ,

where

is the vertical hydraulic conductivity of the 

streambed (LT ),

is the area of the streambed within t

1

A he cell (L ), and

is the thickness of the streambed (L).

2

b

Conductance values were variable depending on the 

hydrofacies of the cell and the area of the stream within 

the cell. The initial conductance values were calculated by 

assuming that the vertical hydraulic conductivity of streambed 

sediments was 1/100th of the horizontal hydraulic conductivity 

of the cell. Streambed thickness was set to 1 ft throughout the 

model.

Springs

Ground-water discharge to springs was simulated with 

the drain (DRN) package of MODFLOW-96 (Harbaugh 

and McDonald, 1996). Data required by the DRN package 

include the spring elevation and conductance for each cell 

that contains a spring. This package was used to represent 

discharge from springs at the contact between the basinfill sediments and volcanic rocks on the west side of the 

Deschutes River between the confluence with Spring River 

and the northern boundary of the model (fig. 10). All spring 

elevations were estimated using a 10-m digital elevation 

model (DEM) data. The conductance parameter has the same 

dimensions (L2

T-1) as the streambed conductance; however, 

the parameter cannot be easily related to measurable physical 

characteristics of the springs. Initial values of conductance 

ranged from 5,000 to 10,000 ft2/d and were adjusted during 

calibration. 

Evapotranspiration

Discharge from the water table by evapotranspiration 

(ET) was simulated using the EVT package of MODFLOW-96 

(Harbaugh and McDonald, 1996). ET can occur where plants 

with roots extending to the water table (phreatophytes) use 

ground water for transpiration or where the water table is 

close enough to land surface to allow direct evaporation from 

bare soil. The EVT package uses the assumption of a linear 

relation between ET discharge and depth to the water table. 

Data required to define this relation are the potential ET rate 

(LT-1), depth below land surface at which ET ceases (extinction 
28 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

depth), and land surface elevation (L). The potential ET rate 

was specified as 1.0 ft/d, the extinction depth was 5 ft, and the 

land surface elevation was specified from DEM data. Only 

cells where the water table was estimated to be within 10 ft of 

land surface were activated in the EVT package (fig. 10). The 

nitrate concentration in ground water used by plants in the ET 

process was specified as zero because the (1) part of modelsimulated ET that is due to plant transpiration (as opposed to 

evaporation from bare soil) could not be determined, and (2) 

nutrient uptake by phreatophytes is species dependent and not 

well constrained.

Transect Model

The transect model was developed for a part of the 

ground-water system near Burgess Road (fig. 10). This area 

was selected for detailed scale modeling because it was 

the site of intensive data collection and study of changes 

in ground-water age and chemistry along a ground-water 

flow path transect (Hinkle and others, 2007a). The 3.5-mi 

long transect model comprises 17 wells at 6 sites (2–4 wells 

per site) aligned with the predominant direction of groundwater flow (fig. 11). In addition to ground-water age and 

chemistry data, hydraulic head, lithologic data, and hydraulic 

conductivity data were available from the transect wells to 

construct and calibrate the transect model. 

Boundaries and Discretization

The boundaries of the transect model (fig. 11) were 

selected to coincide with hydrogeologic boundaries of the 

local flow system. The western boundary of the model is a 

no-flow boundary defined by the ground-water divide that 

coincides with the surface-water drainage divide between the 

Deschutes and Little Deschutes Rivers. The north and south 

boundaries of the model are streamline boundaries that parallel 

the principal direction of ground-water flow from the drainage 

divide to the Little Deschutes River. The Little Deschutes 

River, a discharge area, was simulated as a head-dependent 

flux boundary. The initial conceptual model of the system 

included a deep regional component of upward ground-water 

flow across the lower model boundary. The concept for the 

lower boundary was evaluated by applying a specified flux 

across the lower model boundary based on data from the 

regional ground-water flow model by Gannett and Lite (2004). 

The upward flux estimates from the regional model did not 

allow a good simulated match to observed conditions (head 

and flux) when hydraulic conductivity and recharge values 

were within reasonable ranges based on independent data. 

Little or no upward flow from a deep regional source is needed 

to simulate the observed distribution of head and estimated 

ground-water flux to the Little Deschutes River. Upward 

vertical gradients were measured in well pairs in the area and 

geochemical evidence indicates that ground-water discharge 

to areas near the rivers may have a deep regional component 

(Hinkle and others, 2007a). Seepage estimates and simulations 

using the transect model, however, indicate that the magnitude 

of regional discharge must be small compared to local 

recharge to the shallow alluvial aquifer. The lower boundary of 

the transect model was specified as a no-flow boundary.

The numerical methods used by MODFLOW and 

MT3DMS to solve the equations for ground-water flow and 

solute transport require that the ground-water system be 

divided into discrete rectangular volumes, called cells. The 

rows, columns, and layers of cells that make up the threedimensional array are collectively referred to as the model 

grid. The center of each cell defines a point at which the 

hydraulic head and solute concentration is simulated. The 

simulated head and concentration values are taken to represent 

the average values within the cell.

The longest dimension of the model grid was aligned 

to the principal direction of ground-water flow from west to 

east. The grid consists of 18 rows, 92 columns and 15 layers 

of cells, with each cell having dimensions of 200 ft in the row 

and column directions and 5 ft in the vertical direction. The 

three-dimensional transect model was 3.5 mi in the east-west 

(column) direction, 0.7 mi wide in the north-south (row) 

direction (normal to the principal direction of flow), and 75 ft 

thick in the vertical (layer) direction (fig. 11). 

A three-dimensional hydrofacies model of the fluvial 

and lacustrine sediments within the transect flow model 

was constructed using the methods previously described for 

constructing the hydrofacies model for the study-area model. 

A sectional view through the model along row 8 is shown 

in figure 11C and illustrates the heterogeneity of the fluvial 

sediments. Transition probability geostatistical methods were 

used with the same model parameters (mean lengths, transition 

probabilities, and proportions) as were used for the study-area 

hydrofacies model, however the transect hydrofacies model 

was discretized into cells with the same dimensions as the 

transect simulation model.

Calibration and Sensitivity

The MODFLOW model was first calibrated to steadystate flow conditions using head data from the June 2000 

synoptic well measurement. Data were available from 26 

wells within the transect model area (fig. 11); 17 wells were 

piezometers installed as part of the Burgess Road transect 

study and 9 were existing domestic water-supply wells. 

Observations (either measurements or estimates) 

of ground-water discharge to streams were available to 

constrain model calibration; however, flux between the Little 

Deschutes River and the shallow ground-water system was 

difficult to characterize because of uncertainty in the seepage 

measurements. The apparent gains and losses for the 13-mi 
Nitrate Fate and Transport Simulation Models 29

Figure 11. Plan and section views of the transect model showing simulated water-levels, ground-water travel time, and particle 

paths in the La Pine, Oregon, study area. (Location of the transect model is shown in figure 10.)

OR19-0049_fig 11

1

2 3

4 5

6

1

2 3

4

5

6

1

2 3

4

5

6

A. PLAN VIEW OF TRANSECT MODEL GRID

B. SIMULATED TRAVEL TIME

C. SIMULATED PARTICLE PATHS

4230

4225

4220

4215

4210

4205

4200 4195

4232

4228

4232

4224 4226

4226

4226

4223

4224 4212

4196 4192

4198

4210

4226

1.1 – 11 yrs

1.2 – 27

2.1 – 16 yrs

2.2 – 31

2.3 – 52

2.4 – 53

3.1 – 20 yrs

3.2 – 25

3.3 – 43

4.1 – 24 yrs

4.2 – 25

4.3 – 38 5.1 – 33 yrs

5.2 – 48

5.3 – >60 6.1 – <34 yrs

6.2 – >47 Little Deschutes River Row 8

75 FEET 600 FEET 75 FEET

18,400 FEET

Unsaturated

Lacustrine clay-silt

Volcanics: basalt, andesite

Fluvial gravel

Fluvial sand

Fluvial clay-silt

2

4230

Cell

Cell simulated as stream

Contour of simulated water-table 

 elevation

Well location

Simulated travel time, in years

Transect well location

Measured water-table elevation, 

 in feet, June 2000

Transect well No.

EXPLANATION

A. PLAN VIEW OF TRANSECT MODEL GRID B. SIMULATED TRAVEL TIME

C. SIMULATED PARTICLE PATHS

100-110

110-120

120-140

140-160

160-180

180-200

200-230

230-260

260-310

310-380

700-1,000

380-470

470-600

600-700

0-2

2-5

5-10

10-15

15-20

20-25

25-30

30-35

35-40

40-45

45-50

50-55

55-60

60-65

65-70

70-80

80-90

90-100

Transect well screen location

Simulated particle path

2 Transect well No.

Hydrofacies

Cell simulated as stream

Transect well screen location

2 Transect well

Transect well and screen number with 

1.1 – 11 yrs Chloroflourocarbon-derived ground-water age

4228
30 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

reach downstream of the gaging station (14063000) ranged 

from a 2.6 ft3/s gain in October 1995 to a 1.6 ft3/s loss in 

October 2000 (table 4). These 1995 and 2000 estimates were 

less than the estimated error of measurement (± 3.9 ft3/s and 

±6 ft3/s, respectively). A third seepage measurement in March 

2000 indicated a 12 ft3/s loss in this reach, but this loss was 

based on measurements at relatively high flow conditions 

(about 300 ft3/s) with a larger associated error. For calibration 

purposes, October low-flow data were interpreted to indicate 

that net seepage to or from this reach (1) is small relative to 

flow, and (2) probably varies in direction as well as magnitude 

depending on annual or seasonal head relations between 

shallow ground water and the stream. A key component of 

the conceptual model of the shallow flow system is that, on 

annual or longer time scales, there is net discharge to the 

Little Deschutes River. On this basis, a range of discharge 

to the stream for the transect model area was computed 

using the October 1995 seepage (2.6 ft3

/s) as the minimum 

and the sum of seepage and estimated error of measurement 

(6.5 ft3/s) as the probable maximum for the range. Dividing 

seepage and estimated error of measurement by the length of 

the reach (13 mi) yielded a range of 0.20 to 0.50 ft3/s/mi and 

the resulting range of expected mean annual ground-water 

discharge to the 1.3 mi reach of the Little Deschutes in the 

model area was 0.3 to 0.7 ft3/s (rounded). 

During the calibration process, parameters such as 

horizontal and vertical hydraulic conductivity, recharge, and 

streambed conductance were adjusted, using a trial-anderror procedure, to obtain the best fit between observed and 

simulated head and stream flux values. The transect model 

also was used to evaluate estimates of upward vertical flux 

across the lower boundary from the regional flow model. 

Horizontal hydraulic conductivity (Kh) and vertical 

hydraulic conductivity (Kv) values were assigned to each 

cell according to the hydrofacies represented in the cell. The 

probable range of Kh based on well-yield data and slug tests 

and the calibrated Kh value that yielded the best model fit 

are listed in table 5. The probable range of vertical hydraulic 

conductivity for each hydrofacies was computed assuming that 

the probable range of the horizontal to vertical anisotropy ratio 

was 10 to 10,000 for the hydrofacies in the model area. The 

vertical anisotrophy ratio in the final model was equal to 100 

for all hydrofacies. 

Streambed conductance was estimated on the basis 

of the assumption that the vertical hydraulic conductivity 

in streambed sediments was 1/100th of the horizontal 

conductivity of the hydrofacies in the model cell. A series of 

simulations were made to test model sensitivity to streambed 

conductance values. Conductance values were modified 

by adjusting the ratio of streambed Kv

 to aquifer Kh, over a 

range of 0.001 to 0.1, but the transect model was relatively 

insensitive to this parameter.

Recharge estimates for the transect model area range 

from 1.5 to 2.1 in/yr. The lower end of the range was estimated 

by using a water balance model for the upper Deschutes Basin 

(Gannett and others, 2001); the upper end of the range was 

computed by using CFC-derived age data from the Burgess 

Road transect study (Hinkle and others, 2007a). The best 

model fit was obtained with a slightly higher recharge rate 

of 2.3 in/yr. 

Hydraulic heads and discharge to the stream simulated 

by the transect model were most sensitive to recharge and 

horizontal hydraulic conductivity of the fluvial clay-silt 

and lacustrine clay-silt hydrofacies. The fluvial sediments 

are relatively thinner than the lower permeability lacustrine 

sediments near the river and this tends to reduce the composite 

transmissivity of the aquifer system. The location, thickness, 

and extent of clay-silt hydrofacies within the fluvial sediments 

are relatively unimportant factors controlling ground-water 

flow because the contrast in hydraulic conductivity between 

the gravel, sand, and clay-silt hydrofacies is small compared 

to the contrast between the fluvial and lacustrine hydrofacies. 

Fitting simulated to observed heads was accomplished 

by adjusting values of hydraulic conductivity for each 

hydrofacies, until a best model fit was obtained by trial and 

error. 

Table 5. Values of horizontal and vertical hydraulic conductivity for 

hydrofacies based on field data and model calibration in the La Pine, Oregon, 

study area.

[Abbreviations: Kh, horizontal hydraulic conductivity Kv, vertical hydraulic conductivity]

Hydrofacies

Horizontal hydraulic 

conductivity, Kh

(feet per day)

Vertical hydraulic 

conductivity, Kv

(feet per day)

Probable 

range

Calibrated 

value

Probable 

range

Calibrated 

value

Fluvial clay-silt 1–20 5 0.001–2 0.05

Fluvial sand 10–100 40 0.01–10 .40

Fluvial gravel 75–150 75 0.075–15 .75

Lacustrine clay-silt 0.1–10 10 0.0001–1 .10

Basalt 100–2,500 150 0.01–250 1.5
Nitrate Fate and Transport Simulation Models 31

The calibrated parameter values (table 5) and boundary 

conditions produced a good model fit to head and flux 

observations. Contours of simulated water-table elevation 

compare well with observed water-table elevations at 15 wells 

(fig. 11A) where measured heads represent the water table. 

The root mean square error (RMSE) of the head residuals 

at the 26 well locations was 2.7 ft. Model fit was slightly 

better (RMSE = 2.3 ft) at the 17 transect wells, where head 

measurements were more accurate because (1) measuring 

point elevations at the wells were surveyed, whereas elevations 

at other wells were estimated from topographic maps, and 

(2) short (2-ft) screens in the transect wells provided better 

observations of head within the 5-ft thick model cells, whereas 

screens in other wells were longer and spanned multiple 

model layers. The simulated heads have a slight bias that 

results in a steeper slope of the water table toward the Little 

Deschutes River. The bias probably could be reduced with 

further calibration if the configuration of the lacustrine-fluvial 

sediment boundary were modified; however, there was no 

basis (using existing geologic data) to modify the hydrofacies 

model and the bias was small and did not significantly affect 

average simulated flow velocities.

Simulated Travel Time and Comparison with 

Ground-Water Age

Estimates of ground-water age were used to further 

constrain the calibration of the transect model. Ground-water 

age estimates at 17 wells in 6 locations along the Burgess 

Road transect (fig. 11B, C) were derived from analysis of 

CFCs (Hinkle and others, 2007a). A particle tracking program, 

MODPATH, (Pollock, 1994) was used to determine the 

simulated travel time of ground-water from recharge locations 

at the water table to the well screens where CFC samples were 

collected from the 17 transect wells. The simulated travel 

times were compared with the CFC-derived age estimates to 

evaluate the accuracy of the simulated ground-water velocity 

distribution. 

MODPATH uses the cell-to-cell fluxes simulated by 

MODFLOW to compute velocity field and path lines of 

individual “particles” of water that may be started at any 

location and tracked either forward to their discharge location 

or backward to their recharge location. The only additional 

parameter required by MODPATH (that is not required by 

MODFLOW) is effective porosity. Effective porosity was 

assumed to be 0.3 for each hydrofacies within the model. For 

this analysis, 1,000 particles were placed in each of the 17 

model cells that contained the screen of a transect well. The 

particles were tracked backward to their recharge locations 

at the water table (fig. 11C). Because of limitations of the 

tracking program, particles could not be placed at the exact 

location of the well screens within each cell, but were placed 

at the center of the cell and distributed along a vertical line 5-ft 

long (the thickness of the cell). 

Simulated ground-water travel time distribution is shown 

in figure 11B for row 8 of the transect model. The simulated 

distribution shows a pattern similar to the distribution of 

CFC-derived ages from the 17 transect wells. Generally, 

ground-water travel times are less than 10 years within the 

upper 5–15 ft of the saturated zone and travel times within 

30–50 ft of the water table are less than 50 years. Simulated 

ground-water travel times are notably greater than CFC ages 

at well 4.3, the deepest at site 4. Travel times are difficult 

to predict at this site because the screen for well 4.3 is near 

the upper surface of the fine-grained, lacustrine clay-silt 

hydrofacies. The model predicts that “older” water will occur 

at this location because slow-moving water from the lacustrine 

sediments re-enters the fluvial sediments as it follows flow 

paths toward the discharge areas along the river. At the 

discharge end of the flow path, near the Little Deschutes River, 

older water occurs at shallower depths as it moves upward to 

discharge to the river. The relation between ground-water age 

and depth depicted in the section (fig. 11B) is analogous to the 

relation between dissolved oxygen concentration and depth 

described by Hinkle and others (2007a).

Study-Area Model

Initial estimates for many parameters in the studyarea model were derived from calibration and testing of the 

transect model. The study-area model was first calibrated 

for simulation of ground-water flow using observations of 

head and discharge to streams. As with the transect model 

calibration, simulated heads were compared with waterlevel data from the synoptic measurement of 228 wells in 

June 2000 (fig. 5). Simulated ground-water discharge to 

streams was compared with estimates of discharge from 

gain-loss measurements (fig. 7, table 4). The measured head 

and discharge data were assumed to represent the long-term 

average (steady-state) conditions in the ground-water system. 

Once the study-area flow model was calibrated, the simulated 

ground-water flow directions and velocity were used to 

simulate the transport and fate of nitrate. Transport of nitrate 

from on-site wastewater systems was simulated for 1960–99 

using the estimated annual nitrate loading rates. The transport 

model was evaluated by comparing the statistical distribution 

of simulated nitrate concentration with statistical distributions 

of measured nitrate concentrations from analyses available 

mostly from domestic wells. 

The transport model calibration period was divided into 

40 1-year stress periods. During each 1-year period, the nitrate 

loading rate in each model cell was specified. The specified 

rate was dependent on the number and type (domestic, 

commercial, etc.) of on-site wastewater systems present. 
32 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

Boundaries and Discretization

The model grid for the La Pine study-area model includes 

an area of 247 mi2 (fig. 10). The grid is elongated along a 

northeast-southwest trend that coincides with the orientation 

of the structural basin, major drainages, and the depositional 

strike and dip of the basin-fill sediments. The axes of the grid 

also are orthogonal to the principal directions of ground-water 

flow. There are 276 rows, 100 columns, and 24 layers in the 

model grid. The lateral dimensions of each cell are 500 ft in 

both the x and y directions. The upper 23 layers of the model 

were each 5 ft thick and the bottom layer was 100 ft thick. The 

coordinates of the lower left corner of the grid are 4,609,651 

ft east and 696,295 ft north (Oregon State Plane north) and the 

grid is rotated 15 degrees clockwise.

The fluvial sand and gravel aquifers that supply most 

drinking water to wells in the area lie within the upper 100 ft 

of sediments. Between depths of 10 to 100 ft, depending on 

location in the study area, there is a transition from the coarser 

fluvial sediments to predominately fine-grained lacustrine 

sediments (fig. 2, pl. 1). As previously described, the fluvial 

deposits are heterogeneous and include fine-grained facies of 

clay, silt, and silty-sand. Similarly, coarse facies occur within 

the lacustrine deposits, although they tend to be isolated 

lenses. The simulation model was constructed to include the 

upper 120 ft of sediments so that the model would include the 

entire thickness of fluvial sediments. Each of the 24 model 

layers had a thickness of 5 ft throughout the grid. The grid 

for the study-area flow and transport model was defined to 

exactly match the dimensions of the grid used to prepare the 

three-dimensional hydrogeologic model described earlier in 

this report. 

Long-term well hydrographs (fig. 6) indicate that groundwater flow in the shallow part of the system is in a state of 

dynamic equilibrium. Seasonal and interannual variations 

in water level and ground-water storage are due to climatic 

variation, but no long term trend is evident in mean water 

levels. The June 2000 synoptic water level measurements 

(fig. 5) were assumed to represent the mean annual head 

distribution for the system and were used as the head 

observations used as calibration targets for the steady-state 

study-area and transect flow models. 

The lateral and bottom boundaries of the model were 

selected to include the most densely developed part of the 

study area and the shallow alluvial deposits that form the 

primary aquifers. The lateral and bottom boundaries of the 

model do not coincide with specific hydrologic or geologic 

boundaries such as geologic contacts or ground-water flow 

divides. The potential magnitude of flow across the lateral 

and lower model boundaries was evaluated by model testing 

and sensitivity analysis with both the transect and study area 

models. The upper boundary of the model is the water table, 

which is the top of the saturated part of the ground-water flow 

system. 

The study area transport model used steady-state flow 

velocities to simulate nitrate fate and transport. Nitrate loading 

from on-site wastewater systems was simulated for 1960–99 

in forty 1-year periods. The loading for each period was 

computed using building records and other information as 

previously described (fig. 9). Nitrate loading was specified 

in the transport model using the mass-loading option in the 

source/sink mixing package of MT3DMS (Zheng and Wang, 

1999). The mass loading rate to each cell, in kilograms per 

day, was specified at the water table for each year of the 

simulation. Initial nitrate concentrations were assumed to be 

zero at the beginning of the simulation in 1960.

Based on the conceptual model for denitrification (Hinkle 

and others, 2007a), the transition zone from oxic to suboxic 

conditions was represented as a sharp boundary that does not 

vary significantly with time. Although the redox boundary 

would be expected to migrate as solid-phase electron donors 

are consumed by dissolved oxygen and nitrate reduction, the 

effectiveness of this approach for characterizing current zones 

of nitrate stability and instability is supported by analysis of 

data collected by ODEQ (Hinkle and others, 2007a). Hinkle 

presents field data from the area showing that 95 percent of 

samples with nitrate concentrations greater than 1.0 mg N/L 

taken from wells with known construction were screened at 

least partly within the mapped (fig. 8) oxic zone.

The position of the boundary within the models was 

specified by subtracting the thickness of the oxic part of the 

aquifer (fig. 8) from the water-table elevation (fig. 5) to obtain 

the elevation of the oxic-suboxic boundary. The boundary 

was assumed to be in each vertical stack of cells at the top of 

the first fully suboxic cell. The boundary was implemented in 

MT3DMS by using the ICBUND array to specify a constant 

concentration of zero in suboxic cells (Zheng and Wang, 

1999).

Calibration and Sensitivity

Calibration of Ground-Water Flow

Ground-water fluxes across the lateral and lower 

boundaries of the study-area model initially were specified 

on the basis of the regional ground-water flow model 

developed by Gannett and Lite (2004). The uppermost layer 

of the regional model was 100 ft thick and represented both 

alluvial deposits and basalts. The regional model simulates 

large subsurface flow through the permeable basalts that 

underlie the northern part of the study area (fig. 2). Most 

of the flow remains within the basalts of the regional flow 

system and discharges to the Deschutes River and tributaries 

north of the study area (Gannett and others, 2001). Some of 

the flow through the basalts, however, discharges to springs 

that contribute to the Deschutes, Fall, and Spring Rivers in 

the northwestern part of the study-area model. These springs 
Nitrate Fate and Transport Simulation Models 33

generally are located where high permeability basalts abut 

lower permeability sediments, forcing ground water to 

discharge at land surface (Gannett and others, 2001). 

Where alluvial deposits occur at the lower and lateral 

boundaries of the study area model, the flux simulated 

between the regional and study-area models was generally 

small. This is consistent with results of testing the transect 

model that indicated upward leakage across the lower 

boundary was not required to match the observed head 

distribution and ground-water discharge to streams when using 

reasonable values for recharge and hydraulic conductivity. 

This result indicates that upward leakage is small compared to 

local recharge from infiltration of precipitation.

When the simulated flux across the lower and lateral 

model boundary was included in the simulation, the study-area 

model had numerical stability problems in the northern part 

of the model area, which is underlain by basalt. The stability 

problems were characterized by large fluctuations in simulated 

heads that were probably caused by (1) a large flux through a 

relatively thin section (120 ft) of the basalts, and (2) adjacent 

model cells with large contrasts in hydraulic conductivity 

(high hydraulic conductivity basalts in immediate contact with 

low conductivity alluvial deposits). 

To test the sensitivity of model results to specified flux at 

the lateral and lower model boundaries, a simulation was made 

specifying zero flux at these boundaries. Simulated heads 

and flux to streams were determined to be insensitive to the 

regional model boundary flux where basalts were not present. 

Few data points were available to constrain the heads in the 

basalt uplands and, therefore, even though simulated heads 

were affected in the basalts, it is difficult to assess how well 

or poorly they match actual heads. The simulated discharge 

to springs along the Fall, Spring, and Deschutes Rivers in 

the northern part of the model area also were affected, but 

the model showed much better stability. Nearly all potential 

concerns related to nitrate from on-site wastewater systems 

are focused on parts of the study area underlain by alluvial 

deposits. Simulation of flux through the basalts and of 

discharge from large springs at the basalt-alluvium contact 

was of less importance and did not inhibit the ability of 

the model to simulate conditions in the alluvial deposits; 

therefore, to attain better model stability in the part of the 

model representing basalts, the zero-flux boundary condition 

was used for the lower and lateral boundaries of the study-area 

model.

The contours defining the simulated steady-state 

potentiometric head surface retain the general features of the 

surface defined by contours of measured heads from the 2000 

synoptic measurement (fig. 12). Both simulated and observed 

flow directions (inferred from contours) support other data 

(vertical head gradient and gain-loss measurements) that 

indicate, on an annual basis, ground water discharges to nearly 

all stream reaches in the study area. 

The model simulated the general direction of groundwater flow (northeast) accurately in the area between the 

Deschutes and Little Deschutes Rivers in southeast T20S/ 

R10E and northeast T21S/R10E; however, the model did not 

simulate as steep a horizontal hydraulic gradient as indicated 

by the contours of observed heads shown in figure 12. The 

location of this area between the two rivers and at the focal 

point of ground-water flow for a large part of the study 

area, make it difficult to characterize the mean annual head 

distribution. Seasonal differences in the stage of the Deschutes 

and Little Deschutes Rivers can cause relatively large seasonal 

variations (5–10 ft) in ground-water levels. The period of 

high river stage in spring coincides with the peak recharge 

period for the shallow aquifer. These conditions lead to larger 

horizontal head gradients, particularly in this area, when the 

water levels used to construct the contours in figure 12 were 

measured (June 2000). The simulated recharge and river 

stage conditions represent mean annual conditions, as do the 

simulated heads, whereas the contours shown in figure 12 may 

represent a seasonal extreme. The transient influence of the 

rivers makes it difficult to accurately simulate ground-water 

flow velocity and direction in this area; however, the model 

simulates the mean annual direction of flow and velocity 

reasonably well. On-site wastewater systems are not used in 

this area so it is not likely that the model would be used to 

predict concentration in this area. Long Creek was represented 

in the model using head-dependent flux boundaries for ET and 

streams to account for discharge in the topographically low 

area; the effects of ground-water discharge are reflected in the 

V-shaped simulated and observed contours (fig. 12).

The mean absolute error for the 170 head observations 

was 7.0 ft, or about 5 percent of the range of observed heads. 

The difference between simulated and observed head ranged 

from -20 to 45 ft (fig. 13). The largest differences generally 

were at wells with long open intervals, in which the measured 

head represents a composite head for a large thickness of the 

system. The simulated heads used for comparison in such 

cases were for the model layer closest to the center of the open 

interval of the well. The differences between simulated and 

observed head generally were less than 10 ft in the central part 

of the study area, between the Deschutes and Little Deschutes 

Rivers. Errors in the simulated heads will result in errors in the 

ground-water velocity field that in turn will affect simulation 

of nitrate transport. These errors would affect simulated 

concentrations at small scales, but would not affect the average 

simulated nitrate concentration over large areas. 

Recharge to the water table, which was 98 percent of total 

recharge, was the only specified component of the simulated 

water budget. The mean annual recharge rate of 58 ft3

/s was 

equivalent to an average of 3.2 in/yr within the study area 

model. The other components of recharge and discharge were 

simulated using head-dependent boundary conditions, as 

previously described. The simulated recharge from streams 
34 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

OR19-0049_fig12

Base 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

Publication projection is 

Lambert Conformal Conic

Standard parallels 42o20' and 44o40', 

central meridian -120o30'

Datum is NAD83

R 9 E R 10 E R 11 E

T 21 S

T 22 S

T 23 S

T 24 S

T 20 S

T 19 S

DESCHUTES

L A K E

KLAMATH

Spring

River

Fall

River

Crescent Creek

Creek

Little

Deschutes

DESCHUTES

RIVER

River

Paulina

Long Creek

5 KILOMETERS

0 5 MILES

0

Wickiup

Reservoir

97

97

31

31

Burgess Road

Drive

Finley

Butte

Road

Road

South

Century

Spring

River

Masten Road

SunriverSunriver

La 

Pine

43°55'

50'

45'

121°40' 35' 30' 25' 121°20'

40'

43°35'

Extent of study area model

EXPLANATION

Water-level contours—In 

feet above National 

Geodetic Vertical Datum 

of 1929. Contour interval 

is 20 feet.

Simulated

Measured (June 2000)

4200

4200

Long Prairie

4220

4210

4170

4190

4180

4240

4160

4210

4220

4170

4230

4190

4180

4160

4260

4200

4250

4200

4200

4190

4200

4180

4260

4270

4280

4250

4240

4220

4230

4210

4290

4160

4300

4310

4320

4330

4340

4150

4170

4240

4170

4220

4230

4170

4290

4280

4270

4260

4250

4240

4230

4220

4210

4200

4190

4180

4170

4160

Figure 12. Contours of simulated and observed heads (June 2000) for the La Pine, Oregon, study area.
Nitrate Fate and Transport Simulation Models 35

was small (1.3 ft3/s) and was from a few isolated stream 

reaches where simulated heads were less than the specified 

stage of the stream. Because the simulated system is assumed 

to be at steady state (no long-term change in storage), total 

discharge is equal to recharge. Ground-water discharge to 

streams accounts for 67 percent (39.5 ft3/s) of total discharge. 

Simulated discharge to the Little Deschutes River is 15 ft3/s 

which is within the range of 7 to 20 ft3/s expected on the 

basis of measured discharge (table 4). Simulated groundwater discharge to the Deschutes River was 12 ft3

/s and was 

consistent with measurement data for reaches upstream of 

river mile 199.7. Downstream of river mile 199.7 on the 

Deschutes River and on the Fall River, most ground-water 

discharge emanates from the large springs where basaltic rocks 

are in contact with the lower permeability alluvial sediments. 

Recharge by subsurface inflow to the basaltic rocks that feed 

these springs was not simulated in the study-area model. The 

simulated discharge to ET of 16 ft3/s fell within the estimated 

range of 10–20 ft3/s and was distributed within the floodplain, 

where shallow water-table conditions persist through the dry 

months (fig. 10). 

Comparison of Simulated and Measured Nitrate 

Concentrations

The study-area transport model simulates nitrate 

concentrations in ground water and in ground-water discharge 

to the near-stream environment. Simulated concentrations 

are averages for the 500-ft wide by 500-ft long by 5-ft thick 

model cells. With cells covering nearly 6 acres and minimum 

lot sizes of 0.5 acre, each cell can contain as many as about 

10 homes. Nitrate data collected for this study from closely 

spaced sampling locations near the Burgess Road transect 

model (Hinkle and others, 2007a) indicate that even in 

mature, high-density residential areas, nitrate plumes have 

not coalesced to a great degree, and concentrations are highly 

variable at the scale of an individual model cell. Because 

of the high variability of nitrate concentrations in a cell, 

concentrations at individual wells cannot be simulated with the 

study-area model. The inability to delineate the edges of, and 

concentrations within, individual solute plumes is a limitation 

of transport models at the watershed scale. This limitation 

does not affect this study because the information from 

the model is intended to help understand and predict water 

quality conditions at scales larger than individual plumes 

or wells; however, it does limit the degree that measured 

nitrate concentration data from wells can be used for direct 

comparison with simulated concentrations.

To assess of the ability of the study-area model to 

represent the primary processes that affect nitrate movement 

in the ground-water system, the statistical distribution 

of simulated nitrate concentrations was compared with 

distributions for two sets of measured nitrate concentrations 

from wells. The first measured dataset was from a synoptic 

sampling of 192 wells in June 2000 by ODEQ (Hinkle 

and others, 2007a). Only data from the 109 wells where 

ground water was oxic (dissolved oxygen concentration was 

greater than 0.5 mg/L) were used in the comparison because 

denitrification has been shown to be an important process 

where ground water is suboxic (Hinkle and others, 2007a). 

The second observed dataset was collected under a program 

administered by Oregon Department of Human Services 

Health Division (DHS), which requires that water from 

domestic wells is tested whenever a property is sold. Nitrate 

analyses from 1,572 such tests were available for homes in the 

La Pine area (Rob Keller, ODEQ, written commun., August 

2006). The DHS data were collected from 1989 to 2004. 

Dissolved oxygen concentrations are not analyzed as part of 

the DHS program so it was not possible to discriminate wells 

that pump from the suboxic part of the system. 

The simulated nitrate concentrations used for comparison 

were from the end of the simulation period (1999) and were 

taken from 1,398 cells randomly selected from locations 

where active on-site wastewater systems existed. Only cells 

that contained oxic ground water were selected (because 

suboxic cells would have simulated concentrations of zero 

by default) and more than one cell could be selected from 

more than one layer in the same row/column. The statistical 

distributions of measured nitrate concentrations and the 

simulated concentrations are similar (fig. 14). The maximum 

simulated nitrate concentration was 29 mg N/L, with a mean 

of 2.0 mg N/L, and a median of 0.8 mg N/L, and 10 percent 

Figure 13. Simulated head residuals and observed heads 

(June 2000) from the La Pine, Oregon, study area.

OR19-0049_fig13

OBSERVED HEAD, IN FEET ABOVE NORTH AMERICAN DATUM 1983

4120 4140 4160 4180 4200 4220 4240 4260 4280 4300

RESIDUAL OF SIMULATED HEAD , IN FEET

-30

-20

-10

 0

10

20

30

40

50
36 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

of concentrations greater than 6 mg N/L. The maximum of 

the ODEQ June 2000 synoptic nitrate concentration data 

was 26 mg N/L, with a mean of 1.6 mg N/L, and a median 

of 0.3 mg N/L, and 10 percent of concentrations greater 

than 4 mg N/L. The maximum of the DHS real estate 

nitrate concentration data was 22 mg N/L, with a mean of 

1.6 mg N/L, and a median of 0.5 mg N/L, and 10 percent 

of concentrations greater than 4.5 mg N/L. The primary 

difference in the three nitrate concentration distributions was 

the slightly greater proportion of high values in the simulated 

concentration distribution. This difference is likely due to 

simulated values being sampled from the entire thickness 

of the oxic part of the system, including cells near the water 

table where nitrate loading occurs and concentrations are 

greatest. Samples from the measured datasets were collected 

from wells where the screened intervals typically were below 

the water table and would be less likely to include water 

with high nitrate concentrations. Good agreement between 

the summary statistics of the measured and simulated nitrate 

concentrations (mean, median, 90th percentile, and maximum) 

indicates that the simulated mass of nitrate in the ground-water 

system at the end of the 1960–99 period, is similar to the mass 

indicated by available sample data. This agreement increases 

confidence that the primary processes affecting the fate and 

transport of nitrate in the ground-water system are represented 

in the simulation model. Even though the model does not 

simulate concentrations at individual wells, it is a useful tool 

for assessing the effects of on-site systems on average groundwater nitrate concentrations at the scale required for evaluation 

of management alternatives for protecting ground-water 

quality.

The spatial distribution of simulated nitrate concentration 

at the water table in 1999 is shown in figure 15 and closely 

mirrors the locations of on-site wastewater systems (fig. 1). 

The effect of ground-water movement on nitrate concentration 

is evident where areas of high concentration are elongated 

parallel to the primary directions of ground-water flow, 

such as immediately south of Burgess Road. The effect of 

denitrification on the simulated distribution is evident where 

concentrations sharply decrease along easterly ground-water 

flow paths that terminate at the Little Deschutes River, such 

as in central T21S R10E (fig. 15). The sharp concentration 

gradient is coincident with an area where the oxic part of the 

system decreases in thickness (compare fig. 8). This decrease, 

along with the downward component of advective transport, 

forces a large fraction of the nitrate in the system to be 

transported into the suboxic zone and lost to denitrification.

At the end of the simulation period (1999), the rate 

of nitrate (as N) loading to the ground-water system 

was 82,000 lb/yr (37,000 kg/yr). The simulated rate of 

denitrification in the suboxic part of the system was 31,000 lb/

yr (13,900 kg/yr) and the simulated discharge of nitrate to the 

near-stream environment was 8,000 lb/yr (3,650 kg/yr). The 

remaining 43,000 lb/yr (19,400 kg/yr) was added to storage 

in the shallow ground-water system. Nitrate added to storage 

increased the mean concentration in the ground-water system 

from essentially zero in 1960 to a mean of 2 mg N/L in 1999.

Management Scenario Simulations

Simulation models often are developed with the goal 

of using them for predicting future effects of management 

strategies. The study-area model was initially used in what is 

referred to as a trial-and-error prediction mode. In this mode, 

future scenarios are designed in which the nitrate loading 

input to the model is varied according to a hypothetical 

set of management strategies that could be imposed. The 

locations and rates of loading over time are specified as input 

to the simulation model and the model predicts the resulting 

distribution of nitrate concentrations in the aquifer and the 

discharge of nitrate to the streams. The scenario results 

then are compared to assess whether management strategies 

succeeded in meeting water-quality goals. This is referred to 

as a trial-and-error procedure because often many simulations 

must be made to find management strategies that meet 

water quality goals. The results of the scenario simulations 

are presented here for later comparison to results of the 

simulation-optimization approach.

OR19-0049_fig14

Oregon

Department

of

Environmental

Quality data

n=109

Simulated

 0

5

10

15

20

25

30

n=1,398

Department

of Human

Services data

n=1,572

NITRATE, IN MILLIGRAMS N PER LITER

n=109 Number of observations

90th percentile

75th percentile

Mean

Median

25th percentile

10th percentile

Sample not included in 

 the 10-90 percentile 

 range

EXPLANATION

Figure 14. Measured and simulated nitrate concentrations 

in the La Pine, Oregon, study area, 1999.
Nitrate Fate and Transport Simulation Models 37

Figure 15. Simulated nitrate concentrations near the water table in the La Pine, Oregon, study area, 1999.

OR19-0049_fig15

Base 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

Publication projection is 

Lambert Conformal Conic

Standard parallels 42o20' and 44o40', 

central meridian -120o30'

Datum is NAD83

R 9 E R 10 E R 11 E

T 21 S

T 22 S

T 23 S

T 24 S

T 20 S

T 19 S

DESCHUTES CO.

L A K E C O .

KLAMATH CO.

Spring

River

Fall

River

Crescent Creek

Creek

Little

Deschutes

DESCHUTES

RIVER

River

Paulina

Long Creek

5 KILOMETERS

0 5 MILES

0

Wickiup

Reservoir

97

97

31

31

Burgess Road

Drive

Finley

Butte

Road

Road

South

Century

Spring

River

Masten Road

SunriverSunriver

La 

Pine

43°55'

50'

45'

121°40' 35' 30' 25' 121°20'

40'

43°35'

Extent of study area model

EXPLANATION

Nitrate concentration, in 

milligrams N per liter

0

5

10

20

50

100
38 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

The calibrated simulation model was used to predict 

the effects of eight decentralized wastewater management 

scenarios on ground-water and surface-water quality. Scenario 

number 1 was defined as the “base” scenario which would 

be used to evaluate the effects of a “no action” management 

strategy and serve as a benchmark for evaluating the benefits 

of other management alternatives. Under the base scenario, 

homes built on approximately 4,100 undeveloped lots 

remaining in 2000 used conventional on-site wastewater 

systems. Conventional systems were assumed to produce the 

same effluent concentration (46 mg N/L) and loading as used 

during the 1960-1999 model calibration simulations. The 

remaining seven scenarios were defined by Deschutes County 

using two primary management options to reduce nitrate 

loading: (1) a transferable development credit (TDC) program 

to reduce the number of on-site wastewater systems in the 

area, and (2) advanced treatment on-site wastewater systems. 

Deschutes County has modified the TDC program since the 

scenarios were defined; however the goal of the original TDC 

program was to reduce the number of on-site septic systems in 

the area by shifting, or transferring, development to a receiving 

area served by a community sewer and water system. 

Advanced treatment on-site wastewater systems capable 

of removing nitrogen were extensively field-tested as part of 

the La Pine NDP (Barbara Rich, ODEQ, written commun., 

2003). Based on data from the field testing program, three 

levels of nitrogen reduction performance were evaluated using 

the simulation model. Three reduction levels, expressed as a 

percentage of the nitrate concentration in effluent reaching 

the water table from conventional systems (46 mg N/L), 

were used: 57, 78, and 96 percent. These reductions would 

be equivalent to effluent nitrate concentrations at the water 

table of 20, 10, and 2 mg N/L, respectively. Four additional 

scenarios were defined where advanced treatment systems 

were combined with a TDC program in which development of 

1,500 lots was moved to a receiving area served by centralized 

sewer.

Additional assumptions for the scenarios included:

•	 Development would continue at the 1990–99 rates of 

250 homes per year until full build-out (in 2019);

•	 All new homes would use advanced treatment on-site 

wastewater systems; lots where new homes were built 

each year were randomly determined;

•	 Locations of existing homes where conventional 

on-site wastewater systems were upgraded to advanced 

treatment systems were randomly selected in each 

year; 

•	 Upgrades were done at a rate of 94 per year (based on 

historic rate of system repair permits issued by county) 

until all on-site wastewater systems had been upgraded 

(2056); and

•	 For scenarios involving TDCs, 1,500 lots were 

randomly selected from a pool of 2,600 candidates 

provided by Deschutes County. 

Descriptions of the scenarios and the resulting nitrate 

loading rates are summarized in table 6. 

Each scenario was simulated for 140 years beginning 

in 2000. The historical nitrate loading rates (1960-99) are 

compared with estimated future loading rates (2000-2139) for 

each scenario in figure 16. 

Under the base scenario (number 1, table 6), all available 

lots (as of 2000) were projected to be developed by 2019, and 

figure 16 shows that estimated nitrogen loading continued to 

increase until that time, after which loading remained constant 

at a rate of 147,000 lb/yr. Simulation of the base scenario 

showed that nitrate concentrations continued to increase 

long after the maximum loading rate was reached in 2019. 

This is because the amount of nitrate entering the groundwater system will exceed the amount leaving the system until 

equilibrium is reached. Equilibrium occurs when loading 

is balanced by the sum of the rates of denitrification and 

discharge of nitrate to streams. The simulation model suggests 

that it could take more than 140 years to reach equilibrium for 

scenario 1. At equilibrium, 78 percent of nitrogen entering the 

system (114,000 lb/yr) will be transported to the suboxic part 

of the aquifer and removed by denitrification. The remaining 

22 percent (33,000 lb/yr) will be transported into the nearstream areas adjacent to the Little Deschutes and Deschutes 

Rivers. This should be considered an upper bound on the 

amount of nitrate reaching the rivers because the current study 

area model cannot account for processes (denitrification, 

plant uptake, microbial uptake) that may remove nitrate from 

ground water before it discharges to the rivers. 

The simulated nitrate concentrations for scenario 1 

exceeded 10 mg N/L over large areas prior to equilibrium, 

however at equilibrium, the nitrate concentration near the 

water table averaged more than 10 mg N/L over areas totaling 

about 9,400 acres (table 7, fig. 17). 

The results of the other seven scenarios showed 

improved water-quality conditions for each level of increased 

nitrate loading reduction. Simulations indicate significant 

improvements in overall ground-water quality for all 

performance levels tested if all new homes use advanced 

treatment on-site wastewater systems and all existing homes 

replace conventional on-site wastewater systems with 

advanced treatment on-site wastewater systems. For the 

20 mg N/L nitrate performance level without a TDC program 

(scenario 3) the area where nitrate concentrations exceed the 

10 mg N/L drinking-water standard is reduced by 80 percent 

to about 1,900 acres (table 7, fig. 18) at equilibrium. A 

further reduction in effluent nitrate concentration, to 10 mg 

N/L, (scenarios 5 and 6) reduced the area where simulated 

equilibrium nitrate concentrations exceed 10 mg N/L to 

less than 700 acres. The maximum reduction, to 2 mg N/L, 

(scenarios 7 and 8) resulted in only a few small areas with 

concentrations greater than 10 mg N/L, and those areas were 

related to nonresidential loading (RV parks) not reduced under 

the scenarios. 
Nitrate Fate and Transport Simulation Models 39

Table 6. Summary of eight on-site wastewater management scenarios tested with the study-area model in 

the La Pine, Oregon, study area.

[TDC lots: assumes 1,500 lots eremoved from pool of candidate lots. Locations were randomly selected from 2,600 possible 

lots. Retrofitting: assumes all on-site systems serving homes built prior to 2000 are retrofitted with advanced systems with 

performance equal to systems installed in new homes. Retrofits are made at a rate of 94 per year until completed in 2057. 

Abbreviations: TDC, transferable development credit; lb/yr, pound per year; mg N/L, milligram nitrogen per liter]

Scenario

Effluent 

recharge 

concentration 

(mg N/L)

Percent 

reduction 

from standard 

systems

TDC lots Retrofitting

Peak 

loading

(lb/yr)

Peak 

loading 

year

Final 

loading

(lb/yr)

1 46 0 No No 147,000 2019 147,000

2 46 0 Yes No 125,000 2013 125,000

3 20 57 No Yes 93,000 2019 65,300

4 20 57 Yes Yes 88,100 2013 55,900

5 10 78 No Yes 81,200 2000 34,200

6 10 78 Yes Yes 81,200 2000 29,500

7 2 96 No Yes 80,400 2000 9,200

8 2 96 Yes Yes 80,400 2000 8,300

YEAR

NITRATE LOADING, IN POUNDS PER YEAR X 1,000

1960 1980 2000 2020 2040 2060 2080 2100 2120 2140

1−Base scenario 

2−TDCs only 

3−20 mg/L 

4−20 mg/L + TDCs

5−10 mg/L 

6−10 mg/L + TDCs

7−2 mg/L 8−2 mg/L + TDCs

Historical

 0

40

80

60

20

120

100

160

140

Figure 16. Historical nitrate loading from on-site 

wastewater systems and eight nitrate loading scenarios 

tested with the study-area model. 

See table 6 for scenario descriptions. (mg/L, milligram per 

liter; TDC, transferable development credit.)
40 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

Each scenario in which conventional systems were 

replaced with advanced treatment systems indicated that 

ground-water nitrate concentrations would peak between 

2027 and 2053 (about 25–50 years after loading reduction 

begins) and then decrease to their equilibrium levels. Less 

time is required to reach the peak and subsequent equilibrium 

concentrations for scenarios with greater reductions in loading 

(table 7). 

The TDC program, as defined for these scenarios, 

reduced the number of homes that would use on-site 

wastewater systems by shifting or transferring development 

of new homes to a receiving area with centralized wastewater 

treatment. Using conventional on-site systems, the 1,500 

homes selected for the TDC program in the scenarios would 

contribute about 21,000 lb of nitrogen annually, with each 

home adding 14 lb/yr to the aquifer. This equals nearly 

14 percent of the total base scenario loading of 147,000 lb/yr. 

For other scenarios, in which advanced treatment on-site 

wastewater systems are used, the effect of a program like this 

is diminished. For example, if advanced treatment systems 

reduce effluent concentration from 46 to 10 mg N/L, the 

annual loading per home is reduced from 14 to 3 lb/yr. Under 

this scenario, removing 1,500 homes only reduces total 

loading by 4,500 lb/yr (3 percent). The relative effectiveness 

of the TDC program in reducing loading at various levels 

of assumed advanced on-site system performance is shown 

in figure 16. Development transfer programs (like the TDC 

program) might be most effectively applied in high-sensitivity 

areas, where advanced treatment on-site wastewater systems 

cannot meet loading reduction goals. 

Results of trial-and-error simulations show that the 

capacity to receive on-site wastewater system effluent and 

maintain satisfactory water-quality conditions is variable 

within the study area. The capacity of any area to receive 

on-site wastewater system effluent is related to many factors, 

including the density of homes, presence of upgradient 

residential development, ground-water recharge rate, groundwater flow velocity, and thickness of the oxic part of the 

aquifer.

Each scenario tested with the simulation model was 

limited to management strategies that were applied uniformly 

in the study area. Typically, this is how simulation models 

must be used because as the size and complexity of the waterquality management problem increases, the decision makers’ 

ability to design scenarios with management strategies that 

reflect variability in the loading capacity of individual areas 

is diminished. Uniform management strategies, such as 

requiring all on-site wastewater systems to meet the same 

nitrate reduction standards, may be more costly than variable 

management strategies that account for the variability in the 

nitrate loading capacity of the ground-water system. Uniform 

strategies that are stringent enough to protect water quality 

in some areas may be more than is needed to protect water 

quality in other areas. The simulation model represents 

the hydrogeologic and chemical processes that cause this 

variability, and by adding optimization capability to the model, 

the nitrate loading capacity of the ground-water system can be 

determined for individual management areas. 

Table 7. Summary of model simulation results for eight on-site wastewater management scenarios tested with the study-area model in 

the La Pine, Oregon, study area.

[Acres are based on simulated nitrate concentrations in the uppermost active (saturated) model cells. Equilibrium acreages are based on simulated nitrate 

concentrations at end of simulation (2139). Peak acreages are based on simulated nitrate concentrations in the year when the highest concentrations occur for 

each scenario. Abbreviations: TDC, transferable development credit; mg N/L, milligram nitrogen per liter; lb/yr, pound per year; mg/L, milligram per liter; 

Symbol: –, no data]

Scenario

Effluent 

recharge 

concentration 

(mg N/L)

TDC

Equilibrium Peak

Equilibrium 

nitrate 

loading

(lb/yr)

Reduction in 

loading from 

scenario 1 

(percent)

Acres greater 

than 10 mg/L 

at equilibrium

Reduction 

from 

scenario 1

(percent)

Acres greater 

than 10 mg/L 

at peak

Reduction from 

scenario 1

(percent)

Peak 

concentration 

year

1 46 No 147,000 – 9,398 – 9,398 – 2139

2 46 Yes 125,000 15 7,317 22 7,317 22 2139

3 20 No 65,300 56 1,909 80 1,944 79 2053

4 20 Yes 55,900 62 1,433 85 1,456 85 2049

5 10 No 34,200 77 700 93 1,233 87 2038

6 10 Yes 29,500 80 619 93 1,153 88 2035

7 2 No 9,200 94 115 99 975 90 2027

8 2 Yes 8,300 94 109 99 975 90 2027
Nitrate Fate and Transport Simulation Models 41

OR19-0049_fig17

Base 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

Publication projection is 

Lambert Conformal Conic

Standard parallels 42o20' and 44o40', 

central meridian -120o30'

Datum is NAD83

R 9 E R 10 E R 11 E

T 21 S

T 22 S

T 23 S

T 24 S

T 20 S

T 19 S

DESCHUTES

L A K E

KLAMATH

Spring

River

Fall

River

Crescent Creek

Creek

Little

Deschutes

DESCHUTES

RIVER

River

Paulina

Long Creek

5 KILOMETERS

0 5 MILES

0

Wickiup

Reservoir

97

97

31

31

Burgess Road

Drive

Finley

Butte

Road

Road

South

Century

Spring

River

Masten Road

SunriverSunriver

La 

Pine

43°55'

50'

45'

121°40' 35' 30' 25' 121°20'

40'

43°35'

Extent of study area model

EXPLANATION

Nitrate concentration, in 

milligrams N per liter

0

5

10

20

50

100

Figure 17. Simulated equilibrium ground-water nitrate concentrations near the water table for the base scenario (scenario 1, 

table 6) in the La Pine, Oregon, study area. (Colors indicate maximum nitrate concentration in vertical dimension.)
42 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

OR19-0049_fig18

Base 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

Publication projection is 

Lambert Conformal Conic

Standard parallels 42o20' and 44o40', 

central meridian -120o30'

Datum is NAD83

R 9 E R 10 E R 11 E

T 21 S

T 22 S

T 23 S

T 24 S

T 20 S

T 19 S

DESCHUTES

L A K E

KLAMATH

Spring

River

Fall

River

Crescent Creek

Creek

Little

Deschutes

DESCHUTES

RIVER

River

Paulina

Long Creek

5 KILOMETERS

0 5 MILES

0

Wickiup

Reservoir

97

97

31

31

Burgess Road

Drive

Finley

Butte

Road

Road

South

Century

Spring

River

Masten Road

SunriverSunriver

La 

Pine

43°55'

50'

45'

121°40' 35' 30' 25' 121°20'

40'

43°35'

Extent of study area model

EXPLANATION

Nitrate concentration, in 

milligrams N per liter

0

5

10

20

50

100

Figure 18. Simulated equilibrium ground-water nitrate concentrations near the water table for 20 milligrams N per liter advanced 

treatment on-site wastewater systems (scenario 3, table 6) in the La Pine, Oregon, study area. (Colors indicate maximum nitrate 

concentration in vertical dimension.)
Nitrate Loading Management Model 43

Nitrate Loading Management Model

The study-area simulation model was linked to 

optimization methods to produce the Nitrate Loading 

Management Model (NLMM). This section includes a 

description of how the management model was formulated, 

how it is solved, and how it can be used to evaluate alternative 

strategies for managing nitrate loading to the shallow groundwater system.

Formulation of Nitrate Loading Management 

Model

To use optimization methods, management goals must 

be formulated into a mathematical structure. The structure 

of an optimization problem has three components that must 

be defined: decision variables, an objective function, and 

constraints. 

One inherent value of the optimization approach is that 

its use requires decision makers and stakeholders to quantify 

planning goals and objectives as well as environmental 

and other constraints. The La Pine NLMM was formulated 

through a close collaboration with Deschutes County planners 

and resource managers. The management objective of the 

NLMM, to minimize the reduction from base scenario 

nitrate loading from on-site wastewater systems, reflects the 

goals of: (1) allowing as many platted lots as possible to be 

developed under current development goals and policies, 

and (2) minimizing the required reduction in nitrate loading 

from existing on-site systems. Because costs are associated 

with reducing loading either by not allowing development 

or using advanced treatment on-site wastewater systems, the 

management objective can be simply stated as “minimize 

the cost” of meeting water-quality goals. The water-quality 

goals are constraints on the management model and have 

direct and quantifiable effects on the solution (or cost) and 

their values often reflect regulatory requirements or economic 

and community values. If regulatory standards apply, the 

process of setting water-quality constraints is straightforward. 

However, if less well defined economic or community values 

are to be the basis for water quality constraints, then a tradeoff analysis is a common process for finding the balance 

(or evaluating the relations) between costs and economic, 

aesthetic, health, environmental, and other benefits. The 

management model can assist in this process by quantifying 

the relationships between constraints and management 

objectives.

The NLMM was formulated mathematically to minimize 

the reduction in nitrate loading (or cost of reduction) from 

base scenario loading, subject to constraints on minimum 

reductions in ground-water nitrate concentrations (relative to 

base scenario), minimum reductions in ground-water discharge 

nitrate loading to streams (relative to base scenario), and 

minimum and maximum loading reductions for existing and 

future homes in each management area.

The decision variables in the La Pine NLMM were the 

reduction in nitrate loading (relative to base scenario nitrate 

loading rates) that would be needed to maintain or achieve 

desired water-quality conditions. See table 6 for a description 

of the base scenario. (Nitrate loading rates were specified 

using metric units [kilograms per day] in the NLMM and 

these units are used in the following description of the model 

for consistency.) To define decision variables, the study area 

was divided into 97 management areas. The basic unit for 

the management areas was a 160-acre area (quarter-section) 

based on the Public Land Survey System. With the goal of 

making the base scenario loading in each management area 

approximately equal, as many as four quarter-sections (640 

acres) were combined to form management areas in areas 

with lower lot density. The average number of lots in each 

management area was approximately 100 (fig. 19). The 

model had 194 decision variables, NRi,j, because two decision 

variables were defined for each management area: loading 

reduction (kilograms per day, kg/d) for existing homes and 

loading reduction for future homes. About 350 homes and lots 

were not included in the NLMM because they were in low 

density areas. These homes and lots represented only 7 kg/d of 

the total base scenario loading. 

The objective function of the NLMM was to minimize 

the reduction in loading (from the base scenario), and is given 

as

minimize u NR

NM

j i j

j

NT

i

NM

, , ==

∑∑ 11

where

is the number of management areas,

NT is the number of management area types, which

can be either 1 for existing homes or 2 for

future homes,and

is a dimensionless unit cost factor for nitrate

load

uj

ing reduction from existing or future homes.

(1)

Constraints on Ground-Water Nitrate 

Concentration

The value of the objective function was limited by a set of 

constraints on minimum reductions (relative to base scenario) 

in ground-water nitrate concentrations, minimum reductions 

(relative to base scenario) in ground-water discharge nitrate 

loading to streams, and minimum and maximum loading 

reductions to the aquifer for existing and future homes in each 

management area.
44 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

OR19-0049_fig19

Base 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

Publication projection is 

Lambert Conformal Conic

Standard parallels 42o20' and 44o40', 

central meridian -120o30'

Datum is NAD83

R 9 E R 10 E R 11 E

T 21 S

T 22 S

T 23 S

T 24 S

T 20 S

T 19 S

DESCHUTES

L A K E

KLAMATH

Spring

River

Fall

River

Crescent Creek

Creek

Little

Deschutes

DESCHUTES

RIVER

River

Paulina

Long Creek

5 KILOMETERS

0 5 MILES

0

Wickiup

Reservoir

97

97

31

31

Burgess Road

Drive

Finley

Butte

Road

Road

South

Century

Spring

River

Masten Road

SunriverSunriver

La 

Pine

43°55'

50'

45'

121°40' 35' 30' 25' 121°20'

40'

43°35'

Extent of study area model

EXPLANATION

Management area

Location of homes

Existing (1999)

Future

Ground-water quality 

constraint location

See figure 20

for details

Figure 19. Locations of management areas and ground-water nitrate concentration constraint locations in the Nitrate Loading 

Management Model for the La Pine, Oregon, study area.
Nitrate Loading Management Model 45

Minimum reductions in ground-water nitrate 

concentration constraints were set at 339 locations in the 

simulation model. The minimum reduction values were 

determined by simulating the base scenario equilibrium 

concentration at each location and computing the reduction 

necessary to meet a specified concentration value. For 

example, if the base scenario nitrate concentration were 

calculated by the simulation model to be 23 mg N/L and the 

maximum allowable concentration was 7 mg N/L, then the 

minimum reduction constraint at that location would be 16 mg 

N/L.

Constraints were set at one or two locations in the 

simulation model for each management area. One point was 

selected to be sensitive to loading from existing homes and 

another point was selected to be sensitive to loading from 

future homes in the management area. The most sensitive 

locations for existing and future loading were determined by 

simulating loading only from existing or future homes in a 

management area and determining the location (cell) in the 

model where the highest simulated concentrations occurred. 

If the same location were most sensitive to both existing and 

future loading from a management area, then only one location 

was set for that management area. 

The model simulates three-dimensional flow and 

transport; therefore, nitrate concentrations can vary with 

depth. To account for the variation in concentration with 

depth, shallow and deep constraints were specified at most 

locations. The shallow constraint was specified 5–10 ft below 

the water table (an average of 20 ft below land surface) and the 

deep constraint was specified immediately above the suboxic 

boundary (an average of 50 ft below land surface). Four 

unique constraints (one shallow and one deep for both existing 

and future homes) specified in each management area for each 

of the 194 locations where concentration was constrained 

would have required specified constraints in 388 model cells. 

However, forty-nine locations were eliminated because the 

same location was sensitive to loading from both existing and 

future homes, and some cells were eliminated because the oxic 

ground-water layer was only two cells thick and only one cell 

was used. In all, there were 339 concentration constraint cells 

at 174 locations in the NLMM (fig. 19).

Mathematically, the minimum ground-water 

concentration reduction constraints were specified as

CRmink ≤ CRm,k,

. (2)

The minimum concentration reduction was computed as

CRminm,k = Csqm,k – Cmaxk, (2a)

The effect of loading from the 350 homes and lots not 

included in the NLMM was accounted for by simulating 

equilibrium concentrations throughout the model that resulted 

from base scenario loading at these 350 locations. The 

resulting concentrations at each constraint site were used to 

adjust the CRminm,k values at each site. 

Constraints on Discharge of Nitrate from Ground 

Water to Streams

The simulation model results indicated that nitrate from 

on-site wastewater systems can reach streams through groundwater discharge. Because the addition of nitrate to streams 

could have a deleterious effect on stream quality, the NLMM 

was configured to allow constraints on the amount of nitrate 

discharged to streams. Constraints on minimum reduction 

in discharge of nitrate from ground water to streams were 

specified for 14 reaches on the Deschutes and Little Deschutes 

Rivers. The minimum reduction values were determined 

by simulating the base scenario equilibrium ground-water 

discharge nitrate loading to each reach and computing the 

reduction necessary to meet a specified discharge constraint. 

For example, if the base scenario ground-water discharge of 

nitrate to a reach was determined to be 3 kg/d and the specified 

discharge constraint (maximum allowable discharge of nitrate) 

was 1 kg/d, then the minimum reduction for that reach would 

be 2 kg/d. Mathematically, the minimum reduction in groundwater discharge nitrate loading to stream constraints were 

specified as

DRminr ≤ DRr (3)

The minimum discharge loading reduction was computed as

DRminr = Dsqr – Dmaxr (3a)

Constraints on Reduction of Nitrate Loading

Minimum and maximum nitrate loading reduction 

constraints were specified in the NLMM. Mathematically the 

constraints were expressed as

NRminj ≤ NRj ≤ NRmaxj, (4)

where NRminj andNRmaxj are minimum and maximum 

loading reduction constraints for existing (j=1) and future 

(j=2) homes (units of kg/d). The loading reduction constraints 

establish minimum or maximum loading reduction constraints 

on either existing or future homes.
46 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

Response-Matrix Technique for Solution of 

Nitrate Loading Management Model

The optimization method used to solve the NLMM is 

based on a widely applied technique for solving ground-water 

management problems called the response-matrix technique. 

The assumption that must be satisfied to use this technique 

is that the nitrate concentrations at each constraint site are a 

linear function of the loading rates in each management area. 

By assuming linearity, determining the nitrate concentration 

or ground-water discharge loading rate to streams is possible 

at any location where a constraint is specified by summing 

the contribution of loading in each management area to the 

concentration or discharge at that location. The responsematrix technique is described in detail by Gorelick and others 

(1993) and Ahlfeld and Mulligan (2000). The responsematrix approach has been used to solve ground-water wastemanagement problems similar to the problem evaluated in this 

work by Moosburner and Wood (1980), Gorelick (1982), and 

Gorelick and Remson (1982). 

The assumption that simulated nitrate concentration 

and discharge nitrate loading to streams is a linear function 

of loading to the aquifer was tested through a series of 

simulations. The loading rate to a single management area 

was varied in each simulation and the computed nitrate 

concentrations at several constraint locations (model cells) 

were recorded. When loading rate was plotted against 

concentration for each location, the relation was shown to be 

linear. 

One requirement of this technique is that the effect of the 

sources of nitrate loading on the ground-water velocity field 

is known because the mass loading of nitrate is the product 

of the nitrate concentration and hydraulic loading rate of the 

on-site wastewater systems. Previous workers have assumed 

that for concentrated pollutant sources, the source-water flow 

rate (hydraulic loading rate) has a negligible influence on 

the ground-water velocity field (Gorelick, 1982; Gorelick 

and Remson, 1982). This also is the assumption used for this 

study, where recharge from on-site wastewater systems to, and 

domestic water-supply well withdrawals from, the shallow 

aquifer are nearly equal, and both are small relative to natural 

recharge rates. Therefore, the ground-water velocity field 

without the influence of on-site wastewater systems as sources 

of recharge was used to simulate transport. 

Use of the response-matrix technique requires that 

response functions for ground-water nitrate concentrations and 

ground-water discharge of nitrate to streams are calculated at 

each of the 339 model cells where concentration constraints 

were specified and the 14 stream reaches where ground-water 

discharge of nitrate was constrained. The response functions 

were calculated by making 194 simulations, one for each 

decision variable (97 management areas, each with decision 

variables for loading reduction for existing and future homes). 

The response to base scenario loading from either existing or 

future homes in a single management area was simulated. The 

spatial distribution of loading within each management area 

was at the resolution of a model cell—500 ft in the x- and ydimensions. The purpose of these simulations was to compute 

the change in concentration at each constraint location caused 

by a unit change in loading at the location represented by 

the decisions variable. Calculation of the unit change in 

concentration was made using initial nitrate concentrations 

of zero and simulating equilibrium concentrations. The base 

scenario loading rates (kg/d) for existing and future homes 

in management area i and type j were defined as Nsqi,j. 

Concentration response coefficients, rci,j,m,k, were computed as

rc

C

N i j m k

m k

i j

sq

sq , , ,

,

,

= , (5)

Ground-water discharge nitrate loading response coefficients, 

rdi,j,r, for each stream reach were computed as

rd D

N i j r

r

i j

sq

sq , ,

,

= , (6)

Units of milligrams per liter per kilogram per day were 

used for concentration response coefficients and groundwater discharge nitrate to stream response coefficients are 

dimensionless. The magnitude of the response coefficient 

is directly proportional to the sensitivity of the nitrate 

concentration or discharge nitrate loading at a constraint site to 

loading in a management area. This relation is illustrated for 

the concentration constraint site 31-E-S (fig. 20), which is the 

shallow concentration constraint site selected for its sensitivity 

to loading from existing homes in management area 31. 

Management area 31 is adjacent to the west bank of 

the Little Deschutes River, approximately 1 to 1.5 mi north 

of Burgess Road (fig. 20). The linear response coefficients, 

rci,j,k,m, for several nearby management areas and types of 

homes (existing or future), show that loading in management 

areas 31, 32, 50, 51, 52, and 53 significantly affect the 

nitrate concentration at constraint location 31-E-S (table 7). 

Management areas 48, 49, 63, and 64 have lesser effects. 

Response coefficient values (table 8) indicate that the nitrate 

concentration at site 31-E-S is not a simple function of the 

number of homes in management area 31 or the distance 

from the site to adjacent management areas. Table 8 shows 

that the largest influence on concentration at point 31-E-S is 

from future homes in area 31, followed by existing homes 

in areas 31 and 50, future homes in areas 50 and 51, and 

existing homes in area 51. Response coefficients indicate that 

management area 32, the closest area to site 31-E-S, has less 

influence on concentration than more distant areas because of 

the direction of ground-water flow through the area (toward 

the Little Deschutes River) and the location of site 31-E-S in 

relation to homes in area 32.
Nitrate Loading Management Model 47

Assuming the system is linear, the reduction in nitrate 

concentration, CRm,k , at constraint location m and depth k can 

be calculated with the concentration response coefficients by 

summation of the individual concentration reductions caused 

by reductions in loading to existing and future homes in each 

management area. The summation is written as

CRm k rci j m kNRi j

j

NT

i

NM

, , , , , = . ==

∑∑ 11

(7)

Base 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

Publication projection is Lambert Conformal Conic

Standard parallels 42o20' and 44o40', central meridian -120o30'

Datum is North American Datum of 1983 R 9 E R 10 E

R 10 E

T 21 S

T 22 S

Little

Deschutes

River

0 1 MILE

0 1 KILOMETER

43°

45'

121°35' 32'30” 121°30'

43°

42'

30”

EXPLANATION

Management area

Line of equal simulated 

water table elevation— In 

feet above NAD83. 

Contour interval is 10 feet

Location of homes

Existing (1999)

Future

Nitrate concentration 

constraint location 31-E-S

4220

97

Burgess Road

4200

4260

4250

4240

4220

4230

4210

4220

4210

4200

4190

67

66 56

58

57 61

62

65 54

74

64

63

53

52

51 50

32

31

26

35

55 36

25

49 48

44

46 45

47

43

41

42

30

29

27

Similarly, the reduction in ground-water discharge nitrate 

loading to streams, DRr, can be calculated with the discharge 

loading response coefficients by summation of the individual 

discharge loading reductions caused by reductions in loading 

to existing and future homes in each management area. The 

summation is written as

DRr i rd j rNRi j

j

NT

i

NM

= ==

∑∑ , , , . 11

(8)

Figure 20. Locations of management areas near Burgess Road and management area 31 in the La Pine, Oregon, study area.
48 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

These summations include 194 terms because there 

are 97 management areas (NM) and 2 loading types (NT), 

existing and future homes, in each management area. In each 

summation, however, many terms are equal to zero because the 

response coefficients are zero; this occurs when concentration 

or discharge loading at a constraint location does not affect 

loading in a management area. 

The response matrix was prepared by running one 

simulation for each decision variable in which loading 

was applied only at the locations defined for that decision 

variable (for example, existing homes in management area 

31). The loading rate and resulting concentrations at each of 

the constraint locations were used to compute the response 

coefficient at each constraint location using equations 5 and 6. 

A set of utility programs and scripts were developed to run the 

194 simulations, extract the simulated concentration values at 

constraint sites, compute the response coefficients, and format 

the coefficients into a matrix for input to the optimization 

program.

Table 8. Response coefficients relating the effects of loading in 

nearby areas to the nitrate concentration at a constraint location 

in management area 31 in the La Pine, Oregon, study area.

[Management area: 31—site 31-E-S in figure 20. Response coefficients are 

in milligrams per liter per kilogram per day loading]

Management 

area Type Response coefficient

31 Future 14.85

50 Existing 11.10

31 Existing 9.92

50 Future 8.85

51 Future 3.49

51 Existing 3.07

32 Future 1.71

52 Future 1.56

52 Existing 1.31

53 Existing 0.75

32 Existing 0.74

53 Future 0.50

63 Existing 0.15

63 Future 0.15

64 Existing 0.03

64 Future 0.02

49 Existing 0.01

49 Future 0.01

48 Future 0.01

48 Existing 0.01

Response coefficients are the link between the simulation 

model and the NLMM. The response coefficients are utilized 

by the NLMM by replacing CRm,k (equation 2) with the righthand side of equation 7 and replacing DRr

 (equation 3) with 

the right hand side of equation 8. The constraints for reduction 

in nitrate concentration and discharge loading to streams are 

then written as

CR m k rci j m kNRi j

j

NT

i

NM

min , , , , , ≤ , ==

∑∑ 11

(9)

and

DR r i rd j rNRi j

j

NT

i

NM

min ≤ ==

∑∑ , , , . 11

(10)

Equations 9 and 10 replace equations 2 and 3 in the 

NLMM.

The modified NLMM, defined by equations 1

(objective function), 4 (loading constraints), 9 (concentration 

constraints), and 10 (discharge loading constraints), constitutes 

a linear program. This program was solved using the “What’s 

Best!” optimization program (LINDO Systems, 2003). This 

set of solvers is implemented as an add-in to Microsoft Excel. 

The response matrix and constraint definitions are defined 

in the spreadsheet using special functions provided in the 

“What’s Best!” add-in. The NLMM can mathematically 

search for the minimum nitrate loading reductions for existing 

and future homes in each management area that satisfy the 

constraints on nitrate concentrations in ground water and 

ground-water discharge loading of nitrate to streams. The 

program also identifies management problems that are not 

feasible. This occurs when at least one constraint cannot be 

met with any combination of decision variable values.

Application of Model

In the application of the La Pine NLMM described 

here, the values of constraints were varied in a sensitivity 

analysis to explore relations between constraint values and 

optimal loading rates to the aquifer. The three constraint 

options used in the NLMM were varied and optimal solutions 

were computed to demonstrate how decision makers might 

use the model in a trade-off analysis to determine the effect 

of constraint values on the objective of minimizing nitrate 

loading reductions. The NLMM also was used to evaluate the 

influence of a cost variable that accounts for differential costs 

in reducing nitrate loading from existing and future homes.
Nitrate Loading Management Model 49

Sensitivity of Optimal Solution to Water-Quality 

Constraints

Ground-Water Nitrate Concentration

Ground-water concentration constraints were specified 

at as many as two locations for each management area: one 

location was the most sensitive to loading from existing homes 

in that management area, and one location was most sensitive 

to loading from future homes. In 20 management areas, the 

same location was most sensitive for both existing and future 

homes. A shallow constraint was specified near the water 

table at each these 174 locations (fig. 19). In addition, where 

the oxic part of the aquifer was greater than 10-ft thick (163 

locations), a deep constraint was specified. As described in 

the discussion of the model formulation, the concentration 

constraint values, CRminm,k, were computed (equation 2a) as 

the minimum reduction from base scenario concentration that 

would be required to meet a specified concentration limit, 

Cmaxk

. 

Concentration constraints were varied to show the 

relative sensitivity of optimal loading solutions to constraints 

in different parts of the aquifer system. The total nitrate 

loading for existing and future homes in all 97 management 

areas was summed for each optimal solution to provide a 

basis for comparing the sensitivity of the solutions to values 

of the concentration constraints. The NLMM was solved 

using a range of 1 to 25 mg N/L for the maximum allowable 

concentration values for both the shallow (Cmaxs) and deep 

(Cmaxd) constraint locations. This range was selected because 

it includes most of the range of values that would constrain, 

or bind, the optimal solution. Binding constraints limit the 

amount of loading to the aquifer. Nonbinding constraints do 

not limit the loading to the aquifer system, because of their 

value or location. 

The NLMM was used to compute optimal nitrate 

loading rates for various combinations of shallow and deep 

concentration constraints (fig. 21). Optimal solutions were 

most sensitive to concentration constraints in the shallow part 

of the ground-water system. This result was expected because 

concentrations are greatest at the water table where loading 

occurs, and decrease with depth and distance downgradient 

from the source. The base scenario loading used in the 

NLMM was 190 kg/d (153,000 lb/yr) and, as expected, this 

is the optimal loading solution for the case where there are 

no constraints specified. Curve A (fig. 21) shows the effect 

on optimal loading when concentrations are constrained only 

in the deep part of the system. Deep constraints do not limit 

optimal loading until constraint values are less than 10 mg N/L 

(point C, fig. 21). Specifying a nitrate concentration constraint 

of 1 mg N/L in the deep part of the aquifer limited optimal 

loading to 110 kg/d (point D). Curve B shows the effect on 

optimal loading when concentrations are constrained only 

in the shallow part of the system. Shallow constraints limit 

optimal loading throughout the range, with loading limited 

to 17 kg/d at a constraint value of 1 mg N/L (point E) and 

168 kg/d at a constraint value of 25 mg N/L (point F). If a 

regulatory limit, such as the ODEQ “action level” of 7 mg N/L 

is applied to the shallow constraints (and no deep constraints 

are specified), optimal loading is reduced to 84 kg/d (point 

G)—a 56 percent reduction from base scenario loading. For 

comparison, if no shallow constraints are specified and 7 mg 

N/L is applied to the deep constraints, optimal loading is 

reduced very little, to 183 kg/d (point H). Loading sensitivity 

curves also are shown in figure 21 for shallow concentration 

constraint (Cmaxs) values of 5, 10, 15, 20, and 25 mg N/L; 

these curves were constructed by keeping the shallow 

constraint constant at the indicated value and varying the 

value of the deep constraint. For example, with shallow and 

deep constraints of 5 and 15 mg N/L respectively (point I), 

the optimal loading is 65 kg/d. For reference, the optimal 

loading for shallow and deep concentration constraints of 7 

and 3 mg N/L, respectively, is 84 kg/d and plots at point J in 

figure 21. The effects of more- or less-stringent concentration 

constraints on nitrate loading in the shallow and deep parts of 

the ground-water system can be analyzed using figure 21.

NITRATE CONCENTRATION CONSTRAINT, MILLIGRAMS NITROGEN PER LITER

0 5 10 15 20 25

OPTIMAL NITRATE LOADING, KILOGRAMS PER DAY

 0

50

100

150

200

250

Unconstrained loading = 190 kg/d

Cmaxs = 25

Cmaxs = 20

Cmaxs = 15

Cmaxs = 10

Cmaxs = 5

A Deep constraints only

B Shallow constraints only

C

D

E

F

G

I

J

H

EXPLANATION

Curve A—Relation between optimal loading and values of concentration 

 constraints at deep locations

Curve B—Relation between optimal loading and values of concentration 

 constraints at shallow locations

Points C-J—Optimal loading for combinations of shallow and deep 

 concentration constraints

 C—Deep constraints set to 10 mg N/L, no shallow constraints

 D—Deep constraints set to 1 mg N/L, no shallow constraints

 E—Shallow constraints set to 1 mg N/L, no deep constraints

 F—Shallow constraints set to 25 mg N/L, no deep constraints

 G—Shallow constraints set to 7 mg N/L, no deep constraints

 H—Deep constraints set to 7 mg N/L, no shallow constraints

 I—Shallow constraints set to 5 mg N/L, deep constraints set to 15 mg N/L

 J—Shallow constraints set to 7 mg N/L, deep constraints set to 3 mg N/L

 

Figure 21. Sensitivity of optimal loading solutions to 

ground-water nitrate concentration constraints in the La 

Pine, Oregon, study area. 
50 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

These results can be used to assess trade-offs in 

protection of water quality in the shallow and deep parts of the 

aquifer. Few wells in the La Pine area are open to the upper 

10 ft of the aquifer where the shallow constraints are specified. 

Using these results, decision makers could evaluate the cost, in 

terms of reduced loading, of protecting the shallowest part of 

the system (within 10 ft of the water table) to the same degree 

as the deeper part of the system (30–50 below the water table), 

where most domestic wells obtain water. 

Ground-Water Nitrate Discharge Loading to Streams

Fourteen stream reaches were defined where constraints 

could be applied to ground-water nitrate discharge loading. 

As described in the discussion of the model formulation, the 

discharge loading constraint values, DRminr, were computed 

(equation 3a) as the minimum reduction from base scenario 

discharge loading that would be required to meet a specified 

limit, Dmaxr

. The beginning and end of each reach were 

selected to coincide with confluences or road crossings. 

The purpose of this analysis was to demonstrate how the 

NLMM can be used to conjunctively manage water quality of 

ground-water and surface-water resources using optimization 

techniques. Discharge loading constraints were specified 

by setting a minimum percent reduction from base scenario 

ground-water discharge loading to streams. To determine 

the sensitivity to discharge loading constraints, the percent 

reduction was varied from 0 to 99 percent. Two sets of 

solutions were computed with NLMM: 1) with no constraints 

on ground-water concentrations, and 2) with shallow and 

deep ground-water concentration constraints of 7 and 3 mg 

N/L. Without concentration constraints, the discharge loading 

constraint directly limited the loading to streams as would 

be expected based on the formulation of the model (fig. 22). 

However, with the limitations imposed on loading by 

concentration constraints (7 and 3 mg N/L), discharge loading 

constraints had no effect on loading to streams until the 

constraint values were greater than about 50 percent reduction 

(fig. 22). At values greater than 50 percent, most loading 

reductions required to meet this constraint are needed in the 

management areas adjacent to streams, where, according to 

simulation results, shallow ground-water flow paths through 

the thin oxic part of the ground-water system are connected to 

the streams. With no concentration constraints and discharge 

loading to streams constrained to a minimum reduction of 

97 percent, the total loading to the system is about 100 kg/d. 

By comparison, constraining concentration to 7 and 3 mg N/L 

(shallow and deep, respectively) and constraining discharge 

loading reduction to at least 97 percent requires that total 

loading be reduced to about 40 kg/d (fig. 22). The additional 

reduction in total loading required to meet the ground-water 

concentration constraints, 60 kg/d, reflects the additional 

costs of management strategies designed to protect both 

ground-water quality and stream quality.

These results can be used to assess trade-offs in the 

protection of surface-water quality. To make decisions on 

the values of discharge constraints, the processes affecting 

nitrate as it is transported through the near stream and riparian 

environments need to be better understood (Hinkle and others, 

2007b). The simulation model uses simple assumptions 

regarding the fate and transport of nitrate in this part of the 

ground-water system and the estimates of ground-water 

discharge of nitrate to streams should be considered the upper 

limits of possible discharge. 

Sensitivity of Optimal Solution to Nitrate Loading 

Constraints

Minimum and maximum loading reduction constraints, 

NRmin and NRmax in equation 4, can be specified for existing 

or future homes in the NLMM. The constraints are computed 

as a percentage of the base scenario loading for the existing 

or future homes in each management area. For example, if 

base scenario loading for future homes in a management area 

is 2 kg/d, and the desired loading reduction is a minimum of 

25 percent of status-quo build-out loading, then NRmin would 

be equal to 0.5 kg/d. 

OR19-0049_fig 22

MINIMUM REDUCTION IN DISCHARGE LOADING TO STREAMS, PERCENT

0 20 40 60 80 100

NITRATE LOADING, KILOGRAMS N PER DAY

 0

50

100

150

200

Total loading, without concentration constraints

Total loading, with concentration constraints

Discharge to streams, without concentration constraints

Discharge to streams, with concentration constraints

EXPLANATION

Figure 22. Sensitivity of optimal loading solutions to 

constraints on the minimum reduction in ground-water 

discharge loading to streams in the La Pine, Oregon, study 

area.
Nitrate Loading Management Model 51

In this example, the loading constraints were used to 

evaluate how specifying minimum nitrate reduction standards 

for future homes would affect the overall loading and the 

proportion of loading from existing homes. As formulated, 

the NLMM also offers the ability to specify other constraints 

on loading reduction, such as maximum loading reduction 

for either future of existing homes. This type of constraint 

could be used to reflect the limitations of alternative 

treatment systems to reduce loading by more than a specified 

percentage. 

Optimal solutions were computed in which the minimum 

reduction constraint for future homes varied and all other 

constraints were constant. Ground-water concentration 

constraints of 7 and 3 mg N/L were specified for the shallow 

and deep sites, respectively, and ground-water discharge 

nitrate loading to streams was unconstrained. All existing 

on-site wastewater systems in the study area are assumed to 

discharge effluent with 46 mg N/L nitrate, and this was the 

performance level used to compute the base scenario loading 

of 190 kg/d. Maximum loading reduction constraints were 

96 percent of base scenario loading for both existing and 

future homes; this constraint reflects the assumption that 

the best attainable on-site wastewater system performance 

is 2 mg N/L. The minimum loading reduction constraint 

for future homes was specified on the basis of assumed 

performance standards for on-site wastewater systems that 

ranged from 2 to 46 mg N/L concentration in effluent that 

recharged the ground-water system. Thus, in this series of 

solutions, the 46 mg N/L performance level was the equivalent 

of zero reduction from base scenario loading. The other 

performance levels used were 2, 10, 20, and 30 mg N/L, 

corresponding to minimum loading reductions of 96, 78, 57, 

and 35 percent, respectively (fig. 23). 

Under base scenario conditions, existing and future 

homes will contribute 104 and 86 kg/d, respectively, to the 

total loading of 190 kg/d. The results of this analysis show that 

as the minimum loading reduction constraint for future homes 

is increased from 0 to 96 percent, optimal loading from future 

homes decreases from 38 to 3 kg/d and the associated total 

loading decreases from 84 to 58 kg/d. 

The reduction in loading to the aquifer from future homes 

implemented using a minimum reduction constraint allows 

higher loading rates to be maintained from existing homes 

while still meeting concentration constraints. If no loading 

reduction is required for future homes, loading from existing 

homes will have to be reduced by 56 percent (from 104 to 

46 kg/d). By requiring improved performance of on-site 

wastewater systems in future homes, there is less need for 

loading reduction from existing homes. For example, if a 

96 percent reduction requirement is imposed for future homes, 

loading from existing homes would have to be reduced by only 

47 percent to meet ground-water concentration constraints 

of 7 and 3 mg N/L for shallow and deep parts of the aquifer. 

If it were less costly to reduce loading from future homes 

by installing advanced treatment on-site wastewater systems 

than requiring existing homes to be retrofitted with advanced 

treatment on-site wastewater systems, then using this type of 

constraint in the NLMM would allow planners to incorporate 

specific on-site wastewater system performance standards for 

future homes into the optimal solution.

Sensitivity of Optimal Solution to Cost Factors

The simulation model accounts for the physical and 

geochemical complexities of the ground-water system and 

this information is available to the NLMM through the 

response coefficients that were computed using the study 

area simulation model. Other variables (external to the 

simulation model) also can be important in determining the 

optimal nitrate loading solution. The cost of implementing 

management strategies is the most common external variable 

that affects management decisions. In many optimization 

problems, the objective is to minimize the cost of satisfying 

the constraints on the problem.

OR19-0049_fig 23

MINIMUM REQUIRED REDUCTION IN NITRATE LOADING 

FROM NEW HOMES, IN PERCENT

0 35 57 78 96

OPTIMAL NITRATE LOADING, IN KILOGRAMS N PER DAY

0

20

40

60

80

100

 

Existing homes

Future homes

Figure 23. Sensitivity of optimal solution to minimum 

decentralized wastewater treatment performance 

standards for future homes in the La Pine, Oregon, study 

area.
52 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

Reducing loading from on-site wastewater systems can 

incur significant costs whether by limiting housing density 

or installing and maintaining advanced treatment on-site 

wastewater systems. The cost factor of reducing nitrate 

loading, uj in equation 1 of the NLMM, can be used to account 

for differences in cost for existing and future homes. If the 

cost of nitrate reduction varies among decision variables in 

other ways, such as geographically, cost variables can be 

specified for individual management areas.

The sensitivity of the La Pine NLMM to differences in 

the cost of reducing nitrate loading from existing and future 

homes was evaluated by computing the optimal loading 

solutions for a range of unit cost ratios. Arbitrary cost factors 

were specified that resulted in ratios of existing to future unit 

costs ranging from 0.11 (1:9) to 9.0 (9:1). For example, for 

a cost ratio of 0.11, the cost factors for existing and future 

loading reduction were 1 and 9, meaning that it was 9 times 

more costly to reduce nitrate loading from future homes. To 

define the extremes, two additional solutions were computed 

in which the unit cost was set to zero, first for existing, and 

then for future homes. The only constraints were ground-water 

concentration constraints of 7 and 3 mg N/L for the shallow 

and deep sites. The loading from existing and future homes for 

each cost ratio are shown in figure 24. 

The highest optimal total loading (84 kg/d) occurs when 

the cost ratio (existing:future) is 1.0 (costs are equal). As the 

ratio increases or decreases to favor reductions in future or 

existing homes, respectively, the optimal balance in loading is 

affected because the objective function (to minimize cost) is 

most efficiently reduced by eliminating loading from homes 

that have the lowest cost per unit reduction. As the cost ratio 

increases or decreases from 1.0, the total loading that can 

be maintained decreases because cost variables now act as 

weighting factors that partially determine which management 

areas and home types will have reduced loading. It is not likely 

that the ratio of unit costs would fall outside the range of 0.43 

to 2.33, within which the effect on total loading is relatively 

small. The relative contributions to loading from existing and 

future homes, however, do vary significantly within this range. 

This indicates that wherever possible, true cost factors should 

be incorporated into the NLMM to account for this important 

external variable.

Spatial Distribution of Loading for Optimal 

Solution

The NLMM computes the minimum loading reduction 

needed from existing and future homes in each management 

area to achieve the water-quality goals prescribed by a 

constraint set. To illustrate the variation in optimal loading 

between management areas, optimal reductions in loading 

rates to each management area from both existing and future 

homes were mapped for a specific set of constraints. The 

constraints for this solution included (1) ground-water nitrate 

concentrations could not exceed 7 mg N/L at shallow locations 

and 3 mg N/L at deep locations, and (2) the maximum loading 

reduction possible for either existing or future homes was 96 

percent. The ground-water concentration limits were selected 

to reflect the goals of minimizing loading reductions and 

insuring protection of the part of the aquifer (30 to 50 ft below 

the water table) where most domestic wells obtain water. The 

96 percent limit on loading reduction is consistent with data 

from the La Pine NDP that show that the best-performing 

advanced treatment on-site wastewater systems tested are 

capable of reducing nitrogen to about 2 mg N/L as nitrate in 

effluent leaving the drain field. No constraints were placed on 

discharge of nitrate to streams or minimum loading reduction 

and no differential cost factors were included.

The optimal minimum reductions vary over a broad 

range in the 97 management areas. The breadth of the range 

is a reflection of the complex relations between development 

location and density, hydrogeology, and geochemical 

processes. For the specified constraints, the optimal reduction 

in total loading is 107 kg/d, or 56 percent of the 190 kg/d of 

loading projected for the base scenario. Contributions to buildout loading were reduced from 104 to 46 kg/d for existing 

NITRATE LOADING REDUCTION COST RATIO,

COST OF REDUCING EXISTING HOMES LOADING/

COST OF REDUCING FUTURE HOMES LOADING 

Existing=0

0.110.25

0.430.661.00

1.50

2.334.00

9.00

Future=0

OPTIMAL NITRATE LOADING, IN KILOGRAMS N PER DAY

0

20

40

60

80

100

Higher costs for future homes Higher costs for existing homes

Loading from future homes

Loading from existing homes

EQUAL COST

Optimal total loading

Figure 24. Sensitivity of optimal solutions to relative cost 

difference of nitrate loading reduction for existing and 

future homes in the La Pine, Oregon, study area.
Nitrate Loading Management Model 53

homes and from 86 to 38 kg/d for future homes. The optimal 

reductions in each management area for existing and future 

homes are shown in figures 25 and 26, respectively. 

Overall, the loading from existing homes was reduced 

by the same proportion, 56 percent, as the loading from 

future homes because there was roughly the same number of 

existing and future homes (undeveloped lots) in 2000. Locally, 

the reductions in loading from existing and future homes 

for individual management areas often were quite different. 

For example, in management area 33, where most lots had 

been developed by 2000, loading from existing homes and 

future homes would need to be reduced by 81 and 67 percent, 

respectively (figs. 25, 26). Even fewer undeveloped lots 

were in the adjacent management area 34 and the entire 

loading reduction needed to meet concentration constraints 

(93 percent) would need to come from existing homes. 

To further validate the assumption that the simulated 

nitrate concentrations and ground-water discharge nitrate 

loading to streams were linear functions of nitrate loading to 

the aquifer, the base scenario loading (scenario 1) was reduced 

in each management area by the percentages computed using 

the NLMM (figs. 25, 26). This loading was used as input to 

the study-area simulation model and the equilibrium nitrate 

concentrations and discharge loading to streams simulated by 

the study-area model were compared to the values computed 

by the NLMM at the 339 concentration and 14 discharge 

loading constraint locations and reaches. The simulated 

concentrations and discharge loading were equal to the values 

computed by the NLMM using the response coefficients 

at all constraint locations and reaches. The detailed spatial 

distribution of loading reduction illustrated in this example 

could be used by decision makers to delineate ground-water 

protection zones and set performance standards for on-site 

wastewater systems to achieve the needed loading reductions.

Comparison of Scenario Simulations and Optimal 

Solution

Simulated ground-water nitrate concentrations were used 

to compare the results of the scenario simulations with the 

results of an optimal solution computed using the NLMM. 

For the comparison, concentration values were compiled at 

the 339 NLMM constraint locations (fig. 19) for four of the 

scenarios simulated using the simulation model in a “trialand-error” mode. The four scenarios (1, 3, 5, and 7 in table 6) 

specified on-site wastewater concentrations at the water table 

of 46, 20, 10, and 2 mg N/L, respectively. For this analysis, a 

fifth scenario using a wastewater concentration of 30 mg N/L 

also was simulated. Scenario 1 was the base scenario where 

conventional on-site wastewater systems yield 46 mg N/L 

nitrate at the water table. The other four scenarios are based 

on the assumption that all existing and future homes would 

install advanced treatment on-site wastewater systems that 

reduce loading to the same standard. The optimal solution 

used in the comparison is that in which ground-water nitrate 

concentrations could not exceed 7 mg N/L at shallow sites and 

3 mg N/L at deep sites, and the maximum loading reduction 

possible for either existing or future homes was 96 percent. 

The relative effectiveness of the management alternatives 

used in each scenario was evaluated by comparing waterquality improvements (as measured by the percentage of 

constraint locations with nitrate concentrations greater than 

7 mg N/L) with reduced loading for the scenarios and the 

optimal solution (fig. 27). The base scenario (no loading 

reduction) resulted in equilibrium concentrations greater 

than 7 mg N/L for 46 percent of the constraint locations. 

Using advanced treatment systems that can produce nitrate 

concentrations of 30, 20, 10, and 2 mg N/L, reduced loading 

by 35, 57, 78, and 96 percent and reduced the fraction of 

constraint locations where ground-water concentrations 

exceeded 7 mg N/L to 36, 26, 8, and 1 percent, respectively. 

The optimal solution would reduce total loading 

by 56 percent, and yet had only 3 constraint locations 

(<1 percent) with concentrations greater than 7 mg N/L. These 

three locations were affected by a small amount of unmanaged 

loading from dispersed homes not included in the NLMM. 

Loading reduction for the optimal solution was 106 kg/d, 

which was nearly equivalent to the loading reduction in the 

20 mg N/L scenario (3), in which 26 percent of the sites had 

concentrations greater than 7 mg N/L. The 10 mg N/L and 

2 mg N/L scenarios (5 and 7, respectively) are approximately 

equivalent to the optimal solution in reducing nitrate 

concentrations in ground water. Because they apply uniform 

management strategies, however, loading must be reduced by 

40–69 percent more than would be required under the optimal 

solution to achieve similar results.

The scenarios are based on uniform management 

strategies (for example, all new homes have advanced 

treatment on-site wastewater systems and all existing systems 

are replaced), whereas the optimal solution implements 

loading reductions only where reductions are needed to 

meet quality constraints. Decision makers would unlikely 

be able to implement the optimal solution exactly because it 

would be difficult to have variable on-site wastewater system 

performance requirements across management areas as small 

as those defined in the NLMM. More likely, decentralized 

wastewater treatment would be managed over larger areas that 

have similar nitrate loading capacity in the NLMM solution. 

This approach could increase the overall loading reduction 

required to meet water-quality standards, but still be less costly 

than uniform management strategies.
54 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

OR19-0049_fig25

Base 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

Publication projection is 

Lambert Conformal Conic

Standard parallels 42o20' and 44o40', 

central meridian -120o30'

Datum is NAD83

R 9 E R 10 E R 11 E

T 21 S

T 22 S

T 23 S

T 24 S

T 20 S

T 19 S

DESCHUTES

L A K E

KLAMATH

EXPLANATION

Management area with 

number

Optimal reduction from base 

scenario loading—

 In percent

0-10 

11-35

36-57 

58-78 

79-100 

16

Spring

River

Fall

River

Crescent Creek

Creek

Little

Deschutes

DESCHUTES

RIVER

River

Paulina

Long Creek

5 KILOMETERS

0 5 MILES

0

Wickiup

Reservoir

97

97

31

31

Burgess Road

Drive

Finley

Butte

Road

Road

South

Century

Spring

River

Masten Road

SunriverSunriver

La 

Pine

43°55'

50'

45'

121°40' 35' 30' 25' 121°20'

40'

43°35'

Extent of study area model

77

76

88 89

91

87

86

85

84

83 24

24

81

79 80

78 75

73 71

72

68 69 70

67

66

93

56 58

57

61 62

65 54

74

90

64

63

53 52

51 50 32 31

26

35

55 36 34

33

25

48 49

44 46 45

47

43

41

40 39 38 37

42

30

29 27

28

14 94

13 12

19

22 21

95

17 18

16 15

4

11

23

10

9 8

7

1 6

2 3

5

59

60

92

20

Figure 25. Optimal reduction in nitrate loading from existing homes in the La Pine, Oregon, study area.
Nitrate Loading Management Model 55

OR19-0049_fig26

Base 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

Publication projection is 

Lambert Conformal Conic

Standard parallels 42o20' and 44o40', 

central meridian -120o30'

Datum is NAD83

R 9 E R 10 E R 11 E

T 21 S

T 22 S

T 23 S

T 24 S

T 20 S

T 19 S

DESCHUTES

L A K E

KLAMATH

EXPLANATION

Management area with 

number

Optimal reduction from base 

scenario loading— 

In percent

0-10 

11-35 

36-57 

58-78 

79-100 

16

Spring

River

Fall

River

Crescent Creek

Creek

Little

Deschutes

DESCHUTES

RIVER

River

Paulina

Long Creek

5 KILOMETERS

0 5 MILES

0

Wickiup

Reservoir

97

97

31

31

Burgess Road

Drive

Finley

Butte

Road

Road

South

Century

Spring

River

Masten Road

SunriverSunriver

La 

Pine

43°55'

50'

45'

121°40' 35' 30' 25' 121°20'

40'

43°35'

Extent of study area model

77

76

88 89

91

87

86

85

84

83 24

24

81

79 80

78 75

73 71

72

68 69 70

67

66

93

56 58

57

61 62

65 54

74

90

64

63

53 52

51 50 32 31

26

35

55 36 34

33

25

48 49

44 46 45

47

43

41

40 39 38 37

42

30

29 27

28

14 94

13 12

19

22 21

95

17 18

16 15

4

11

23

10

9 8

7

1 6

2 3

5

59

60

92

20

Figure 26. Optimal reduction in nitrate loading from future homes in the La Pine, Oregon, study area.
56 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

Limitations and Appropriate Use of 

Models

The transect and study-area simulation models were 

developed to generate a better understanding of the fate 

and transport of nitrate from on-site wastewater systems at 

multiple scales. The study-area model also may be used to 

help evaluate alternative options for management of nitrate 

loading from on-site wastewater systems. Limitations of 

the modeling software, assumptions made during model 

development, and results of model calibration and sensitivity 

analysis all are factors that constrain the appropriate use of 

these models and highlight potential future improvements.

A simulation model is a means for testing a conceptual 

understanding of a system. Because ground-water flow 

systems are inherently complex, simplifying assumptions must 

be made in developing and applying model codes (Anderson 

and Woessner, 1992). Models solve for average conditions 

(for example, head or nitrate concentrations) within each cell 

using parameters which are interpolated or extrapolated from 

measurements, and (or) estimated during calibration. Practical 

limitations on model size, and hence minimum cell size, are 

imposed by the size and speed of available computers. More 

commonly, however, it is the availability of data to define the 

system that limits the scale and accuracy of the model. In light 

OR19-0049_fig27

46 mg/L

systems

(Scenario 1)

30 mg/L

systems

Optimal

solution

20 mg/L

systems

(Scenario 3)

10 mg/L

systems

(Scenario 5)

2 mg/L

systems

(Scenario 7)

NITRATE LOADING REDUCTION, 

IN KILOGRAMS PER LITER

 0

 50

100

150

200

PERCENTAGE OF SITES GREATER THAN 

7 MILLIGRAMS PER LITER

 0

 10

 20

 30

 40

 50

Figure 27. Comparison of loading and water quality 

between optimal and nonoptimal management scenarios 

for the La Pine, Oregon, study area.

of this, the intent in developing the simulation models was not 

to reproduce every detail of the natural system, but to portray 

its important characteristics in sufficient detail to provide a 

useful tool for testing the conceptual model and evaluating 

alternative management options. 

Simulation Models

The study-area simulation model is a decision-support 

tool for evaluating the effects of wastewater management 

alternatives on ground-water and surface-water quality at the 

neighborhood to watershed scale. The study area and transect 

models are not capable of simulating nitrate concentrations 

at individual wells; however, the transect model (which has 

more than twice the lateral resolution of the study-area model) 

has sufficient detail to approximately simulate the location of 

nitrate plumes. 

The ground-water flow system was assumed to be 

at steady-state, meaning that the velocity and direction of 

ground-water flow did not change with time. Water-level 

variation occurs seasonally and over the long term in response 

to stresses like climatic variation. The variation can change 

the velocity, and possibly direction, of ground-water flow over 

periods ranging from hours to years depending on the cause; a 

change in river stage might affect the system for hours to days 

whereas an extended drought might have effects that last for 

months to years. These changes in the flow system could have 

effects on the fate and transport of nitrate not represented by 

the simulation models. The simulation models are designed 

to evaluate the long-term effects of options for management 

of nitrate loading. The models should not be used to evaluate 

short-term changes without considering the possible effects 

of changes in the ground-water velocity distribution from the 

steady-state conditions represented in the models. 

The location of the boundary between the oxic and 

suboxic parts of the ground-water system was mapped based 

on dissolved oxygen concentrations in 256 wells sampled as 

part of a synoptic sampling of private wells by ODEQ and 

Deschutes County in June 2000. Because denitrification is 

assumed to occur at the oxic-suboxic boundary and nitrate 

concentration below the boundary (in the suboxic zone) is 

specified as zero, simulated nitrate concentrations near and 

below the boundary are sensitive to location. Uncertainty in 

the boundary location will result in uncertainty in simulated 

nitrate concentrations. The distribution of wells used to map 

the boundary was generally good for a study area this size, 

however, the boundary location is less certain in some areas. 

For example, there were fewer wells available to constrain the 

location of the boundary near the margin of the model area and 

near streams. In these areas, model results should be evaluated 

with respect to the effects of uncertainty on simulated nitrate 

concentrations.
Summary and Conclusions 57

The ground-water discharge to evapotranspiration process 

is simulated by the study-area model and accounts for the 

mass of water lost from the system where deep-rooted plants 

extract ground water for transpiration and where ground 

water is shallow enough to be evaporated from bare soil. 

Plants also may take up nutrients dissolved in ground water; 

however, the rate of uptake is highly variable and poorly 

understood in non-agricultural settings. Nutrients and other 

solutes are not removed by evaporation and this process 

results in concentration of solutes in ground water. For this 

study, there was no basis for partitioning the mass of ground 

water discharged by ET into its transpiration and evaporation 

components and it was assumed that no nitrate was taken up 

with the mass of water discharged by ET. This assumption 

may bias simulated nitrate concentrations toward high values 

in areas where ET is a significant part of ground-water 

discharge.

Management Model

Because the NLMM was developed using optimization 

methods with the study-area simulation model, the NLMM is 

subject to the same limitations listed for the study-area model. 

However, additional factors should be considered when using 

the management model that relate to how the management 

problem is formulated. 

The sensitivity analysis of the NLMM presented in this 

report illustrates how closely optimal solutions are tied to the 

definition of the management problem. The NLMM solutions 

were shown to be highly dependent on the value of the 

maximum nitrate concentration constraint and on the number, 

location, and depth of specified constraints. Assignment of the 

constraints is an important part of developing a strategy for 

protecting ground-water resources. 

The management problem for this study was formulated 

with the objective of minimizing the amount of reduction 

in nitrate loading that would be required to meet specified 

water-quality goals within management areas. Management 

area boundaries were defined using the township and section 

lines of the Public Land Survey System (PLSS) and included 

areas ranging from 160 to 640 acres. The management-area 

boundaries do not coincide with the hydrologic, geologic, 

and geochemical boundaries that control the nitrate loading 

capacity of the system. The loading capacity for some 

management areas may be strongly controlled by loading in 

part of the area close to where constraints were specified. 

Large differences in computed optimal reduction requirements 

can occur across management-area boundaries even though 

there may be little difference in lot densities, recharge, depth 

of the suboxic zone, or other factors that affect loading 

capacity. Users of the NLMM need to be cognizant of the 

effects of problem formulation on results and interject 

knowledge of on-the-ground conditions when using model 

results to support management decisions.

Summary and Conclusions

Ground-water is an important resource in the rural 

communities of southern Deschutes and northern Klamath 

County, near La Pine, Oregon. The primary aquifer, and 

only source of drinking water to about 14,000 residents, 

comprises alluvial sand and gravel deposits within 100 ft of 

land surface. Nearly 60 percent of residential lots are less 

than 1 acre and almost all homes use on-site wastewater 

disposal systems. Nitrate concentrations greater than the 

U.S. Environmental Protection Agency drinking water MCL 

of 10 mg/L were discovered in the oldest developed part of 

the area in the late 1970s and elevated concentrations have 

subsequently been detected in more recently developed 

areas. In 2000, nitrate concentrations greater than 4 mg N/L 

were detected in 10 percent of domestic wells sampled by 

Oregon Department of Environmental Quality. Because of 

concern for the vulnerability of the ground-water resource, 

the Oregon Department of Environmental Quality and 

Deschutes County, in cooperation with the U.S. Geological 

Survey, conducted a study to develop a better understanding 

of the hydrologic and chemical processes that affect the 

movement and fate of nitrogen within the shallow aquifers of 

the La Pine region. Simulation models were used to test the 

conceptual understanding of the system and were coupled with 

optimization methods to provide a management model that 

can be used to efficiently evaluate alternative approaches for 

managing nitrate loading from on-site wastewater systems.

The geologic, hydrologic, and geochemical frameworks 

for the conceptual and numerical models were developed 

using several data sources including previous hydrogeologic 

and water-quality studies in the area, an associated, largescale field experiment evaluating advanced treatment on-site 

wastewater systems, literature for similar studies in other 

areas, and extensive field data collection for this study. 

The primary aquifer in the study area is composed of 

complexly interbedded fluvial silt, sand, and gravel deposits. 

A three-dimensional hydrofacies model of the fluvial system 

was created with transition probability geostatistical methods 

using parameters derived from analysis of two-dimensional 

lithologic sections and lithologic data from more than 400 

drillers’ logs. Five hydrofacies were included in the final 

model: clay-silt, sand, gravel, lacustrine clay-silt, and basalt. 

Mean annual ground-water recharge to the alluvial aquifer 

is 3.2 in/yr, primarily from infiltration of precipitation and 

snowmelt. Ground-water discharges to streams, springs, and 

wells, and by evapotranspiration. The water-table generally 

is within 5–20 ft of land surface and varies seasonally over 

a range of a few feet in response to recharge and changing 

stream stage.
58 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

On-site wastewater systems are the only significant 

source of anthropogenic nitrogen to shallow ground water 

in the study area. Low recharge rates and ground-water 

flow velocities have, for now, generally restricted nitrate 

occurrence to discrete plumes within 20–30 ft of the water 

table. Concentrations of nitrate typically are low in deeper, 

older ground water due to the nature and timing of nitrate 

loading and transport, and to loss by denitrification. Ground 

water in the study area evolves from oxic to increasingly 

reduced conditions with increasing depth below the water 

table. Suboxic conditions are achieved in 15–30 years, and the 

transition zone from oxic to suboxic ground water is narrow. 

Nitrate is denitrified near the oxic-suboxic boundary. Nitrate 

loading from residential, commercial, and other sources 

using on-site wastewater systems was estimated for 1960–99 

using county building records, census data, monitoring data 

from field studies of on-site systems, and literature values. 

Adjusted for seasonal residency, residential loading estimates 

ranged from 12 to 14 lbs/yr per home between 1960 and 

1999. During this period total nitrogen loading increased from 

3,900 to 81,000 lb/yr. Nitrogen loading increased by 17,000 

to 98,000 lb/yr between 1999 and 2005. When all approved 

lots are developed (projected to occur in 2019 at current 

building rates), nitrogen loading is estimated to reach nearly 

150,000 lb/yr.

 Three-dimensional numerical simulation models were 

constructed at transect (2.4 mi2) and study-area (247 mi2) 

scales to simulate the fate and transport of nitrate within the 

shallow ground-water system. The transect model was used 

to test conceptual models at the site of a detailed geochemical 

investigation along a 3.5-mi long flow path within the 

study area. The study-area model was constructed at a scale 

appropriate as a planning tool for prediction of average nitrate 

concentrations in neighborhoods and subdivisions. 

Calibration of the models was constrained by data 

that included measured heads, measured and estimated 

ground-water discharge to streams, time-of-travel estimated 

from chlorofluorocarbon age dates, and measured groundwater nitrate concentrations. Eight scenarios representing 

nitrate-loading management strategies were simulated 

for the 140-year period, 2000–2139. A base scenario was 

simulated which assumed existing and future homes would 

continue to use conventional on-site systems and nitrogen 

loading would reach the projected maximum of nearly 

150,000 lb/yr in 2019. Under this scenario, simulated nitrate 

concentrations continue to increase until the rate of nitrate 

loading to the aquifer system is balanced by nitrate losses to 

denitrification and ground-water discharge to the nearstream 

environment of the Deschutes and Little Deschutes Rivers. 

At equilibrium, average nitrate concentrations near the water 

table exceed 10 mg N/L over areas totaling 9,400 acres. 

Other scenarios were simulated that evaluated the effects of 

reduced loading on water quality. Scenarios in which nitrate 

loading was reduced by 15–94 percent overall resulted in 

reductions of 22–99 percent in the area where average nitrate 

concentrations near the water table exceed 10 mg N/L at 

equilibrium. Simulated ground-water ages agree with groundwater age data and show that the system is slow to respond 

to changes in nitrate loading due to low recharge rates and 

ground-water flow velocity. Consequently, reductions in 

nitrate loading will not immediately reduce average nitrate 

concentrations and the average concentration in the aquifer 

will continue to increase for 25–50 years depending on the 

amount and timing of loading reduction. The time required for 

average concentrations to decrease is, in part, also due to the 

assumption that replacement of existing on-site wastewater 

systems would take place over approximately 50 years.

Results of the scenario simulations showed that there is 

variable capacity to receive on-site wastewater system effluent. 

The capacity of an area to receive on-site wastewater system 

effluent is related to many factors, including the density 

of homes, presence of upgradient residential development, 

ground-water recharge rate, ground-water flow velocity, and 

thickness of the oxic part of the aquifer.

The study-area simulation model was used to develop a 

decision-support tool by incorporating optimization methods. 

The resulting model, the Nitrate Loading Management Model 

(NLMM), was formulated to minimize the reduction from 

estimated base scenario loading that would be needed to 

maintain ground-water nitrate concentrations or ground-water 

discharge of nitrate to streams below specified levels. The 

NLMM uses the response matrix approach to find the optimal 

(minimum) loading reductions in each of 97 management 

areas that will meet the specified water-quality constraints. 

The sensitivity of the optimal solutions (loading reductions) 

to water-quality and other constraints was evaluated by 

altering constraint values. The sensitivity of optimal solutions 

to constraints allows decision makers to assess tradeoffs 

between higher levels of water quality protection and the cost 

of reducing nitrate loading. The sensitivity analysis of the 

NLMM showed that optimal (minimum) loading reductions 

were most sensitive to the constraints on ground-water nitrate 

concentration in the shallow part of the oxic ground-water 

system, within 5–10 ft of the water table. 
References Cited 59

Acknowledgments

We gratefully acknowledge the cooperation and 

assistance of the hundreds of residents of southern Deschutes 

and northern Klamath Counties who gave us permission to 

collect geologic and water-level data from their wells, allowed 

us access to their property, and provided other information 

that proved invaluable to this study. This study also benefited 

greatly from the dedicated efforts of many Deschutes County 

Community Development Department employees; special 

thanks go to Barbara Rich, Dan Haldeman, Catherine Morrow, 

Todd Cleveland, Peter Gutowsky, and Tim Berg. Oregon 

Department of Environmental Quality employees Bill Mason, 

Diane Naglee, Robert Bagget, and Dick Nichols, also provided 

valuable assistance during the study. Joe Miller, Nicholas Zerr, 

and Justin Rohrbaugh, all formerly of the U.S. Geological 

Survey, also are acknowledged for their valuable assistance in 

collecting and compiling data for this investigation. 

References Cited

Ahlfeld, D.P., and Mulligan, A.E., 2000, Optimal management 

of flow in groundwater systems: San Diego, Calif., 

Academic Press, Inc., 185 p.

Allen, J.E., 1966, The Cascade Range volcanic-tectonic 

depression of Oregon, in Transactions of the Lunar 

Geological Field Conference, Bend, Oregon, August 1965: 

Oregon Department of Mineral Industries, p. 21-23.

Anderson, C.W., 2000, Framework for regional, coordinated 

monitoring in the middle and upper Deschutes River basins, 

Oregon: U.S. Geological Survey Open-File Report 00-386, 

81 p.

Anderson, M.P., 1987, Field studies in groundwater 

hydrology—A new era: Reviews in Geophysics, v. 25, no. 

2, p. 141-147.

Anderson, M.P., and Woessner, W.W., 1992, Applied groundwater modeling: San Diego, Calif., Academic Press, Inc., 

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.

Boyd, T., 1996, Estimated recharge in the Upper Deschutes 

Basin: Portland, Oreg., Portland State University, M.S. 

thesis, 102 p.

Bower, H., and Rice, R.C., 1976, A slug test for determining 

hydraulic conductivity of unconfined aquifers with 

completely or partially penetrating wells: Water Resources 

Research, v. 12, p. 423-428.

Brutsaert, W.H., 1982, Evaporation into the atmosphere—

Theory, history, and applications: Dordrecht, Holland, D. 

Reidel Co., 299 p.

Carle, S.F., 1999, T-PROGS—Transition probability 

geostatistical software: Davis, Calif., University of 

California, 76 p.

Carle, S.F., and Fogg, G.E., 1996, Transitional probabilitybased indicator geostatistics: Mathematical Geology, v. 28, 

no. 4, p. 453-476.

Carle, S.F., and Fogg, G.E., 1997, Modeling spatial variability 

with one- and multi-dimensional continuous Markov chains: 

Mathematical Geology, v. 29, no. 7, p. 891-917.

Carle, S.F., LaBolle, E.M., Weissman, G.S., VanBrocklin, 

D., and Fogg, G.E., 1998, Geostatistical simulation of 

hydrostratigraphic architecture—A transition probability 

/ Markov approach in Concepts in hydrogeology and 

environmental geology: Society for Sedimentary Geology 

Special Publication, no. 2, p. 147-170.

Century West Engineering, 1982, La Pine aquifer management 

plan: Bend, Oreg., Century West Engineering, 597 p.

Cole, D.L., 2006, Groundwater quality report for the 

Deschutes basin, Oregon: Oregon Department of 

Environmental Quality Laboratory and Water Quality 

Division, 54 p., accessed July 2007 at http://www.deq.state.

or.us/lab/techrpts/groundwater/dbgroundwater/dbgwreport.

pdf

Couch, R., and Foote, R., 1985, The Shukash and La Pine 

Basins—Pleistocene depressions in the Cascades of 

Central Oregon [abs]: Eos Transactions of the American 

Geophysical Union, v. 66, no. 3, p. 24.

Daly, C., and Gibson, W., 2002, 103-year high-resolution 

precipitation climate data set for the conterminous United 

States: Corvallis, Oreg., Spatial Climate Analysis Service, 

accessed May 2002 at ftp://ftp.ncdc.noaa.gov/pub/data/

prism100/.

Fogg, G.E., 1986, Groundwater flow and sand body 

interconnectedness in a thick multiple aquifer system: Water 

Resources Research, v. 22, no. 5, p. 679-694.
60 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

Frimpter, M.H., Donohue, J.J., and Rapacz, M.V., 1990, A 

mass-balance nitrate model for predicting the effects of land 

use on ground-water quality: U.S. Geological Survey OpenFile Report 88-493, 56 p.

Gannett, M.G., and Lite, K.E., 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.

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 WaterResources Investigation Report 00-4162, 77 p.

Gettings, M.W., and Griscom, A., 1988, Gravity model studies 

of Newberry Volcano, Oregon: Journal of Geophysical 

Research, v. 93, p. 10,109-10,118.

Goode, D.J., and Konikow, L.F., 1990, Apparent dispersion in 

transient ground-water flow: Water Resources Research, v. 

26, no. 10, p. 2,339-2,351.

Gorelick, S.M., 1982, A model for managing sources of 

groundwater pollution: Water Resources Research, v. 18, 

no. 4, p. 773-781.

Gorelick, S.M., Freeze, R.A., Donohue, David, and Keely, J.F., 

1993, Ground-water contamination—Optimal capture and 

containment: Boca Raton, Fla. Lewis Publishers, 385 p.

Gorelick, S.M., and Remson, I., 1982, Optimal management of 

groundwater pollutant sources: Water Resources Research, 

v. 18, no. 1, p. 71-76.

Harbaugh, A.W., and McDonald, M.G., 1996, User’s 

documentation for MODFLOW-96, an update to the U.S. 

Geological Survey modular finite-difference ground-water 

flow model: U.S. Geological Survey Open File Report 

96-485, 56 p.

Herrero-Bervera, E., Helsley, C.E., Sarna-Wojcicki, A.M., 

Lajoie, K.R., Meyer, C.E., McWilliams, M.O., Negrini, 

R.M., Turrin, B.D., Donnelly-Nolan, J.M., and Liddicoat, 

J.C., 1994, Age and correlation of a paleomagnetic episode 

in the western United States by 40Ar/39Ar dating and 

tephrochronology—The Jamaica, Blake, or a new polarity 

episode?: Journal of Geophysical Research, v. 99, p. 

24,091-24,103.

Hinkle, S.R., Böhlke, J.K., Duff, J.H., Morgan, D.S., 

and Weick, R.J., 2007a, Aquifer-scale controls on the 

distribution of nitrate and ammonium in ground water near 

La Pine, Oregon, USA: Journal of Hydrology, v. 333, p. 

486-503. 

Hinkle, S.R., Morgan, D.S., Orzol, L.L., and Polette, D.J., 

2007b, Ground water redox zonation near La Pine, 

Oregon—Relation to river position within the aquiferriparian zone continuum: U.S. Geological Survey Scientific 

Investigations Report 2007-5239, 30 p.

Hinkle, S.R., Weick, R.J., Johnson, J.M., Cahill, J.D., Smith, 

S.G., and Rich, B.J., 2005, Organic wastewater compounds, 

pharmaceuticals, and coliphage in ground water receiving 

discharge from onsite wastewater treatment systems near La 

Pine, Oregon—Occurrence and implications for transport: 

U.S. Geological Survey Scientific Investigations Report 

2005-5055, 98 p.

Hutton, S.S., Barber, N.L., Kenny, J.F., Linsey, K.S., Lumia, 

D.S., and Maupin, M.A., 2004, Estimated use of water in 

the United States in 2000: U.S. Geological Survey Circular 

1268, 46 p.

Johnson, N.M., and Dreiss, S.J., 1989, Hydrostratigraphic 

interpretation using indicator geostatistics: Water Resources 

Research v. 25, no. 12, p. 2501-2510.

Jones, Lesley, 2003, Characterization of select water quality 

parameters within the upper Deschutes and Little Deschutes 

study areas: Upper Deschutes Watershed Council Bend, 

Oregon: 2003, 127 p.

LINDO Systems, Inc., 2003, What’s best! User’s manual: 

Chicago, Ill, LINDO Systems Inc., 291 p.

Lite, K.E., and Gannett, M.W., 2002, Geologic framework 

of the regional ground-water flow system in the Upper 

Deschutes Basin, Oregon: U.S. Geological Survey WaterResources Investigations Report 02-4015, 44 p.

Lohman, S.W., 1972, Ground-water hydraulics: U.S. 

Geological Survey Professional Paper 708, 70 p. 

MacLeod, N.S., and Sherrod, D.R., 1992, Reconnaissance 

geologic map of the west half of the Crescent 1 degree by 

2 degrees Quadrangle, central Oregon: U.S. Geological 

Survey Miscellaneous Investigations Map I-2215, scale 

1:250,000.

MacLeod, N.S., Sherrod, D.R., Chitwood, L.A., and Jenson, 

R.A., 1995, Geologic map of the Newberry volcano, 

Deschutes, Klamath, and Lake Counties, Oregon: U.S. 

Geological Survey Miscellaneous Investigations Map 

I-2455, scales 1:62,500 and 1:24,000.

Moosburner, G.J., and Wood, E.F., 1980, Management model 

for controlling nitrate contamination in the New Jersey Pine 

Barrens aquifer: Water Resources Bulletin, v. 16, no. 6, 

p. 971-978.
References Cited 61

Oregon Department of Environmental Quality, 1994, 

Statewide groundwater monitoring program, La Pine area 

groundwater investigation, Deschutes County, Oregon: 

Portland, Oreg, Oregon Department of Environmental 

Quality, v. II, variously paginated.

Oregon Department of Environmental Quality, 2004a, La Pine 

demonstration project fact sheet: Oregon Department of 

Environmental Quality, 2 p.

Oregon Department of Environmental Quality, 2004b, 

Water Quality—La Pine demonstration project: Oregon 

Department of Environmental Quality, accessed May 2004, 

at http://www.deq.state.or.us/wq/lapine/siterptcriteria.asp

Oregon Department of Human Services, 2005, Drinking 

water program—Inventory list for all drinking water 

systems: Oregon Department of Human Services, accessed 

October 7, 2005, at http://170.104.158.45/index.html

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

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

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

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

Robinson, T.W., 1958, Phreatophytes: U.S. Geological Survey 

Water-Supply Paper 1423, 85 p.

Sherrod, D.R., and Pickthorn, L.G., 1989, Some notes on 

Neogene structural evolution of the Cascade range in 

Oregon, in Muffler, L.J.P., Weaver, C.S., and Blackwell, 

D.D., eds., Geology, geophysics, and tectonic setting of the 

Cascade Range: U.S. Geological Survey Open-File Report 

89-178, p. 351-368.

Tchobanoglous, G., and Burton, F.L., 1991, Wastewater 

engineering—Treatment, disposal, and reuse, 3rd ed.: 

Metcalf & Eddy, Inc., New York , McGraw-Hill, Inc.

U.S. Bureau of the Census, 1960, Profile of selected housing 

characteristics—1960, geographic area, Table DP-1, La 

Pine, Oregon: U.S. Bureau of the Census, accessed (2001) 

at http://quickfacts.census.gov/qfd/.

U.S. Bureau of the Census, 1970, Profile of selected housing 

characteristics—1970, geographic area, Table DP-1, La 

Pine, Oregon: U.S. Bureau of the Census, accessed (2001) 

at http://quickfacts.census.gov/qfd/.

U.S. Bureau of the Census, 1980, Profile of selected housing 

characteristics—1980, geographic area, Table DP-1, La 

Pine, Oregon: U.S. Bureau of the Census, accessed (2001) 

at http://quickfacts.census.gov/qfd/.

U.S. Bureau of the Census, 1990, Profile of selected housing 

characteristics—1990, geographic area, Table DP-1, La 

Pine, Oregon: U.S. Bureau of the Census, accessed (2001) 

at http://quickfacts.census.gov/qfd/.

U.S. Bureau of the Census, 2000, Profile of selected housing 

characteristics—2000, geographic area, Table DP-1, La 

Pine, Oregon: U.S. Bureau of the Census, accessed (2001) 

at http://quickfacts.census.gov/qfd/.

U.S. Environmental Protection Agency, 1980, Onsite 

wastewater treatment and disposal systems design manual: 

Cincinnati, Ohio, U.S. Environmental Protection Agency 

report EPA-625/1-80-012, 391 p.

U.S. Environmental Protection Agency, 1994, Ground-water 

modeling compendium, 2nd ed.: U.S. Environmental 

Protection Agency report 500-B-94-004, 175 p.

Vorhis, R.C., 1979, Transmissivity from well data: National 

Water Well Association, Well Log, v. 10, no. 11, p. 50-52.

Weissmann, G.S., Carle, S.F., and Fogg, G.E., 1999, Threedimensional hydrofacies modeling based on soil surveys 

and transitional probability geostatistics: American 

Geophysical Union, Water Resources Research, v. 35, no. 6, 

p. 1761-1770.

Winter, T.C., LaBaugh, J.W., and Rosenberry, D.O., 1988, 

The design and use of a hydraulic potentiomanometer 

for direct measurement of differences in hydraulic head 

between groundwater and surface water: Limnology and 

Oceanography, v. 33, no. 5, p. 1209-1214. 

Zdanowicz, C.M., Zielinski, G.A., and Germani, M.S., 1999, 

Mount Mazama eruption—Calendrical age verified and 

atmospheric impact assessed: Geology, v. 1999, no. 27, 

p. 621-624.

Zheng, C., and Bennett, G.D., 1995, Applied contaminant 

transport modeling—Theory and practice: New York, John 

Wiley and Sons, 656 p.

Zheng, C., and Wang, P.P., 1999, MT3DMS—A modular 

three-dimensional multispecies transport model (release 

DoD_3.50.A) Documentation and User’s Guide: U.S. Army 

Corps of Engineers, Contract Report SERDP-99, 220 p.
62 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

This page intentionally left blank.
Appendix A 63

Appendix A. Vertical Hydraulic Head Gradient Data from Measurements Made 

on the Deschutes and Little Deschutes Rivers, in the La Pine, Oregon, Study 

Area, October 23– November 4, 2000. 

Figure A1. Locations of vertical head gradient measurement sites in the La Pine, Oregon, study area. 

OR19-0049_figA1

Base 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

Publication projection is 

Lambert Conformal Conic

Standard parallels 42o20' and 44o40', 

central meridian -120o30'

Datum is NAD83

R 9 E R 10 E R 11 E

T 21 S

T 22 S

T 23 S

T 24 S

T 20 S

T 19 S

DESCHUTES

L A K E

KLAMATH

Spring

River

Fall

River

Crescent Creek

Creek

Little

Deschutes

DESCHUTES

RIVER

River

Paulina

Long Creek

5 KILOMETERS

0 5 MILES

0

Wickiup

Reservoir

97

97

31

31

Burgess Road

Drive

Finley

Butte

Road

Road

South

Century

Spring

River

Masten Road

SunriverSunriver

La 

Pine

43°55'

50'

45'

121°40' 35' 30' 25' 121°20'

40'

43°35'

Extent of study-area model

EXPLANATION

Gradient measurement site 

with number

L02-0N

D192-5N

L02-0N

L05-4F

L07-8N

L09-8N

L12-5N

L14-8F

L16-7N

L19-1N

L21-2N

L23-0N

L24-9F

L26-5N

L28-1N

L29-5N

L32-1N

L34-2N

L37-6N

L38-6N

L40-8N

L42-6N

D194-7N

D196-8N

D198-8N

D200-6N

D201-9F

D203-8N

D207-0N

D209-4N

D211-7N

D212-7N

D218-4N
64 Evaluation of Approaches for Managing Nitate Loading from On-Site Wastewater Systems near La Pine, Oregon

OR19-0049_TableA1

Table A1. Vertical hydraulic head gradient data from measurements made on the Deschutes and Little Deschutes Rivers, in the La Pine, Oregon, study area, October 23–November 4, 2000.[Nutrients: Yes, ground-water samples collected and anlayzed for nitrogen by Oregon Department of Environmental Quality, using methods described in Hinkle and others, 2007a. Depth: average stream depth at time of measurement. Width, stream width at time of measurement. HW_depth: average depth at high water mark. HW_width: width at high water. SURFACE_LITH: streambed lithology at surface (see lithology codes). UNDERLYING_LITH: streambed lithology 0.5–1.0 foot below surface. MEAN_GRADIENT: vertical hydraulic gradient between stream and sediments 1.0 foot below steambed surface. Gradients are in units of foot per foot, positive upward. Lithology codes: CL, clay; CS, claystone; GC, gravel, coarse; GF, gravel, fine; GR, gravel, medium; ML, silt; SA, sand, medium; SC, sand, coarse; SF, sand, fine; TU, tuff. Abbreviation: ft, foot. Symbol: >, greater than]

Site 

identification 

No.

Latitude Longitude Date Nutrients Depth (ft) Width (ft) HW_depth (ft) HW_width (ft) SURFACE_LITH UNDERLYING_LITH MEAN_GRADIENT Remarks

Deschutes River

D192-5N 43.854.04 -121.453.41 10-26-00 3.0 120 5.0 150 GFSA 0.020 

D194-7N 43.838.06 -121.462.77 10-26-00 4.0 120 6.0 150 SASF SCGF .040 

D196-8N 43.823.58 -121.471.84 10-26-00 4.0 100 7.0 140 SAGR SA .045 

D198-8N 43.817.41 -121.491.07 10-26-00 2.0 195 4.0 235 SCSA SCGF .025 

D200-6N 43.809.13 -121.497.73 11-04-00 3.0 100 5.0 120 GFGR .050 Hard (tuff?) below 1 ft.

D201-9F 43.796.98 -121.500.63 11-04-00 3.0 60 5.0 100 SCSA GF .010 

D203-8N 43.791.82 -121.507.65 11-04-00 4.0 80 6.0 100 SA .025 

D207-0N 43.772.66 -121.518.44 11-04-00 2.0 110 5.0 150 GRSA GCML .060 

D209-4N 43.771.66 -121.538.71 10-26-00 2.5 90 4.0 130 GCTU TU .010 Gravel on bedrock.

D211-7N 43.763.98 -121.567.13 10-25-00 2.0 105 3.5 125 GCTU GRSA .031 Gravel on bedrock (0.8 ft 

penetration), iron rich.

D212-7N 43.760.80 -121.580.55 10-25-00 2.0 106 4.0 120 GCTU TU .035 Gravel/cobbles on bedrock.

D218-4N 43.732.73 -121.622.24 10-24-00 2.5 128 4.5 140 GCTU TU .045 gravel/cobbles on bedrock.

Little Deschutes River

L02-0N 43.841.25 -121.445.92 10-27-00 3.0 46 6.0 60 SASC GR .033 Blocks of tuff on rt bank.

L05-4F 43.820.99 -121.451.21 10-27-00 3.0 42 5.0 50 SC SAML .028 

L07-8N 43.801.55 -121.454.10 10-27-00 Yes 2.5 42 4.5 46 SC MLGR .010 

L09-8N 43.791.56 -121.459.23 10-27-00 Yes 2.0 44 4.5 52 SC SFML .020 Depths 1, 2, 3 ft.

L12-5N 43.778.04 -121.471.14 10-27-00 Yes 1.5 47 3.5 57 SCML MLSF .035 Depths 1, 2, 3 ft.

L14-8F 43.762.16 -121.480.19 10-27-00 2.0 35 4.0 55 GFCS SAML .050 Three measurements at 1-ft depth.

L16-7N 43.749.33 -121.488.87 10-26-00 Yes 2.0 31 3.0 50 SAGR SCML .015 

L19-1N 43.737.57 -121.490.42 10-26-00 Yes 2.0 40 3.5 50 SCGF SFML .015 Measurements at 1, 2 ft.

L21-2N 43.725.54 -121.491.91 10-26-00 Yes 2.0 40 3.0 45 SC CSML .010 No yield at 1 ft, gas; this is 2 ft 

measurement.

L23-0N 43.714.71 -121.494.13 10-26-00 2.0 32 4.0 36 SASC SFML .020 

L24-9F 43.703.25 -121.499.02 10-23-00 Yes 1.5 57 3.5 77 GFSA .017 15 measurements, 3 per depth 

(1–3 ft).

L26-5N 43.691.59 -121.503.07 10-25-00 2.0 37 4.0 47 SCGF MLSA .010 Measurements at 1, 2 ft.

L28-1N 43.685.07 -121.513.44 10-25-00 2.5 36 4.0 45 SCGF MLSF .000 Measurements at 1, 2 ft.

L29-5N 43.680.99 -121.520.43 10-25-00 2.0 35 3.5 40 SCGF CL -.003 Clay at >1 ft, no yield at 2 ft.

L32-1N 43.674.35 -121.542.20 10-25-00 2.0 40 3.5 46 SCGR ML .010 

L34-2N 43.669.07 -121.563.59 10-25-00 2.0 28 3.0 30 GF CL .025 Clay >1 ft, no yield at 2 ft.

L37-6N 43.650.22 -121.591.21 10-24-00 1.5 47 3.5 52 GR CL .075 Clay >2 ft, no yield at 3 ft.

L38-6N 43.640.24 -121.590.31 10-24-00 2.0 38 3.5 40 GRCS ML .025 Silt >1 ft, no yield at 2 ft.

L40-8N 43.628.71 -121.585.89 10-24-00 2.0 50 4.0 55 GFSA SAML .010 Silt >1 ft, no yield at 2 ft.

L42-6N 43.617.56 -121.600.17 10-24-00 1.5 40 3.0 50 SCGR ML .010 Silt >1ft, no yield at 2 ft.
Manuscript approved for publication, September 29, 2007

Prepared by the USGS Publishing Network, 

Publishing Service Center, Tacoma, Washington

Bob Crist

Bill Gibbs

Debra Grillo

Jackie Olson

Bobbie Jo Richey

Sharon Wahlstrom

For more information concerning the research in this report, contact 

Director, Oregon Water Science Center

U.S. Geological Survey, 2130 SW 5th Avenue

Portland, Oregon 97201

http://or.water.usgs.gov
Morgan and others—Evaluation of Approaches for Managing Nitrate Loading near La Pine, Oregon—Scientific Investigations Report 2007-5237